一维空间中的吃豆子问题

mathe

我觉得是因为模型是错误的，所以计算结果会和试验结果不同

mathe

23楼方法比如离散成100个点，g,F,W的规模都在5000个点左右。于是积分变成求和（边界可以加上1/2的权值），给定g求F,W变成求线性方程组而给定F,W求g是一个线性规划问题。但是同时作为变量是非线性规划，很难有很好的方法。于是我们可以试验迭代.如果迭代过程不收敛，可以试着使用松弛因子法

mathe

推导了一下，发现对于稳定状态下，最优的时候还是选择g(a,b)=(a+b)/2 我觉得计算机模拟发现(a+b)/2不是最优解的原因是开始的时候，密度分布函数还没有稳定，我们不妨试着把前面部分数据抛弃以后再看看

mathe

我们可以认为使用某个策略g(a,b)经过充分长时间后，分布密度函数为p(a,b,w),其中a,b是两个豆子的坐标0<=a<=b<=1,而w是wayne的坐标。 于是平均长度为 $\int_0^1\int_a^1 C(a,b)dbda$ 其中 $C(a,b)=\int_0^a (a-w)p(a,b,w)dw+\int_a^{g(a,b)}(w-a)p(a,b,w)dw+\int_{g(a,b)}^b(b-w)p(a,b,w)dw+\int_b^1(w-b)p(a,b,w)dw$ 我们注意到上面表达式中对于给定的a,b,第一项和最后一项是常数，而中间两项取决于g,而对于给定a,b有唯一的g(a,b)使上面表达式最小，而这个表达式关于g的导数是 $(2g-a-b)p(a,b,g)$ 所以$g={a+b}/2$时取到最小值,这时二阶导数取值为$2p(a,b,(a+b)/2)>0$，所以是极小值点而不是极大值点

mathe

我现在用Pari/Gp大概计算了一下。但是离散计算时，网格分割的多一点，Pari/Gp就堆栈溢出了。 现在在Windows下只能计算到15×15的网格，所以精度估计有点低。 zgg的表示方法不错，也就是假设给定g(a,b)策略后，稳定下来后，每次吃掉一个豆子时，余下一个豆子c和wayne的位置w可以有一个分布函数p(c,w),而新出来的豆子位置u是满足均匀分布的，可以先不予理会。 而由对称性我们知道p(c,w)=p(1-c,1-w)，所以我们可以光考虑0

mathe

而得出$p(c,w)$以后，在给定$g(a,b)=(a+b)/2$时计算平均移动位置比较简单。我们同样只考虑$0c$时 $d(c,w)=\int_0^{2w-c}(c-w)du+\int_{2w-c}^w(w-u)du+\int_w^c(u-w)du+\int_c^1(c-w)du=-c^2+2wc-w^2+c-w$ 在$2w

mathe

对应Pari/Gp代码，不过由于n最大只能使用15,计算结果精度不高 比如 (09:15) gp > average_weight(10,200) %34 = 0.2571836267388383841198304306 (09:17) gp > average_weight(11,200) %35 = 0.2541804007241987873686495312 (09:17) gp > average_weight(12,200) %36 = 0.2513520585091040215503666987 (09:17) gp > average_weight(13,200) %37 = 0.2492168669439747355725463563 (09:17) gp > average_weight(14,200) %38 = 0.2471598401292312838397486659 (09:17) gp > average_weight(15,200) %39 = 0.2455242992227459358290402694

