wOOL's Blog

Lethal Sweety . Sunny Yawn . Deicidal Jujube . …

Archive for 十一月 2009

不相关且不独立的2组随机变量的构造 连续&离散

without comments

今天在课本上瞄到这样一句话, 随手写了2组均为(0,1)分布的随机变量, 发现如果不相关的话一定就是相互独立的, 于是我就纠结了, 后来才知道, 对于只有2个取值的2组随机变量, 不相关和独立是等价的, 但是通常情况下: 独立=>不相关 但是 不相关 =/=> 独立, 下面是一些例子

连续型:
设X为均值为0的服从正态分布的随机变量
设Y=X^2
由Cov(X,Y)=E(XY)-E(X)E(Y)=E(X^3)-0=E(X^3)=0 可知 X, Y不相关
但是
f(y|x) =/= f(y)
所以X, Y也不独立
类似的例子WIKI上也有: Normally distributed and uncorrelated does not imply independent

离散型:
这个就有些复杂了, 刚才说过, 如果随机变量必须有3种或者以上的取值才能构造出来的, 下面是一个例子

\begin{array}{ccccc}\\ X,Y & 2 & 1 & -1\\5 & 0.085 & 0.0475 & 0.0675\\ 2 & 0.090 & 0.0325 & 0.0775\\4 & 0.225 & 0.1200 & 0.2550\end{array}

E(X)=3.8
E(Y)=0.6
E(X)E(Y)=2.28
E(XY)=2.28
所以 不相关
P(X=5|Y=2)=0.0845/0.4=0.2124 =/= P(X=5)=0.2
所以 不独立
有一篇文章专门介绍这种离散型随机变量的构造, 文章在这里不相关且不独立二维离散型随机向量的构造方法
文中附带了一个生成程序, 但是是一个试卷生成软件用的代码, 也不复杂, 下面是我改写的Python版

  1. from numpy import array,dot,min,mat
  2. import random
  3.  
  4. def gcd(a,b):
  5.     while b > 0: a,b = b, a%b
  6.     return a
  7.    
  8. x2=random.randint(1,3)
  9. x1=random.randint(3,5)
  10. x3=random.randint(3,5)
  11. q=gcd(gcd(x3-x2,x1-x3),gcd(x1-x3,x2-x1))
  12. y1=random.randint(1,3)
  13. y2=random.randint(-1,1)
  14. y3=random.randint(-1,1)
  15. r=gcd(gcd(y3-y2,y1-y3),gcd(y1-y3,y2-y1))
  16. t=random.randint(1,2)
  17. s=random.randint(1,2)
  18. a1=random.randint(1,2)
  19. a2=random.randint(1,2)
  20. a3=5-a1-a2
  21. xi=0.2*array([[a1,a2,a3]])
  22. xi=mat(xi)
  23. b1=random.randint(1 ,2)
  24. b2=random.randint(1 ,2)
  25. b3=5-b1-b2
  26. eta=0.2*array([[b1],[b2],[b3]])
  27. eta=mat(eta)
  28. X=array([x1,x2,x3])
  29. Y=array([y1,y2,y3])
  30. H=array([[(y3-y2)/r,(y1-y3)/r,(y2-y1)/r],\
  31.         [t,0.0,-t]])
  32. G=array([[s,(x3-x2)/q],\
  33.         [-s,(x1-x3)/q],\
  34.         [0.0,(x2-x1)/q]])
  35. GH=dot(G,H)
  36. etaxi=dot(eta,xi)
  37. c=min(GH)
  38. d=min(etaxi)
  39. b=-d/c
  40. A=b*GH+etaxi
  41. print A
  42. print x1,x2,x3
  43. print y1,y2,y3
  44. raw_input("\nXD")
  45.  

Written by wOOL

19/11/2009 at 22:13

Posted in Calvados

所学过的一些数值方法及其Python实现

without comments

