wOOL's Blog

Lethal Sweety . Sunny Yawn . Deicidal Jujube . …

Archive for the ‘Calvados’ Category

rpy2 – 结合Python与R

without comments

Python被称作胶水语言, 她和N多语言都能结合的很好, 在Python中引入R使用的模块是rpy2, 可以从http://rpy.sourceforge.net/下载
安装后尝试
import rpy2.robjects as robjects
如果提示找不到R.exe, 需要在Windows环境变量中添加R_HOME, 值为R的安装目录, 我的是C:\Program Files (x86)\R\R-2.9.2
如果提示没有win32api模块, 需要安装pywin32
大致上, 有2种调用R的方法
其一是直接使用R代码, 比如随机选取10个标准正态分布的变量
R Code
rnorm(10)
Python

  1. import rpy2.robjects as robjects
  2. r = robjects.r
  3. y = r(‘rnorm(10)’)

注意这里得到的结果总会是一个RVector(尽管有时候只有一个值), 有点类似于python的list, 她当然也可以转换为list
其二是调用R的函数, 依然是上面的例子
Python

  1. import rpy2.robjects as robjects
  2. r = robjects.r
  3. y = r.rnorm(10)

下面的程序, 以标准正态分布为例, 演示了用小样本方差去估计总体方差时, 虽然使用Bessel’s correction后的无偏估计\frac{1}{n-1}\sum_{i=1}^n\left(X_i-\overline{X}\,\right)^2的确在均值方面更加接近真实值, 但是用通常的距方法得到的有偏估计\frac{1}{n}\sum_{i=1}^n\left(X_i-\overline{X}\,\right)^2, 却可以提供更小的MSE, Wiki中有更加精确的数学表达

  1. import rpy2.robjects as robjects
  2. r = robjects.r
  3. ss = r.rnorm
  4.  
  5. def test(n):
  6.     y = ss(n)
  7.     avg = sum(y)/n
  8.     s=0
  9.     for i in y: s=s+(i-avg)**2
  10.     v1=s/n
  11.     v2=s/(n-1)
  12.     return v1,v2
  13.  
  14. def test0(m,n):
  15.     v1,v2=[],[]
  16.     for i in range(m):
  17.         v1.append(test(n)[0])
  18.         v2.append(test(n)[1])
  19.     avg1=sum(v1)/m
  20.     avg2=sum(v2)/m
  21.     s=0
  22.     for i in v1: s=s+(i-avg1)**2
  23.     s1=s/m
  24.     s=0
  25.     for i in v2: s=s+(i-avg2)**2
  26.     s2=s/m
  27.     return avg1,avg2,s1,s2
  28.  
  29. t1,t2=0,0
  30. for k in range(10):
  31.     v=test0(100,5)
  32.     if abs(1-v[0])<abs(1-v[1]): t1=t1+1
  33.     if v[2]<v[3]: t2=t2+1
  34.  
  35. print "有%d次,有偏估计的期望更接近真实值" % t1
  36. print "有%d次,有偏估计的方差更小" % t2

Written by wOOL

05/12/2009 at 05:32

Posted in Calvados

Tagged with , , , , ,

Google API for LaTeX to Image

without comments

Here is a Google Chart API which is used to render your LaTeX code to an image.
http://chart.apis.google.com/chart?cht=tx&chs=1×0&chf=bg,s,FFFFFF00&chco=000000&chl=$$$e^{i\pi}%2B1%3D0$$$
Use the HTML IMG tag wraping above url and it looks
e^{i \pi} + 1 = 0
For someone who doesn’t use WordPress XD

Written by wOOL

01/12/2009 at 20:43

Posted in Calvados

Tagged with

不相关且不独立的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

吃… 拉面 的 问题

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

德·梅齐里亚克的砝码问题

with one comment

碎碎念 关于睡前看实分析有助于睡眠一事的确合理 尤其是窝在被窝里 但是 隔壁的老外如果叫得很基情的话 就… 很难想象3个男人在一起为什么会发出那种声音哦 于是我起床 于是我更新blog

关于这个梅氏砝码的问题, 貌似出自1624年出版的Problèmes plaisants et délectabled qui se font par les nombres一书, 作者是法国数学家G·B·德·梅齐里亚克(Gaspard Bachet de Méziriac，1581 – 1655)

问题描述如下

一个商人有一个40磅的砝码，由于跌落在地而碎成4块。后来，称得每块碎片的重量都是整磅数，而且可以用天平和这4块来称从1至40磅之间的任意整数磅的重物。
问这4块砝码片各重多少？

这个问题我所能理解的最为简单的解法如下: 对于每一个砝码, 都有3种状态, 即 放在天平左盘 放在天平右盘 不放在天平上, 那么这个问题就可以被理解为 找到一个三进制的计数系统, 那么每个砝码的重量必须是3的n次幂, n=0, 1, 2, 3 可知 砝码重量为1 3 9 27 和恰好是40且可以组成1~40所有的数

