随机游走中的概率问题

mathe

A和B一开始站在同一个地方，他们不停地猜拳，A赢了就前进1米，B赢了就前进$pi$米（他们朝同一个方向前进）直到A前进到B的前面为止，求A走到B前面的概率.
转自http://tieba.baidu.com/f?kz=42429501
显然，这个概率很难用解析表达式表示出来，那么现在问题是如何用数值计算的方法找到一个比较好的近似值（比如精确到小数点后面10位数？？？）

无心人

pi是多少啊？

用计算机模拟吧

mathe

随机模拟，也就是模特卡罗法精度太低了，而且很难保证。
现在通过http://bbs.emath.ac.cn/thread-332-1-1.html中的结论，我们可以构造出一种比较合适的解决方案。

mathe

我们先将问题推广到A赢了前进1米，B赢了前进x米(x>1).
容易看出来，这个问题A能够超过B的概率p(x)是关于x的单调减函数。而我们现在要计算$p(pi)$.
为此，我们可以选择两个非常接近$pi$的有理数a和b,$a>pi>b$,分别计算$p(a)$和$p(b)$,那么必然有$p(a)<p(pi)<p(b)$
在a,b都很接近$pi$时，这两个结果应该可以很好的逼进$p(pi)$.
所以问题转为对于有理数$n/m$,计算$p(n/m)$
而这个问题可以转化为A赢了前进m米，B赢了前进n米的问题(n>m,(m,n)=1)

mathe

对于A前进m，B前进n的问题，我们可以知道A和B之间的间隔始终是整数。
而且任何时候A能够超过B的概率之和他们间相对距离有关系。我们记当A在B后面k米时，A能够超过B的概率为$q(k), (k>=1)$
而且我们可以定义对于$k<=0, q(k)=1$
于是我们有递推式$q(k)=1/2 q(k-m)+1/2 q(k+n), (k>=1)$
这个递推式对应的特征方程为$x^{m+n}-2x^m+1=0$. 记$y=1/x$,那么y满足方程$y^m=2-y^{-n}$
根据 http://bbs.emath.ac.cn/thread-332-1-1.html中的结论，正好有m个范数大于1的y满足上面的方程，
(这里也可以看出为什么我前面会有那样的猜想，只有这样才合理，我们这里只有m个初试条件，所以通项公式里面应该只有m个候选系数）
也就是特征方程有m个范数小于1的x。此外特征方程还有n个解的范数不小于1。
假设特征方程m个范数小于1的解为$x_1,x_2,...,x_m$,范数不小于的解为$x_{m+1},...,x_{m+n}$
由此它的通解形式为
$q(k)=a_1 x_1^k +a_2 x_2^k+...+a_{m+n} x_{m+n}^k$
由于q(k)有界 ($0<=q(k)<=1$)，而且容易证明，当k趋向无穷时,q(k)趋向0,然后根据：http://bbs.emath.ac.cn/thread-354-1-1.html
我们可以得到
$a_{m+1},a_{m+2},...,a_{m+n}$都是0.
所以$q(k)=a_1 x_1^k+a_2 x_2^k +... +a_m x_m^k$
而$a_1,a_2,...,a_m$我们可以通过$q(0)=1,q(-1)=1,...,q(-m+1)=1$使用待定系数法解出。
而在这个计算过程中，我们记$y_i=1/x_i$,所以$|y_i|>1$
那么对应线性方程组左边的矩阵正好时范得蒙行列式，我觉得应该可以有解析表达式，只是还没有进一步去分析。
在得到$a_1,a_2,...,a_m$的解以后
我们可以计算出$q(n)=a_1*x_1^n+a_2*x_2^n+...+a_m*x_m^n$
而最终所求的概率为$1/2+{q(n)}/2$
而上面问题的算法时间复杂度最大为O(m^3)（解m阶线性方程组的时间复杂度）。而如果解出解析解，应该有更加低的时间复杂度。

mathe

算了一下，范得蒙矩阵对应的线性方程组可以在O(m^2)内解决。
所以只要m不是太大，计算不是问题，倒是浮点运算容易溢出，所以不得已我才用gmp进行计算。
我分别采用m=22/7, 355/113, 333/106来计算，得到的概率都是0.543643。所以这个就应该$pi$倍移动对应的解了。
只是现在gmp里面对于浮点运算中三角函数不支持，所以前面求$x_i$部分我还是采用原先迭代方法，需要三角函数运算，所以就没有采用gmp计算。所以现在精度也不能无限提高，只能到这个精度了。

无心人

支持高精度超越函数的数学工具少的很

mathe

呵呵，印象中Mathematics支持。
还有Linux下面的bc工具就支持，估计perl也支持。只是这两个东西实在太慢了一些

无心人

Mathematics性能也不咋地

mathe

精确一点应该是0.54364331158
333/106计算结果为
0.543643311581
355/113计算结果为
0.543643311576
不过就是担心前面解方程的根的精度有点不够。