统计学之路

正则简单笔记

发表于： 2010年05月28日分类：统计之路 2 Comments 标签: regular expression, 正则表达式, 笔记.

1. 元字符 (Metacharacter)

匹配字符串的开始(位置): ^

匹配字符串的开始(不受处理多行选项影响): \A

匹配字符串的结束(位置): $

匹配字符串的结束(不受处理多行选项影响): \z

匹配字符串的开始/结束(位置): \b

匹配字符串的结束/行尾(不受处理多行选项影响): \Z

匹配数字/字母/下划线/汉字: \w

匹配(阿拉伯)数字: \d

匹配空白符(空格/全角空格/制表符/换行符): \s

匹配横向制表符: \t

匹配竖向制表符: \v

匹配回车: \r

匹配换行符: \n

匹配换页符: \f

匹配任意字符(除换行符\n): .

前述表达式重复任意次: *

2. 元字符转义 (查找元字符本身)

\. \* \\ \^ \$ 等.

3. 字符类 (自定义的元字符): [ ]

e.g. [13579] 匹配1或3或5或7或9中的任意一个数字.

4. 限定符 (规定表达式重复次数)

重复至少0次: *

重复至少1次: +

重复0次或1次: ?

重复n次: {n}

重复n次或更多次: {n,}

重复m到n次: {m,n}

e.g. \d{2,9} 匹配2-9位阿拉伯数字.

5. 分枝: |

p.s. 分枝条件按照从左到右的顺序测试. 如果已满足了某个分枝, 则匹配结束. 故使用分枝条件时, 一定要反复斟酌各个条件的书写顺序.

6. 反义

对于元字符, 把元字符里的字母小写换为大写即表示反义. 对于字符类, [^13579] 匹配除了1或3或5或7或9以外的任意字符. 一个不错的例子:

<a[^>]+> 匹配用尖括号括起来的以a开头的字符串

7. 分组: ( )

和编程写公式时多加括号防止写错有点类似. 不过这里的现实意义强得多, 可以用来写出比较复杂的表达式.

8. 贪婪/懒惰

在正则表达式中, [最先开始的匹配] 拥有最高优先权. 其次才是默认的贪婪匹配原则. 与之对应的懒惰匹配原则:

重复0次或1次: ??

重复至少1次: +?

重复m到n次: {m,n}?

重复n次以上: {n,}?

重复任意次: *?

9. 零宽断言

(1) 零宽度正预测先行断言: (?=exp)

断言此位置之后能匹配表达式exp

e.g. \b\w+(?=ly\b) 匹配以ly结尾的字符串的前面部分.

(2) 零宽度正回顾后发断言: (?<=exp)

断言此位置之前能匹配表达式exp

e.g. (?<=\bpre)\w+\b 匹配以pre开头的字符串的后面部分.

(3) 零宽度负预测先行断言: (?!exp)

断言此位置之后不能匹配表达式exp

e.g. \b((?!exe)\w)+\b 匹配不包含连续字符串exe的字符串.

(4) 零宽度负回顾后发断言: (?<!exp)

断言此位置之前不能匹配表达式exp

e.g. (?<!010)\d{8} 匹配前面不是010的8位数字.

10. 注释