当然上述答案是比较难想的(精妙的解答往往在于此, 想出来很难, 理解起来却容易), 下面有一个比较容易想到的解法, 摘自100 Great Problems of Elementary Mathematics:Their History and Solution

天平的两个秤盘可区别为砝码盘和称量盘，在砝码盘上只放砝码，而在称量盘上放重物外还可附加砝码。若想设法用最少块数的砝码去称量，就要把砝码也放到称量盘上。
假如任意取出几块砝码放在盘上，例如，在一个盘上放 5 磅砝码和 10 磅砝码各一块，另一个盘上放 1 磅、3 磅、4 磅的各一块，那么这些砝码便使前一个秤盘偏重 7 磅。
我们只考虑重物和砝码均为整数，也就是说，重物和砝码的重量均为整数磅。
假如有一系列砝码 A ，B ，C，…，把它们适当地分放在两个盘上，就能称出从 1 到 n 的所有整数磅的重物。如果有一块新砝码P ，它的重量p 超过原有砝码的重量总和n，超过 量为原有砝码重量的总和加 1：
p – n = n + 1，
或者
p = 2n + 1，
那么，把砝码P 加入砝码组A 、B 、C、…之后就能称出从 1 至p + n = 3n + 1 的所有整数磅的重物。
事实上，原有砝码组足以称出所有从 1 至n 磅的重物，为了称出 1 个p + x 或p – x 磅的重物，这里x 表示从 1 到n 的任一个数，把砝码P 放在砝码盘上，再把砝码A ，B ，C，… 分别放在两个盘上，使砝码盘或称量盘上的重量偏重x 磅。
此法成立后，这个题目就容易解答了。
为了使两个砝码A 和B 能称出最多重量，A 必须是 1 磅，B 必须是3 磅，这两个砝码能称出 1，2，3，4 磅的重物。
如果选第三块砝码 C，使它的重量
c = 2 × 4 + 1 = 9 （磅）， 那么用A ，B ，C 三块砝码就能称出从 1 至c + 4 = 9 + 4 = 13 磅的所有重物。
最后，如果选第四块砝码D ，使它的重量
d = 2 × 13 + 1 = 27 （磅），
那么，这四块砝码A ，B ，C，D 便能称出从 1 至27 + 13 = 40 磅的所有重物。
结论：这个砝码的四块碎片的重量分别为 1，3，9，27 磅。

基情貌似结束了 我去睡觉了 困啊…

Written by wOOL

27/10/2009 at 10:03

Posted in Calvados

Tagged with ,

SAS & R

with 3 comments

“喏, 这便是行业头号软件SAS了, 24万淫民币, 一个子都不能少”
“…”

关于对比 I
安装文件 SAS 3G R 30MB
核心代码 SAS JAVA R C
价格 SAS 24万 R Free

那我为什么要用SAS呢? 答曰: 看上去比较专业 吐槽曰: R看上去更加阴霸

关于对比 II
导入CSV文件 显示数据 显示二维图 进行线性回归计算
SAS代码
proc import datafile="d:\doc\data.csv" out=mydata dbms=csv replace;
getnames=no;
proc print data=mydata;
proc plot data=mydata;
plot VAR2 * VAR1;
proc reg data=mydata;
model VAR2 = VAR1;
run;

R代码
mydata = read.csv("d:\ \doc\ \data.csv", header=FALSE)
mydata
plot(mydata)
mylm = lm(V2~V1, data=mydata)
summary(mylm)

我恨分号我恨分号我恨分号, R更加UNIX, 大小写敏感, 在win下需要双斜线, 很莫名其妙的一点, 如果csv文件的结尾不是一个空行的话, R会有一个warning说最后一行不完全, 但不影响数据获取, SAS则没有这个问题

关于对比 III
折磨一下吧
数据量10 million, CSV文件体积244 MB(其实我开始生成的是100 million条, 2.38 GB的庞然大物, 结果发现我电脑吃不消), 这些数据大概是一个y=0.123456789x+987.654321的回归线, 每个值的误差在[-13,13]内

SAS的表现是7s读取完后瞬间算出y=0.12346x+987.65698
而R在吃掉800MB读了N分钟后读完数据然后计算, 内存占用瞬间飙升1.X G然后提示超过最大可用内存, 罢工了, 杯具啊

Written by wOOL

05/10/2009 at 08:06

Posted in Calvados

Tagged with ,

CUDA & Monte Carlo Method

with one comment

当年在选笔记本的时候, 一定要买NV卡的原因就是NV有CUDA, 而我的陈旧的GeForce 7300并不支持这个功能
简单的说, CUDA是一个C/C++语言的编译器, 由她编译出来的程序可以利用GPU强大的浮点运算能力, 提高运行效率
目前比较常见民用级的应用主要集中于视频压缩
当然, 在金融计算这样的领域, 它的价值应该得到更好的体现, 一堆GPU连起来起来比N堆CPU连起来还要强N多的样子

