距离测度学习

上个学期快结尾时,我和小南决定找点东西做做。左看看右看看,觅了个距离测度学习来研究研究,这个在matlab中已经有了相应的包,只是功能并不是全,于是我们觉得把这个还算有点意义的东西写成R包,算作一个短期的学习。

不过人总是有惰性的,上个学期末我写了个距离测度学习的总结草稿,本打算寒假时候开始R包工作,但是谁知寒假毫无动力,堕落的只剩下吃喝睡觉。伴随着这个新学期的开始,毕业论文现在还没有眉目,距离测度学习这个东西小南却已经在github上开了个头,哪有半途而废之理。于是经过一天的梳理,我把这个总结的文档又重新完善了一下,虽然没有任何实例,但是对写R包还是很有帮助。现在赶紧挂出来,也算给自己一个督促吧。毕竟还有个毕业论文等着我呢。

距离测度学习文档总结

Mandelbrot Set

人懒的时候容易被一些事情所烦,需要找些事情放松放松。看看图形是个不错的选择,比如经典的Mandelbrot集。虽然R画的没有绚烂,但是还能看。下面是搜集的稍作修改的一段代码和图形。

1.ggplot2的图形,循环迭代较多,运行确实很慢。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
library(ggplot2)
max_iter = 25
cl = colours()
step = seq(-2, 0.8, by = 0.005)
points = array(0, dim = c(length(step)^2, 3))
t = 0
for (a in step) {
    for (b in step + 0.6) {
        x = 0
        y = 0
        n = 0
        dist = 0
        while (n < max_iter & dist < 4) {
            n = n + 1
            newx = a + x^2 - y^2
            newy = b + 2 * x * y
            dist = newx^2 + newy^2
            x = newx
            y = newy
        }
 
        if (dist < 4) {
            color = 24
        }
        else {
            color = n * floor(length(cl)/max_iter)
        }
 
        t = t + 1
        points[t, ] = c(a, b, color)
    }
}
df = as.data.frame(points)
ggplot(data = df, aes(V1, V2, color = cl[V3])) + geom_point(size = 0.5) +
    opts(panel.background = theme_blank(), panel.grid.major = theme_blank(),
        panel.grid.minor = theme_blank(), axis.ticks = theme_blank(),
        axis.text.x = theme_blank(), axis.text.y = theme_blank(),
        axis.title.x = theme_blank(), axis.title.y = theme_blank(),
        legend.position = "none")