回想起来, 我在初中的时候, 由于懒, 计算爱用计算器搞定, 但是计算器给的是小数呀, 为了把答案还原成分数(这样看上去才想自己算的), 还专门去学了一下如何把一个无限循环小数变成分数, 当然那样做还是很麻烦, 于是想到了把得到的结果再用计算器依次乘从2, 3…这样的素数(虽然当时不知道原因, 但是显然这个方法十分有效), 以至于后来10以内的最简分数的小数值都记住的差不多了, 大概对于我这样滴人, 如果一个方程的解是(1,2,3), 那么断然不会去接受(1.000001,1.999998,3.000003)这样的答案吧, 不过数值分析就是给出我们这样的答案的, 而悲剧在于数值分析是我这学期的一门课… 我当然承认这门课很有用, 抛弃之前对于小数的厌恶来说, 甚至很好玩, 尤其是Why it works需要的推理都蛮有趣, 而且我淘到一本好书, 基本上常见的数值分析算法都给出了Python代码, 于是我就有些热情了

0. 我接触到的第一个数值分析方法: 牛顿法求平方根
我显然不是这个学期学的这个东西, 看到它是因为一篇文章《John Carmark密码：0×5f3759df》
1597465647 1597463007

1. 解线形方程组(System of linear equations)
1.0 Gaussian elimination
高斯消元法, 这个有点唬人的名字其实就是我们初中就学过的加减消元法来着, 代码实现非常简单:

  1. from numpy import array,dot
  2. def gaussElimin(a,b):
  3.     n = len(b)
  4.     for k in range(0,n-1):
  5.         for i in range(k+1,n):
  6.             if a[i,k] != 0.0:
  7.                 lam = a [i,k]/a[k,k]
  8.                 a[i,k+1:n] = a[i,k+1:n] – lam*a[k,k+1:n]
  9.                 b[i] = b[i] – lam*b[k]
  10.     for k in range(n-1,-1,-1):
  11.         b[k] = (b[k] – dot(a[k,k+1:n],b[k+1:n]))/a[k,k]
  12.     return b
  13.    
  14. A=array([[0.003,59.14],[5.291,-6.13]])
  15. B=array([[59.17],[46.78]])
  16. print gaussElimin(A,B)

这里有一个补充的知识, 就是如何在比较少的有效数字下, 取得尽可能精确的结果, 比如如果我们只取4位有效数字, 上面那个方程的结果就会变成(-10.01, 1.001), 正确的是(10,1), 误差还是很大的. 但是如果我们把2个方程的顺序交换一下, 依然只取4位有效数字, 得到的结果就是(10.00,1.00). 原因非常简单, 舍弃一个比较小的数字的后面几位, 造成的影响, 必然比舍弃比较大的数字的后面几位小. 比如上面的例子中, 如果不调换顺序, 我们首先计算的是5.291/0.003=1763.66666… 由于只能取4位有效数字, 这个结果便成为了1764, 误差0.3333…. 如果我们调换一下顺序, 计算0.003/5.291=0.000567000567000567…. 取4位有效数字, 0.0005670, 误差只有0.000000000567000567….远小于0.33333…. 为了使方程变成这样, 我们要做的就保证: 在进行第i次高斯消元时, 第i个式子中的x_i的系数是最大的, 实现代码如下

  1. from numpy import array,dot,max,argmax,zeros,float64
  2.  
  3. def swapRows(v,i,j):
  4.     temp = v[i].copy()
  5.     v[i] = v[j]
  6.     v[j] = temp
  7.  
  8. def gaussPivot(a,b,tol=1.0e-9):
  9.     n = len(b)
  10.     s = zeros((n),float64)
  11.     for i in range(n):
  12.         s[i] = max(abs(a[i,:]))
  13.     for k in range(0,n-1):
  14.         p = int(argmax(abs(a[k:n,k])/s[k:n])) + k
  15.         if abs(a[p,k]) < tol:
  16.             print "Matrix is singular"
  17.         if p!=k:
  18.             swapRows(b,k,p)
  19.             swapRows(s,k,p)
  20.             swapRows(a,k,p)
  21.         for i in range(k+1,n):
  22.             if a[i,k] != 0.0:
  23.                 lam = a [i,k]/a[k,k]
  24.                 a[i,k+1:n] = a[i,k+1:n] – lam*a[k,k+1:n]
  25.                 b[i] = b[i] – lam*b[k]
  26.     for k in range(n-1,-1,-1):
  27.         b[k] = (b[k] – dot(a[k,k+1:n],b[k+1:n]))/a[k,k]
  28.     return b
  29.        
  30. A=array([[0.003,59.14],[5.291,-6.13]])
  31. B=array([[59.17],[46.78]])
  32. print gaussPivot(A,B)