(?#注释内容)

复化Gauss-Legendre求积

发表于： 2010年05月28日分类：数字之美 0条评论标签: 勒让德, 复化, 数值积分, 高斯.

高斯型求积公式是数值积分中一个比较成熟的想法. 速度快, 精度高, 公式优雅.

一般Gauss-Legendre求积公式

对于一般的Gauss型求积公式

$\int_a^b f(x) \rho (x)\,dx \approx \sum_{k=0}^{n} A_k f(x_k)$
取权函数 $\rho (x)=1$ , 积分区间设定为 $\left[-1, 1\right]$ , 则得到一般的Gauss-Legendre求积公式
$\int_{-1}^{1} f(x) dx \approx \sum_{k=0}^{n} A_k f(x_k).$
取Legendre多项式的零点作为Gauss点, 通过计算得到Gauss点个数 $n=2$ 和 $n=3$ 时的求积公式
$\int_{-1}^{1} f(x)\,dx \approx f(-\frac{1}{\sqrt{3}})+f(\frac{1}{\sqrt{3}})$
以及
$\int_{-1}^{1} f(x)\,dx \approx \frac{5}{9} f(-\frac{\sqrt{15}}{5})+\frac{8}{9}f(0)+\frac{5}{9}f(\frac{\sqrt{15}}{5}).$

变量替换方法

上述分析中, 积分区间固定为 $\left[-1, 1\right]$ , 实际应用时做变量替换

$x=\frac{a+b}{2}+\frac{b-a}{2}t.$
将被积区间 $\left[a, b\right]$ 化为 $\left[-1, 1\right].$

复化Gauss-Legendre求积公式

将被积区间m等分, 记 $h=\frac{b-a}{m}$ , $x_k=a+kh, k=0, 1, 2 \ldots, m.$ 作变换

$x=\frac{x_k+x_{k-1}}{2}+\frac{h}{2}\,t,$

在每个小区间上应用Gauss-Legendre公式, 累加即得复化Gauss-Legendre求积公式

$\int_a^b f(x)\,dx=\frac{h}{2}\sum\sum A_k f(x_k).$

不妨设

$G(t)=f(\frac{x_k+x_{k+1}}{2}+\frac{h}{2}\,t),$

则有:

Gauss点个数 $n=2$ 时,

$\int_a^b f(x)\,dx\approx\frac{h}{2}\sum G(-\frac{1}{\sqrt{3}})+G(\frac{1}{\sqrt{3}}),$

Gauss点个数 $n=3$ 时,

$\int_a^b f(x)\,dx\approx\frac{h}{2}\sum \frac{5}{9}G(-\frac{\sqrt{15}}{5})+\frac{8}{9}G(0)+\frac{5}{9}G(\frac{\sqrt{15}}{5}).$

总结复化Gauss-Legendre求积过程如下:

1. 分割区间, 记录区间端点值；
2. 通过查表或求解非线性方程组, 在所有小区间上, 将Gauss系数和Gauss点的值代入变量替换后的公式；
3. 将所有区间的结果累加, 即得到整个区间上的积分近似值.

针对Gauss点个数 $n=2$ 和 $n=3$ 的复化Gauss-Legendre求积公式编写的一个简单的MATLAB函数 compgauss() 如下:

^?Download compgauss.m

function [ ] = compgauss(a, b, n)
% Composite Gauss Integration
% Equation Type: n=2, n=3
% Coded by Nan.Xiao 2010-05-25
% Step.1 Divide Interval
% Step.2 Calculate
% Step.3 Sum Results
format long
f = @(x) exp(x).*sin(x);
h=(b-a)/n;
xk=zeros(n+1,1);
xk(1,1)=a;
xk(n+1,1)=b;
fk1=zeros(n,1);
fk2=zeros(n,1);
for i=1:n-1
    xk(i+1,1)=a+h*i;
end
for j=1:n
    fk1(j)=f((xk(j)+xk(j+1))/2+(h/2)*(-1/sqrt(3)))+...
        f((xk(j)+xk(j+1))/2+(h/2)*(1/sqrt(3)));
end
for r=1:n
    fk2(r)=(5/9)*f((xk(r)+xk(r+1))/2+(h/2)*(-sqrt(15)/5))+...
        (8/9)*f((xk(r)+xk(r+1))/2+(h/2)*(0))+...
        (5/9)*f((xk(r)+xk(r+1))/2+(h/2)*(sqrt(15)/5));
end
mysum1=h*sum(fk1)/2;
mysum2=h*sum(fk2)/2;
disp('Result of 2 Nodes:')
disp(mysum1);
disp('Result of 3 Nodes:')
disp(mysum2);
end

总结

1. Gauss求积公式较一般的机械求积公式的进步之处, 在于其针对插值型求积公式进行了改进. Gauss公式将插值节点设为未知, 成功地将代数精度由 $n$ 次提高到 $2n+1$ 次. 计算过程简单, 速度快, 达到要求精度所需步骤较少.

2. (复化) Gauss 求积的一个问题在于, 增加Gauss点个数, 继续求解Gauss系数和Gauss点值时, 需要解一系列非线性方程组. 其解析解较难求得. 而取数值解作为公式中的参数, 不如解析解理想.

时间换空间/空间换时间？

发表于： 2010年05月24日分类：数字之美 0条评论标签: MATLAB, Romberg's Method, 时间换空间, 空间换时间, 算法复杂度, 龙贝格算法.

经常看到网上有人在讨论算法时用时间换空间/空间换时间的说法. 可惜神人们使用的例子都很专业, 再加上一两句晦涩而潇洒的写法, 我等凡人基本上就很难看懂了. 今天在实现数值积分中的龙贝格 (Romberg) 算法时, 算是遇到了一个这方面有点关系的问题, 简单记录一下.

首先上一段800字的原理说明:

……

由于作者心里很清楚读者将会人肉Skip上面一段, 所以这里直接人肉省略其中的797字.

如果执意想看一下, 参考 Wikipedia 好了.

下面是一个很好懂的 MATLAB 实现:

^?Download romberg.m

function [ ] = romberg(aa, bb, kk, iter)
% Romberg Integral Method
% Coded by Nan.Xiao 2010-05-23
% Step 1. Calc TrapeStep
% Step 2. Calc Romberg Method
% Step 3. Extract Diagnal & Output
function [ ] = trapestep(a, b, k)
    format long
%     function y = f(x)
%         y=(2/sqrt(pi))*exp(-x);
%     end
%     Examples Can Also Be:
%     function y = f(x)
%         y=x*sin(x);
%     end
%     Or
    function y = f(x)
        y=x*sqrt(1+x.^2);
    end
k1=2^(k-1);
temp=zeros(k,k1);
mysum=zeros(k,1);
t=zeros(k+1,1);
for i=1:k
    for j=1:(2^(i-1))
    temp(i,j)=f(a+(((2*j-1)*(b-a))/(2^i)));
    end
end
for p=1:k
    mysum(p,1)=sum(temp(p,:));
end
t(1,1)=0.5*(f(a)+f(b));
for q=2:k+1
    t(q,1)=(0.5*t(q-1,1))+(((b-a)*mysum(q-1,1))/(2^(q-1)));
end
end
    format long
trapestep(aa, bb, kk);
diagonal=zeros(kk,1);
tt=zeros(kk,kk);
for v=1:kk
    tt(v,1)=t(v,1);
end
for r=2:kk
for s=2:r
    tt(r,s)=((4^(s-1)*tt(r,s-1))-tt(r-1,s-1))/(4^(s-1)-1);
end
end
save('tt.txt', 'tt', '-ascii', '-double'); % Save the Table
% Extract Table Diagonals
for w=1:kk
    diagonal(w,1)=tt(w,w);
end
disp('Romberg Table Diagonals:')
disp(diagonal);
% Output Final Result
cc=2;
while abs(diagonal(cc,1)-diagonal(cc-1,1))>=iter || abs(diagonal(cc+1,1)-diagonal(cc,1))>=iter
    cc=cc+1;
end
disp('The Final Result is:');
disp(diagonal(cc,1));
end

与网上其他程序不同, 我在这里没有做太多的判断, 将二分的次数k作为一个参数手动设定, 是为了观察计算的具体过程. 同时, 为了方便移植, 仅使用最基本的矩阵操作语句, 没有使用 MATLAB 中的 feval() 函数, 要改变被积函数, 直接修改函数体. 最终目的是, 我要取得计算产生的含有 $T_m^{(k)}$ 的表格, 以方便完成笔头的作业. 懒到一定程度了是不是?
继续阅读’时间换空间/空间换时间？’

小径分岔的花园

发表于： 2010年04月25日分类：R 3 Comments 标签: RgoogleMaps, R语言, 卫星地图, 可视化, 地震, 震源, 震级, 高维数据.

前天晚上发出一篇质量平平的博文后, 太云学长提出一个建议(其实算是中肯的批评): 用点的大小表示震级的大小. 其实以前也不是没想过这个问题. 只是不怎么懂R中的向量操作, 加上自己比较懒, 所以就 … 昨天没去教室学习(极度堕落啊), 正好宅在寝室里没事, 就画出来玩玩.

我们首先绘制图1, 点的半径大小代表震级的高低. 点的半径越大, 震级越高. (数据附文后)

图1 震级的大小

然后绘制图2, 点的透明度高低代表震源的深度, 点的透明度越靠近完全不透明 (完全不透明的红色为#FF0000FF), 表示震源深度越浅.

图2 震源的深度

继续阅读’小径分岔的花园’

青海地震分布与京津地震传言

发表于： 2010年04月23日分类：统计之路 6 Comments 标签: 京津, 传言, 可视化, 因果, 地震, 玉树, 青海.

2010年4月14日, 我国青海省玉树藏族自治州发生了破坏性地震. 用R简单地看了一下2010年2到4月青海地区(含西藏/甘肃等地)的地震分布情况.

图1 2010年2月青海地区地震分布

继续阅读’青海地震分布与京津地震传言’

编译器与Debug的传奇：女牛人Grace Murray Hopper小传

发表于： 2010年04月2日分类：网络3C 5 Comments 标签: compiler, debug, Grace Murray Hopper, 传记, 发明, 编译器, 起源.

这两天读《UNIX痛恨者手册》时看到一句被引用的牛人体语录：“标准的伟大之处在于它可以有很多。—— Grace Murray Hopper.”Google之，意外挖出了一位女性牛人，谨作一文以记之。以下内容系根据已有资料编辑整理而成。

Grace Murray Hopper

引言

1992年1月7日，华盛顿阿灵顿国家公墓，美国海军为在元旦凌晨睡梦中安然去世的退休海军女军官格蕾丝·穆雷·赫柏 (Grace Murray Hopper) 举行了隆重的葬礼。海军仪仗队和众多肃穆的海军官兵按照海军的礼仪向这位令人尊敬的长者作最后的告别。千千万万的美国人则通过电视转播观看了葬礼的实况。四年后的1996年1月6日，美国海军在缅因州的巴斯港 (Bath, Maine) 为它新建造的一艘阿利·伯克级驱逐舰举行了隆重的命名仪式，把它命名为“赫柏号”。这是第二次世界大战以后第一次、整个美国海军历史上第二次以一位女性的名字命名一艘战舰。

启蒙

Grace Hopper (1906–1992), 本姓Murray, Hopper为夫姓。1906年12月9日生于美国纽约一个海军世家，其祖父军衔曾达少将。她的外祖父则是一名高级土木工程师，常常带着她去上班，她也十分高兴地去帮着扶红白相间的测量杆，这培养了她对于几何学和数学的兴趣。Grace的父亲因患动脉硬化导致双腿截肢，长期住院，这使得作为长女的她从小就更加懂事和勤奋。

Grace回忆她小时候最喜欢上的课是数学课，特别是几何课。因为在几何课上，她可以把铅笔盒里所有彩色的笔全部拿出来用。虽然她是个女孩子，可是各种量角器、计算尺她都喜欢拿来玩，研究它们的原理和作用。她还做过一些很像男孩子做的事情：她曾经在六、七岁的时候，把家里所有的钟都拆开，但是没有一个成功装回去，因此还受到了严厉的处罚。