在Windows 7平台下搭建CUDA开发运行环境还是非常简单的(相对于Linux), 首先在NV CONTROL PANEL -> SYSTEM INFORMATION -> COMPONENTS 下看看CUDA的Driver版本号是多少, 我的是显卡是Geforce GT 240M, 目前可以用的驱动带的CUDA Driver版本是2.2, 好像多数桌面显卡都有2.3驱动支持, 在确认有CUDA支持的驱动装好后, 依次安装对应版本的CUDA Toolkit 和 CUDA SDK (http://www.nvidia.com/object/cuda_get.html), 然后安装Visual Studio, 需要注意的是, 虽然Visual C++ Express可以支持CUDA, 但是由于其默认不支持64bit编译, 所以需要下载安装.Net Framework和替换, 所以还是推荐安装Visual Studio, 我安装的版本是2008 Pro, 安装的时候务必选择Full安装, 以保证64bit compiler装上, 安装完后安装CUDA VS Wizard (http://sourceforge.net/projects/cudavswizard/), 这样可以在New Project里面看到CUDAWinApp, 至于代码高亮, 发现目前只支持VS 2005, 反正我一般都是在Notepad++里写代码, 然后扔VS里调试编译.

简单安装流程:
1.安装最新显卡驱动并查看CUDA DRIVER驱动版本
2.下载对应的CUDA Toolkit 和 CUDA SDK
3.Full安装Visual Studio
4.安装CUDA VS Wizard

注: 安装好CUDA SDK后可以运行C:\ProgramData\NVIDIA Corporation\NVIDIA CUDA SDK\bin\win64\Release中的程序以测试是否安装了正确的驱动和Toolkit

搭建好CUDA的平台后自然要测试一下, 刚好拿最近在看得Monte Carlo Method试试, 几乎所有的Monte Carlo Method入门介绍都会提到用这种方法计算pi, 其C语言的核心代码非常简单

int i, pi=0;
double x, y, area;

for (i=1; i<=1000000; i++){
x= rand()/16384.0-1;
y= rand()/16384.0-1;
if ((x*x + y*y)<1) pi++;
}

area=4.0*(double) pi/(double) i;
printf("%f \n", area);

CUDA是支持标准的C语句的, 但是它有一些自己的.h文件, 比如cutil_math.h (对应math.h)
进行100万次循环, 纯C和CUDA运行时间分别为56ms和51ms
貌似没有明显差别, 下次打算拿更复杂一些的代码来测试一下, HOHO~

Written by wOOL

05/10/2009 at 03:21

Posted in Calvados

Tagged with ,

[Puzzle]Names in Boxes

without comments

地球上的某日:
我: 唉, 今天忘记带眼镜了, 看不到PPT
桶: 啊, 你眼镜多少度的呀
我: 250(事实)
桶: 不会这么低吧?
我: 我真的是250(干脆滴)
众: 扑哧~~~(笑…)

——————-正文滴分割线———————–
有一个监狱有100个囚犯, 他们的名字被写在一张纸上, 放在外形相同的盒子中, 这些盒子有1~100的编号, 它们被整齐在摆在一个房间里, 现在这些囚犯一个一个地进入这个房间, 进入房间后他们可以打开至多50个箱子, 如果其中有自己的名字, 则表示他”成功”了. 只有所有囚犯都”成功”了, 他们才可以被全部释放, 否则则会被人道毁灭.
1. 所有的囚犯依次进入房间且不可在打开箱子之后和其他囚犯进行任何交流
2. 对于盒子只有打开是允许的, 不能移动他们, 不能拿走里面写有名字的纸

问是否有一种策略, 使他们有30%以上的机会可以被释放呢?

直观上这的确是很难的, 按照最普通的方法, 每个人随机取50个箱子, 拿到自己的那个的概率是50/100=1/2, 他们能被释放的概率只有(1/2)^100, 无疑远小于30%, 但是利用如下的策略却可以达到30%的概率:

制作一个表, 把囚犯的名字随机编号1~100, 每个囚犯带着这张表, 他们中的某个人进入房间之后, 首先打开编号为自己名字对应的编号的那个箱子, 如果箱子里写的不是自己的名字, 则按照这个名字, 找到表上对应的编号, 然后再去打开此编号的箱子, 以此类推, 直到找到自己的名字或者已经打开了50个箱子为止.

解释
一般化: 取n=50
2n个名字可以组成的排列有(2n)!个, 假设这些囚犯需要打开k(k>n)次才能找到各自的名字, 这种情况下可以组成的排列如下:
从2n个名字中取k个, 有(2n)!/(k! * (2n-k)!)种可能, 这些排列中有一个(即囚犯自己的名字)位置是确定的, 剩下的(k-1)个名字可以自由排列有(k-1)!个, 剩下的(2n-k)个也可以自由排列, 又有(2n-k)!个, 把这个3个数相乘是(2n)!/k, 占所有可能性(2n)!的1/k
那么囚犯至多需要n次就可以找到各自名字的概率
1-1/(n+1)-1/(n+2)-……-1/(2n) = 1 -(ln 2n – ln n) = 1 – ln2 = 0.31

Written by wOOL

24/09/2009 at 08:11

Posted in Calvados

Tagged with ,

« Older Entries