1.1 LU decomposition
LU分解法, 其实这个不是一种解方程得方法, 而是一种对系数矩阵进行分解的方法, A=LU, A系数矩阵, L是一个下三角矩阵, U是一个上三角矩阵, Ax=b LUx=b, 令Ux=y 则 Ly=b. 之所以要做这样很蛋疼得分解, 是因为如果我们有很多方程组, 他们的系数矩阵, 也就是A都是相同的, 但是b不同, 那么适用这种方法, 我们只需要对于A分解一次, 通过Ly=b计算出y是很方便的(因为是L是三角矩阵, 相当于已经作好高斯消元的情况), 再通过Ux=y也可以很方便的计算出x, 计算L, U是很恶心的, 我笔算算错N次, 不过交给计算机实现起来就很XX了, Python代码:

  1. from numpy import array,dot,zeros,eye
  2. def LUdecomp(a):
  3.     n = len(a)
  4.     for k in range(0,n-1):
  5.         for i in range(k+1,n):
  6.             if a[i,k] != 0.0:
  7.                 lam = a [i,k]/a[k,k]
  8.                 a[i,k+1:n] = a[i,k+1:n] – lam*a[k,k+1:n]
  9.                 a[i,k] = lam
  10.     return a
  11. def LUshow(a):
  12.     n=len(a)
  13.     U=zeros((n,n))
  14.     for i in range(n):
  15.         for j in range(n):
  16.             U[i][j]=a[i][j]
  17.             if j<i: U[i][j]=0
  18.     L=a-U+eye(n,n,0,dtype=int)
  19.     return L,U
  20. def LUsolve(a,b):
  21.     n = len(a)
  22.     for k in range(1,n):
  23.         b[k] = b[k] – dot(a[k,0:k],b[0:k])
  24.     for k in range(n-1,-1,-1):
  25.         b[k] = (b[k] – dot(a[k,k+1:n],b[k+1:n]))/a[k,k]
  26.     return b
  27. A=array([
  28.     [1,1,0,3],
  29.     [2,1,0-1,1],
  30.     [3,-1,-1,2],
  31.     [-1,2,3,-1]
  32.     ])
  33. B=array([4,1,-3,4])
  34. LU=LUdecomp(A)
  35. X=LUsolve(A,B)
  36. print "L=%s\n U=%s" % LUshow(LU)
  37. print X

1.2 Jacobi method
1.3 Gauss–Seidel method
1.4 Successive over-relaxation

2. 矩阵运算
2.1 求特征值 Power Method
待补充

  1. from numpy import array,dot
  2. from math import sqrt
  3.  
  4. s = array([[3.0, 0.0, 1.0], \
  5.            [2.0, 2.0, 2.0], \
  6.            [4.0, 2.0, 5.0]])
  7.            
  8. v = array([1.0, 1.0, 1.0])
  9. for i in range(100):
  10.     vOld = v.copy()
  11.     z = dot(s,v)
  12.     zMag = sqrt(dot(z,z))
  13.     v = z/zMag
  14.     if dot(vOld,v) < 0.0:
  15.         sign = -1.0
  16.         v=-v
  17.     else: sign = 1.0
  18.     if sqrt(dot(vOld – v,vOld – v)) < 1.0e-6: break
  19. lam = sign*zMag
  20. print "Number of iterations = %d" %i
  21. print "Eigenvalue = %f" % lam
  22. raw_input("\nPrint press return to exit")

(下面这个是课上内容对应的版本)

  1. from numpy import array,dot
  2. from math import sqrt
  3.  
  4. s = array([[3.0, 0.0, 1.0], \
  5.            [2.0, 2.0, 2.0], \
  6.            [4.0, 2.0, 5.0]])
  7.  
  8. v = array([1.0, 1.0, 1.0])
  9. for i in range(100):
  10.     vOld = v.copy()
  11.     z = dot(s,v)
  12.     zMag = max(z)
  13.     v = z/zMag
  14.     if sqrt(dot(vOld – v,vOld – v)) < 1.0e-6: break
  15. lam = zMag
  16. print "Number of iterations = %d" %i
  17. print "Eigenvalue = %f" % lam
  18. raw_input("\nPrint press return to exit")

Written by wOOL

18/11/2009 at 21:20

Posted in Calvados

行列式 矩阵 特征向量 特征值 的 几何解释

