合并石子问题(1)

mathe

这个题目比较有意思的是，最后一步融合时，两个石子分别从k个和n-k个石子融合过来的概率都相等，也就是k!=n-k时都是2/(n-1)，而对于k=n-k,概率为1/(n-1)
假设开始2N个1的状态最后到达k的概率为$p_k(2N)$,$s_k(2N)=\sum_{h=2}^kp_h(2N)$
于是我们有$p_2(2N)=N/{2*N-1}+1/{2*N-1}\sum_{h=1}^{N-1}(p_2(2h)(1-s_2(2N-2h))+(p_2(2N-2h)(1-s_2(2h)))$
$p_k(2N)=1/{2*N-1}\sum_{h=1}^{N-1}(p_k(2h)(1-s_k(2N-2h))+(p_k(2N-2h)(1-s_k(2h))+p_{k-1}(2h)p_{k-1}(2N-2h))$

wayne

哈哈，我想到一个更高效的方法了，我再试试。
看能不能算出
p(102)
p(200)
嘿嘿

mathe

P[400]=1353552604751171454515088436564590520732246019178768149045833000434288749
76607157339342685426434476202342477151608501836327284200201873205568270639676082
06269943972300528019372479612341831147648852565783623358428454588653473723128313
57296519476096885322635218376803108286224769999017087295090804939570045365015056
11654114045590346853040539647191696811120962979691828881038934199237852714005622
72531939209144563570051830553971340808728209383785966872058381282534569902989722
06259865563624017291368022499673051854454363611408596730659473397763582736542798
88280753884780153088916226011182930739265721972318847621868961938707725432967987
31064909663357210780926432281070063758536924263167815254540841805656467858230365
1256440491751088759489268869809/581603970588045117509153218767440541546942959836
65710735650133446091206392857914868534376162568532488496798542439173474358241847
15232482107689655889260279485245169600765614815127664554253855433815842289847901
94373898808965365007231243803069001729357211027325310936221847680801585011740248
44019652090273250905226483087641863009482750054995585612988607781346937120070117
06957362943577416224061583731682498411705522044648519465377477393519374709585871
11083805108399840694686568047016259007381046716961390046121305452294401404745617
63114378149666737671276928213165903606975005398921823133546785482311473516062484
47985203754139367392564231726218497139446106629452714325518390427231424845458377
4275833058030153732698153135061147622764110565185546875~=2.327275385315249012281
025172

mathe

猜测极限值为:
2.327529250867897115016845370858722697323295012330077167951168532096701369
84224874490409159837237145849799749962293598707192569416304117941915450365977242
56874754841025996893012765855243738442437294836105368197937889847238602724879852
52335527910003828971033397229102093373458076117416591392462939161159770200826180
54493594698886027481478198201439774861651420722863594763799813848894513040556801
44129169407889087161969548630443765703805969547479879753094079154936685896625531
69461968419707802284382698923163831891654672262050878385476180806017736487247429
94601993333469893330446121050047491462521064133923730201624903226479353813224480
48244641201408238088983000156915876233821465814998801990000081439981213542811450
88556922705609290868335521744859787275185306230675164989877723148617316840036252
73539068488506366016005202399157871380032162949400030258937276933557411375759585
99398418038013500463235835356464831941386586680681473979477751105832864727366348
97682271762674982400102481123209142838148294634

mathe

道理如下。我们可以认为对于充分大的2N，其最终达到2,3,...,等数值的概率分别为$p_2,p_3,...$,这些值也会稳定下来，于是$p_2+p_3+....=1$
由于我们已经得到对于充分大的2N,
其中有$N/{2N-1}~=1/2$的概率是最后一步合并是两边数都是奇数个数合并过来的，最后达到2。
另外，分别有$1/{2N-1}$的概率最后一步两个是(2k,2N-2k),对于大部分k,我们知道其最终达到各个数的概率也应该近似是$p_2,p_3,....$，所以我们可以认为有大概$1/2$的概率(${N-1}/{2N-1}$)是将两个分布都是$p_2,p_3,...$的状态合并
于是我们可以估计,对于极限分布
$p_2=1/2+1/2*2*p_2*(p_3+p_4+...)=1/2+p_2*(1-p_2)$,得到$p2=sqrt(1/2)$
同样，我们可以得出对于k>2,有
$p_{k+1}=1/2p_k*p_k+1/2*2*p_{k+1}*(p_{k+2}+....)=1/2p_k*p_k+p_{k+1}*(1-p_2-p_3-...-p_{k+1})$
由此可以计算出极限分布概率
可以利用递推公式
$p_2=s_2=sqrt(1/2)$
$p_{k+1}={sqrt(s_k*s_k+2*p_k*p_k)-s_k}/2$
$s_{k+1}=s_k+p_{k+1}$
然后最后计算期望极限为
$\sum_{k=2}^{infty}k*p_k$
另外，我们可以估计出对于充分大的k，由于s_k大概为1，于是我们可以近似有$p_{k+1}~={p_k^2}/2$,所以可以知道，计算公式收敛非常快

wayne

23# mathe
额，
mathe够绝的。
这让我等小生 怎么才能望你的项背啊

wayne

25# mathe
精彩。
我待会腾出时间，看看我的新改进的想法能算到什么程度

mathe

猜出极限值只是跨出万里长征的第一步。
现在可以证明p2的极限存在了。假设2N个石子最后变成两个的极限为$p_2(N)$
假设$lim_{N->infty}p_2(N)$不存在，由于概率显然有界，那么必然有下极限a和上极限A，而且显然a和A都不小于$1/2$,由于函数$1/2+x-x^2$在$1/2<=x<=1$减，我们可以得出
$a>=1/2+A-A^2,A<=1/2+a-a^2$
另外由于容易证明$|1/2+x-x^2-sqrt(1/2)|=|(x-sqrt(1/2))|*|(1-sqrt(1/2)-x)|<|x-sqrt(1/2)|$
容易得出
$|a-sqrt(1/2)|>=|A-sqrt(1/2)|>=|a-sqrt(1/2)|$等号只有$A=a=sqrt(1/2)$取到

mathe

对于$k>=3$,我们如果能够证明$p_k(2N)>=2\sum_{h=k+1}^{infty}p_h(2N)$,那么也就应该能用类似方法来证明极限存在了,这个不等式应该不难

wayne

23# mathe
我的代码也出来了，设计了新的数据结构，但算到p[100] 还是不给力，因为没有递归存储小的结果，所以还有很大的提升空间，
， 不折腾了，还是学学mathe前面贴的代码。
mathe的程序太精彩了，打印结果简直是刷刷的直接滚屏。

1. n = 50; x = {{{{1, n - 2}, {2, 1}}, 1}};
   
2. g[x_] := If[Length[x] > 1,
   
3.    If[x[[2, 2]] == 1, {x[[1]]},
   
4.     x + {{0, 0}, {0, -1}}], {x[[1]] + {0, -2}, {x[[1, 1]] + 1, 1}}];
   
5. Do[x = Table[{ii[[1, 1]], Total[ii[[All, 2]]]}, {ii,
   
6.     GatherBy[
   
7.      Flatten[Table[
   
8.        Join[Table[{Table[{jj[[1, 1]], Total[jj[[All, 2]]]}, {jj,
   
9.             GatherBy[
   
10.              Sort@Select[
    
11.                Join[DeleteCases[tmp[[1]],
    
12.                  ii], {ii + {0, -2}, {ii[[1]] + 1, 1}}], #[[2]] >
    
13.                  0 &], First]}], Binomial[ii[[2]], 2]*tmp[[2]]}, {ii,
    
14.           Select[tmp[[1]], #[[2]] > 1 &]}],
    
15.         Table[{Sort@
    
16.            Join[DeleteCases[DeleteCases[tmp[[1]], ii[[2]]], ii[[1]]],
    
17.             g[ii]], (ii[[All, 2]] /. List -> Times)*tmp[[2]]}, {ii,
    
18.           Subsets[tmp[[1]], {2}]}]], {tmp, x}], 1], First]}], {n - 2}]

复制代码

1. ans = Partition[Flatten[x], 3];
   
2. Total@Table[ii[[1]]*ii[[3]], {ii, ans}]/Total[ans[[All, -1]]]

复制代码