1. index2to1(a,b)=
2. {
3. a*(a-1)/2+b+1
4. }
6. index2wn(n,a,b)=
7. {
8. if(a<b, index2to1(a,b),
9. index2to1(n-a,n-b))
10. }
12. index1to2(n)=
13. {
14. local(a,b);
15. a=floor((sqrt(8*n-7)+1)/2);
16. b=n-a*(a-1)/2-1;
17. [a,b]
18. }
20. tcount(n)=
21. {
22. n*(n+1)/2
23. }
25. buildm(n,g)=
26. {
27. local(t,M,i,b,s);
28. t=tcount(n);
29. M=matrix(t,t);
30. for(u=1,t,
31. i=index1to2(u);
32. b=floor(g[u]);
33. if(b<i[1],
34. for(v=0,b-1,M[u,index2to1(i[1],v)]=1.0/n);
35. M[u,index2to1(i[1],b)]=(g[u]-b)/n,
36. for(v=0,i[1]-1,M[u,index2to1(i[1],v)]=1.0/n);
37. for(v=i[1],b-1,M[u,index2to1(n-i[1],n-v-1)]=1.0/n);
38. M[u,index2to1(n-i[1],n-b-1)]=(g[u]-b)/n
39. );
40. if(b<i[2],
41. for(v=0,b-1,M[u,index2to1(i[2],v)]=1.0/n);
42. M[u,index2to1(i[2],b)]=(g[u]-b)/n,
43. for(v=0,i[2]-1,M[u,index2to1(i[2],v)]=1.0/n);
44. for(v=i[2],b-1,M[u,index2to1(n-i[2],n-v-1)]=1.0/n);
45. M[u,index2to1(n-i[2],n-b-1)]=(g[u]-b)/n
46. );
47. );
48. M
49. }
51. create_g(n)=
52. {
53. local(g,t,i);
54. t=tcount(n);
55. g=vector(t);
56. for(u=1,t,
57. i=index1to2(u);
58. g[u]=(i[1]+i[2])/2.0
59. );
60. g
61. }
63. build_sm(n)=
64. {
65. buildm(n,create_g(n))
66. }
68. normalize(v)=
69. {
70. local(mv,l);
71. l=length(v);mv=0.0;
72. for(u=1,l, mv+=v[u]);
73. for(u=1,l, v[u]/=mv);
74. v
75. }
76. cost(a,b)=
77. {
78. if(a>2*b, a-a^2/2+b^2-b, 2*a*b+a-a^2-b-b^2)
79. }
80. build_sv(n, i)=
81. {
82. local(m,v);
83. m=buildm(n,create_g(n))^i;
84. v=vector(tcount(n),x,1);
85. normalize(v*m~)
86. }
88. average_weight(n,i)=
89. {
90. local(v,t,s,c,w);
91. v=build_sv(n,i);
92. t=tcount(n);s=0.0;
93. for(u=1,t,
94. h=index1to2(u);
95. s+=v[u]*cost(h[1]/(n+0.0),h[2]/(n+0.0));
96. );
97. s
98. }

复制代码

mathe

发现可以通过增加内存来扩展点的数目 (09:17) gp > allocatemem(40000000) (09:19) gp > average_weight(20,200) %41 = 0.2393806053175264654561361443 (09:21) gp > allocatemem(400000000) (09:23) gp > average_weight(30,100) %44 = 0.2331910035260221266557374996 由此可以看出这种近似计算误差还是很大，只能得出一个近似的估计值

mathe

而另外一种更好的方法是通过符号运算得出p的近似表达式，我们可以使用迭代式 $p_{n+1}(c,w)=\int_0^{g(w,c)}p_n(c,x)dx+\int_0^w p_n(w,x) dx + \int_{1-g(w,c)}^{1-w} p_n(1-w,x)dx$ 而我们比如可以选择$p_0(c,w)=1$. 而我上面的数值计算发现迭代次数对精度影响好像不是很大，所以比如我们可以算出$p_100(c,w)$,对于$g(a,b)=(a+b)/2$这应该是一个关于c,w的100次多项式，系数在$100^2$项数量级，所以问题应该不大，只是最好直接使用支持符号运算的软件，Pari/Gp就不大适合了

KeyTo9_Fans

我等小辈功力不足，看来还是需要mathe大师亲自出马才有进展。 "由此可以看出这种近似计算误差还是很大，只能得出一个近似的估计值" 虽然如此，但已经可以得到比试验结果的精确程度更高的结果了。 记 $a_10=0.2571836267388383841198304306$ $a_11=0.2541804007241987873686495312$ $a_12=0.2513520585091040215503666987$ $a_13=0.2492168669439747355725463563$ $a_14=0.2471598401292312838397486659$ $a_15=0.2455242992227459358290402694$ $a_16=...$ ...... $a_30=...$ 通过分析该数列，可以得到精确程度为$10^{-6}$以上的极限值，比试验结果的精确程度高。