一维空间中的吃豆子问题

mathe

找了一下Pari/Gp里面有intformal函数可以对多项式积分，不过里面要求所有符号有个默认顺序，而积分必须对默认符号积分，其实起来不是很方便，所以下面的代码使用前必须保证gw,gc变量没有定义，或者gw定义在gc之前，也就是调用函数reorder()可以看到gc没有定义或者gw,gc都定义了但是gw在前面

2. symbolp(n)=
3. {
4. local(p,q,p1,p2,p3);
5. global(gw,gc,gt1,gt2);
6. p=1;
7. for(u=1,n,
8. q=intformal(p,gw);
9. q=subst(q,gc,gt1);
10. q=subst(q,gw,gt2);
11. p1=subst(q,gt1,gc);
12. p1=subst(p1,gt2,(gc+gw)/2.0);
13. p2=subst(q,gt1,1.0-gw);
14. p3=subst(p2,gt2,1.0-gw);
15. p2=subst(p2,gt2,1-(gc+gw)/2.0);
16. p3=p3-p2;
17. p2=subst(q,gt1,gw);
18. p2=subst(p2,gt2,gw);
19. p=p1+p2+p3;
20. );
21. p
22. }
23. getdist(n)=
24. {
25. local(p,d1,d2,d3,d4);
26. global(gw,gc);
27. p=symbolp(n);
28. d1=p*(-gc^2/2.0+gw^2+gc-gw);
29. d2=p*(-gc^2+2*gw*gc-gw^2+gc-gw);
30. d3=intformal(d1);
31. d4=intformal(d2);
32. d1=subst(d3,gw,gc/2.0);
33. d2=subst(d4,gw,gc);
34. d4=subst(d4,gw,gc/2.0);
35. d2=d2-d4;
36. d1=d1+d2;
37. d1=intformal(d1);
38. 2.0*subst(d1,gc,1.0)
39. }

复制代码

然后计算结果如下: (16:06) gp > getdist(2) %17 = 0.2114583333333333333333333333 (16:07) gp > getdist(10) %18 = 0.2203686807254583318975768873 (16:07) gp > getdist(20) %19 = 0.2204074196271448826430741941 (16:07) gp > getdist(30) %20 = 0.2204077318388700882261725891 (16:07) gp > getdist(40) %21 = 0.2204077356077135199957228968 (16:07) gp > getdist(50) %22 = 0.2204077356585942289207798388 (16:07) gp > getdist(60) %23 = 0.2204077356593274980124179385 (16:07) gp > getdist(70) %25 = 0.2204077356593385397246878697 (16:33) gp > getdist(99) %1 = 0.2204077356593387142224517526 (16:10) gp > getdist(100) %26 = 0.2204077356593387141426361209 也就是我们可以得出g(a,b)=(a+b)/2时的高精度结果，也就是每秒4.537次,同Fans的模拟结果匹配

mathe

14#里面zgg定义的gp(x1,x2)=(x1+x2+q/2)/(q+2)好像有问题，比如在极端情况x1~=x2时，gp(x1,x2)应该也取x1或x2,但是他这里定义的函数不满足，所以结果应该不对。不知道是不是他定义的gp函数不是我定义的g函数？ 而对于Fans前面定义的情况，由于对应的g函数太复杂，符号计算太难。最好我们改成线性的，至少弄成折线，比如(a+b)<1时，g(a+b)=c*a+(1-c)*b,不然g(a+b)=(1-c)*a+c*b,其中c是一个大于1/2的数，这样定义的g就应该偏向于选择离1/2远的点。反之如果c小于1/2，那么偏向于选择离1/2近的点。 当然，由于g表达式更加复杂了，对应的符号计算代码也要修改了

mathe

发现只要g是分段，那么符号计算的结果就会越来越复杂，不是很合适

mathe

由于对于一般的g,直接积分计算很困难，我们可以考虑对p进行多项式拟合。也就是假设p是关于w,c的n次多项式(O(n^2)个系数），然后对p采样（采样数目大于待定系数），代入等式得到线性方程组，当然方程数目大于变量数目，Ax=b,然后解方程A'Ax=A'b应该可以得出一个比较好的p(w,c)

mathe

上面多项式可能不如改成三角函数，不然计算可能会比较变态。唯一比较麻烦的是在g比较复杂时，最后一步平均移动距离的计算公式可能比较复杂，实在不行，这一步可以采用数值积分。

mathe

而34#中得出g(a,b)=(a+b)/2时最优应该是错误的方法。理由只能是p和g是相互约束的，给定g即唯一确定p,同样给定p就唯一确定g,所以我们不能使用假设p是给定的条件求出g的最优值。 而34#方法得出错误的结论说明32#中迭代的方法也不能使用，因为它也必然会收敛到34#中的结果。

mathe

又弄了一个10^8的豆的，感觉上需要验证随机数的可靠性和考虑double的误差积累了。所以请Fans也有空算一下权重p的极小值吧。呵呵。 源码如下，我用vc2005：#include "stdafx.h" #include #include #include t ... zgg___ 发表于 2013-2-19 15:45

19#中经常会出现越过一个豆子吃另外一个豆子，去除这种情况，可以得出一秒可以超过5.9次

zgg___

这些天抽空对25层和30层的方法又检查了一下，觉得有的结论是对的，呵呵。 先说下我在30层为什么认为存在问题？ 主要有二：Q1、25层中对平均吃豆时长的计算结果和试验结果相差大，比如说，对p=0时，即w总去吃最近的豆的情况，此时按照30层公式8的结果得0.2185013458，而试验结果是0.220+（和试验结果差距这么大，一定是有问题的呢，这个不光是Fans的试验，我自己的程序也是这个结果，呵呵。）Q2、是用25层公式7算出的曲线和30层代码第3行的结果是不一致的，不一致性在30层最后输出的图里可以看到，因为都是M的计算，所以和试验都搭不上边，也一定有问题的呢，呵呵。综合这两点，我也对“模型”产生过怀疑，呵呵，但是，我并不觉得“模型”有什么问题，可能在计算中会有计算错误，但是“模型”还是比较要命的，呵呵。 对于Q2的问题，发现了问题所在（正如mathe在42层说的），是因为在计算中没有考虑到g有可能落在x1和x2之外，比如说，当出现1/2

KeyTo9_Fans

听说$19#$的策略会越过一个豆子吃另外一个豆子，立马打个补丁重新运行。 目前已运行$48$小时，结果为： $p=0.478\pm0.005$ $t=0.2191031\pm0.0000005$ 其中$p$是$1/2$的权重，$t$是吃一个豆的平均时长。

wayne

49# KeyTo9_Fans 这是不是 意味着zgg的就是终极答案了？