渲染的图形还是非常漂亮的,如下(点击图形看大图,否则看不出颜色渲染的效果

2. 动态看图形变化,此转于此处,代码简洁易理解。不过代码做出来的图与原网页展示的图形不同,颜色怎么调都没有他们的漂亮,主要是迭代次数还不够,也不够精细,直接看原链接的吧,搭配颜色是个细致活,伤脑细胞。

爱的声音

今天是母亲节。我是记不住这个节日的,可能本身这个节日在并不算做传统节日的缘故吧。长这么大,我一共陪母亲度过两个母亲节。一个是高三的时候,一个高四的时候。

那两年高考的最后几个月母亲过来陪读,很辛苦,很煎熬。但那是我第一次有了那么多时间与母亲交流,也是从那时候起,我才开始懂事,开始关心母亲,懂得一个小小的表达母亲都会很感动。每天的散步闲谈,让我深深体会到了父母对我的无条件爱,母亲说的有句话我一直刻在在心底“:只要你好,不管有没有出息,过的健康快乐,我和你爸就什么都好,都安心。”每次我一想起这句话,生活学习中的压力顿时烟消云散,母亲诚挚的话语,让我懂得了生活意义。虽然我无数次的质疑生活的意义,但是我从未放弃母亲告诉我的两样东西:爱与理想。

还有个场景我一回想就会热泪盈眶:高考前夕,烈日炎炎下,一群家长焦急等待自己孩子下课,手里都提着自己亲手做的可口饭菜,母亲因为身材娇小,而且怕热,每次都不得不站在人群后面,然后踮着脚,仰着头,远远的眺望,满脸的焦急与期盼。一想起这个场景,我心里都发酸,母亲为我倾注了那么多爱,而我却回报的如此至少,内心充满了愧疚与感动。

一晃,三年就过去了,如此之快。记得第一次我给母亲送了百合花,愿她健康快乐;第二次我送给她一只小乌龟,祝愿她长命百岁。我是一个不善于表达情感的人,满肚子的话,到了嘴边都会被憋回来。那两次送礼物的时候,我只是傻傻的笑,连句祝福不敢说,男孩子的倔强让我说不出那些肉麻的话。至今我也没对母亲说一句“妈,我爱你!”我知道,没有这些话母亲都很感动了,如果说了,怕我们两个都会哭出来,都是天蝎座的人,感情很丰富,很敏感,很脆弱,受不住这样的爱的表达。

三年的每个这个星期,虽然没法再送个我精心准备的礼物,但我都会或早或晚的打电话回去问候下母亲,今天也是一样。不过但是今天看到这个网页时候,我想到了一种表达爱的方式。我要一种比较特别方式送给母亲我一直不敢说的那几个字。不多言说,附图附代码。

1
2
3
4
5
6
7
8
9
10
11
12
library(seewave)
library(tuneR)
woaini = readMP3("我爱你.mp3")
woaini1=resamp(woaini, g = 2000)
muqinjie = readMP3("母亲节快乐.mp3")
muqinjie1 = resamp(muqinjie, g = 2000)
iloveyou = readMP3("i_love_you.mp3")
iloveyou1 = resamp(iloveyou, g = 2000)
op = par(mfrow = c(3,1),mar = c(4.5,4,2,2))
oscillo(woaini1, f = 2000, from = 3, to = 6,cexlab = 0.75,colwave = "red",title = "妈-我爱你")
oscillo(muqinjie1, f = 2000, from = 2, to = 5,cexlab = 0.75,colwave="gold",title = "妈-母亲节快乐")
oscillo(iloveyou1, f = 2000, from = 0, to = 4,cexlab = 0.75,colwave="lightblue",title = "i love you so much")

图解:

  • 最上是我用家乡话(信阳罗山话)说的:“妈,我爱你!”(罗山话中,“你”发音类似于“恩”,故振幅小)
  • 中间一幅依然用家乡话:“妈,母亲节快乐!”(录得不好,貌似有杂音)
  • 最下面的那副是英文:“I love you so much!”

其实从图的走势还是可以看出我羞涩的内心,爱的表达,对我来说,竟然真不是一件容易的事情,但是对许多男生来说,又何尝不是呢?男人有着内心的倔强,都懂得~

提高Monte Carlo模拟定积分计算的精度

前几天在COS主站上看到一篇谈Monte Carlo求定积分的方法的文章。记得去年我学matlab时,求定积分便用到了这两种方法,而在上概率论课时,老师也略微涉及到相关的内容。当时的感觉时,定积分这种一般用确定型算法计算的东西,竟然还可以用随机的方法来实现,确是一种非常巧妙的思想,(Monte, Carlo是摩纳哥的一个赌城,二战时Von, Neumanm一伙人研究原子弹时弄出来的方法。)于是被震撼,也没有再深究其中更详细的内容。

但是昨天一看到这篇文章,便想起寒假时候看的那本《计算方法》中关于Monte, Carlo平均值方法的论述,而最近在上计算方法的课,正纠结与一堆无聊的误差计算中,于是便再次翻开此书看看,竟发现了一些一直没太注意的很重要的章节。

具体如下:

  • 平均值方法是根据弱大数定律证明得来的。构造如此

    I(f)\approx 1/N\sum_{i=1}^{N}f(X_i)=I_N(f)


    均方误差

    E|e_N|^2=E(I_N(f)-I(f))^2=\frac{1}{N}Var(f(X))


    而由Schwartz不等式知

    E|e_N|\leq \sqrt{E|e_N|^2}=\sqrt{Var(f)/N}

    ,从中便可知 MC 法具有半阶收敛性,且收敛阶不依赖问题的维数,这样在计算高维积分的时候,比较优越。但是收敛速度还是很慢的,是一个比较坏的收敛速度。因此在有其他高效率、高精度的确定型算法时,就不要采用MC法,那是在没有更好的算法时更好的抉择,比如超高维积分。
  • 由均方误差可知,要提高 MC 法的精确度,那么久必须得增大N,减小VAR(f),也就是说,要增加样本量,同时改善抽样方法的选择,这时便有很多方法选择,书中提到了四种方法,我就谈谈我的理解和试验。
    1. 分层抽样。核心就是将原来的 [0,1] 区间 M 等分后在各个小区间上分别再以均匀分布抽样 N/M 个点,然后以

      I_N=1/N \sum_{k=1}^{M}\sum_{i=1}^{N_k}f(X_i^{(k)})


      理论证明知这样方差减少,具体证明式子较多就不再写,否则就是一本教材的翻版了。
      模拟\int_{0}^{1}exp(-x^2/2)/\sqrt{2pi}dx,具体实验结果如下 

      1
      2
      3
      4
      5
      6
      7
      
      M<-100
      N<-1000
      A=matrix(as.numeric(gl(M,N/M))-1,ncol=N/M,byrow=T)
      B=replicate(M,runif(N/M))
      x=1/M*t(B)+A/M
      y=exp(-x^2/2)/sqrt(2*pi)
      I=mean(y)

      结果为0.3413287,可知只取了1000个样本点的情况下,精确度已经十分好了,而一般的均值模拟
      达到此精确度则需要10^7个点。

    2. 对偶变量法。
      这是一种特殊的针对于定义域具有一定对称性的函数性质的技巧。基于这样的命题:如果f(x)是单调的,则Cov(f(X),f(1-X))\leq0,这里X\sim U[0,1]。则有

      I_N=1/(2N)\sum_{i=1}^{N}(f(X_i)+f(1-X_i))

      ,
      EI_N=I(f),
      经过证明知 Var(I_N)\leq \frac{1}{2N}Var(f),方差显然减少。 

      但是这种方法只是对于定义域有一定对称性的函数,如果对一种函数先验了解不多,而且对称性不知,则运用会出现失误,所以具有一定局限性。
      模拟\int_{0}^{1}exp(-x^2/2)/\sqrt{2pi}dx的代码如下:

      1
      2
      3
      
      a=runif(10^6)
      y=exp(-a^2/2)/sqrt(2*pi)+exp(-(1-a)^2/2)/sqrt(2*pi)
      b=sum(y)/(2*10^6)

      结果为 0.3413487,精确度也得到了提高。

    3. 重要性抽样:核心就是构造一个新的积分

      I(f)=\int_{0}^{1}f(x)dx=\int_{0}^{1}\frac{f(x)}{p(x)}p(x)dx

      其中p(x)[0,1]的概率密度,\int_{0}^{1}p(x)dx=1,p(x)\geq 0。这样得到类似于最初最简单的那种均值法的新的形式积分近似 I(f)\approx 1/N\sum_{i=1}^{1}f(Y_i)/p(Y_i)Y_i为遵循p(y)的独立同分布随机变量(i.i.d),而之所以这样构造,可以根据这幅图来理解,a importance sampling chart\int_{0}^{a}+\int_{b}^{1}f(x)dxI(f) 很少比例(即图像的边缘部分),尽可能的不让随机点去计算这些值,这样的“拟合”积分较好(面积),否则很多的随机点在 [a,b] 外,为计算小比例的东西做贡献,这样必定精度较差,所以我们要期望尽可能的值出现在峰值处。我们称其为“重要性”抽样。理论分析,

      Var_X(f)=\int_{0}^{1}(f-i(f))^{2}dx=\int_{0}^{1}f^{2}dx-I^{2}(f)

      ,

      Var_Y(f/p)=\int_{0}^{1}(f/p)^{2}pdy-I^{2}(f)=\int_{0}^{1}f^2/pdy-I^{2}(f)

      .
      只要选取适当的 p(y) 就可以使方差减少。
      当然不可能是Var(f/p)=0,这在图的模拟便可知,p(y)不可能与f(x)相合的同时,还有\int_{0}^{1}p(x)dx=1。但是p难以选择,需要对问题先验有了一定的了解程度才能适当选取。比如对于正态分布求\int_{0}^{1}exp(-x^2/2)/\sqrt(2pi)dx,我不太好选p,就没有做相关的模拟了。
    4. 控制变量法。
      核心也是构造,数学式子I_N(f)=i/N\sum_{i=1}^{N}(f-g)(X_i)+I(g).
      需要I(g)已知,但是这是同重要性抽样一样,都不知道,也需要对问题先验有了一定的了解程度,才能好的选择,使I(g)已知,并且有Var(f-g)\leq Var(f)。相关模拟我也没有再做

总之,正如书中所言,减小误差都只是具有指导思想的意义,具体问题,往往需要我们对被积函数有先验的了解,才能够减小误差,提高精度。确实如此,我的例子都是比较特殊或者简单的,如果对于复杂的问题,而且对该函数比较怪异少见,那么在以上的方法中寻找什么p(x),都是一项非常艰巨或者不可能完成的任务。因此,我下一结论,奇技淫巧之所以巧无外乎我们有了个千里眼和顺风耳,瞎碰瞎撞弄出一个巧妙的解法那是很难的,大部分都是前人解法突然激发灵感得来(此时此刻竟想起了贝叶斯来...)。所以对于专门解决高维积分强有力算法——Metropolis算法我是敬佩有加,这一算法2000年惊被评为20世纪十大算法之一。后面产生的什么模拟退火,拟Monte Carlo都与他有莫大的渊源。而伟大的东西终究充满了高深而深刻的东西,我目前楞死没有看太明白,已经被符号给堆死了。但是我想过段时间在查阅相关资料应该可以弄明白吧。

最后再扯句从中学到了的随机向量的构造——正态分布随机变量——Box-Muller法
以前学的是用n个均匀分布随机数由中心极限定律构造,但是那个要得到较好的精度需要N很大,故需要精确的还是用上面的那个方法更比较有保证,咱看公式也心安嘛。