with one comment

总结一些东西, 主要是在二维和三维空间中对行列式 矩阵 特征向量 特征值这4个概念进行直观的阐述

行列式: 一个面积/体积/超体积…
阐述oo
矩阵: 一个坐标系/一种线性算子
阐述xx
特征向量: 这样的向量在被某矩阵作用后, 只有其长度变换了, 方向不会发生逆向之外的改变
阐述xx
特征值: 上述特征向量长度变换的尺度, 若其为负值, 则说明方向改变了
阐述oo

To be continued

很早之前听说过一个试验, “拿出一张纸和一支笔，在纸上画一个圆形。注意，慢慢画，别一口气画完整个圆，最后要留下一个小小的缺口。放下笔，看它一眼，过一会儿再看它一眼。你是不是有一种冲动：总是不自觉地想把这个缺口补上，让它变成一个完整的圆？如果不完成这个圆，是不是感到有点儿不舒服？” 现在我要把它用我自己身上了, 我决定一篇BLOG如果只是写了一部分, 我也要把它PO出来, 让自己尽快地补完它, 而不是过一段时间就不想再写了

Written by wOOL

18/11/2009 at 20:29

Posted in Calvados

戏言

with 2 comments

我最近看的一本轻小说的名字 以碎碎念和吐槽为亮点 半个月没更新blog 是因为我正在勤奋的更新twitter 今天终于找到了一个用起来很爽的广告过滤工具 CHROME完美了…
今天重点是帖这个工具

http://ruzanow.ru/extensions/updates.xml

可以基于页面代码过滤 也可以基于URL过滤 很强大 对页面打开速度没什么大的影响
不过不同于PRIVOXY这样前置的软件 它并不能节约流量
但是基于CHROME对页面DIV TD之类页面元素方便的选择功能 这个过滤用起来实在是太爽了
明天有2012看 6.X胖子 不错不错 木哈哈

Written by wOOL

16/11/2009 at 04:24

Posted in Absinthe

Tagged with

吃… 拉面 的 问题

without comments

我在万圣节的晚上饿着肚子来计算一个关于拉面的问题真是杯具啊

一乐拉面馆, 一共有4种口味的拉面, 老板每天会选1种供应, 鸣人每天只去吃1次拉面, 问平均多少天鸣人才能尝遍所有4种拉面呢?

一个数学化的描述: 假设第N天吃到了第4种拉面, 求E(N), E(N) = \sum _{n=4}^{\infty } n P(N=n)

下忍级解法:
显然的最少需要4天
需要4天:
所有可能数 4^4
符合要求的可能数 4*3*2*1
P(N=4)= 4*3*2*1 / 4^4
需要5天:
所有可能数 4^5
符合要求的可能数: 4*(3^4 – 3 – 3*(2^4 -2))
P(N=5)= 4*(3^4 – 3 – 3*(2^4 -2)) / 4^5
首先从4个中选择一个做”第4种”, 剩下3种, 前面4天总的可能有 3^4 种, 去掉只吃了1种拉面的情况 3 种, 去掉只吃了2种拉面的情况: 从3种拉面中选择这2种, 4天中吃2种拉面的可能为2^4 -2 (总可能2^4 去掉只有1种拉面的2种)
…
需要n天:
所有可能数 4^n
符合要求的可能数: 4*(3^(n-1) – 3 – 3*(2^(n-1) -2))
P(N=n)= 4*(3^(n-1) – 3 – 3*(2^(n-1) -2)) / 4^n

最后整理得到的求和公式为\sum _{n=4}^{\infty } \frac{3n\left(3^{n-2}-2^{n-1}+1\right)}{4^{n-1}}

瞄一眼发现第无穷项貌似是0, 这个级数有可能收敛, 但是我饿了, 所以直接丢Mathematica里, 解得答案25/3

影级别的解法:
如果某一事件发生的概率为p, 那么平均1/p天发生一次
第一天吃第1种面, 其后
其后吃到第2种拉面概率为3/4, 那么平均4/3天吃到第2种拉面
其后吃到第3种拉面概率为2/4, 那么平均2天吃到第3种拉面
其后吃到第4种拉面概率为1/4, 那么平均4天吃到第4种拉面

1 + 4/3 + 2 + 4 = 25/3

Written by wOOL

01/11/2009 at 04:05

Posted in Calvados