合并石子问题(2)

KeyTo9_Fans

《合并石子问题(1)》已经是10年前的坟贴了：https://bbs.emath.ac.cn/thread-3978-1-1.html

本贴是《合并石子问题(2)》，对合并规则作了一些改动：

有$N$颗体积为$1$的石子。

我们每次随机抽取$2$颗石子，合并成$1$颗石子。

每颗石子被抽到的概率均等。

假设要合并的$2$颗石子的体积分别为$a$和$b$。

如果$a=b$，那么可以无损合并，合并后的体积为$(a+b)$；

否则大石子会把小石子吸收掉，合并后的体积为$\max\{a,b\}$。

一共合并$(N-1)$次，最后剩下$1$颗石子。

记这颗石子的体积的期望值为$E(N)$。

(1)对于较小的$N$，$E(N)$的精确值是多少？

(2)对于较大的$N$，$E(N)$大约等于$N$的多少次方？

mathe

这个结果只能取2的幂。每种合并过程可以看成一棵二叉树，只是不同的树的权重应该是不一样的

mathe

根据链接内容，假设最后一步两个石子的合并过程，这两个石子分别是从h和N-h棵石子合并过来的概率为1/(N-1)
所以得出$p_k(N)=\frac1{N-1}\sum_{h=1}^{N-1}(2p_k(h)s_{k-1}(N-h)+p_{k-1}(h)p_{k-1}(N-h))$
好像总是倒数第三和第四大的两个非零数出现的概率最高，两者概率加起来通常超过90%

mathe

2:2
3:2
4:2.666667
5:3.333333
10:4.519929
100:19.381907
1000:82.584912
10000:352.778718
100000:1506.949790
可以看出石子数目每增加10倍，E(N)大概增加4.27倍。所以大概$N^0.63$

KeyTo9_Fans

我没有精确计算，直接抽样模拟，结果如下：

1.          3 2.000000
   
2.          9 4.231674
   
3.         27 8.502773
   
4.         81 16.931536
   
5.        243 33.662674
   
6.        729 66.887726
   
7.       2187 132.846680
   
8.       6561 264.103027
   
9.      19683 525.009766
   
10.      59049 1043.968750
    
11.     177147 2078.093750
    
12.     531441 4132.750000
    
13.    1594323 8197.000000
    
14.    4782969 16348.000000
    
15.   14348907 32600.000000
    
16.   43046721 64384.000000

复制代码

$3^n$颗石子合并后，好像是$n$较小的时候，合并后的体积比$2^n$大，$n$较大的时候，合并后的体积比$2^n$小

mathe

1. h=3 w=2.000000
   
2.         m[2]=1
   
3. h=9 w=4.231746
   
4.         m[2]=0.0031746
   
5.         m[3]=0.937302
   
6.         m[4]=0.0595238
   
7. h=27 w=8.502055
   
8.         m[2]=8.32014e-20
   
9.         m[3]=0.00991275
   
10.         m[4]=0.922374
    
11.         m[5]=0.0677133
    
12. h=81 w=16.932300
    
13.         m[2]=8.44582e-96
    
14.         m[3]=3.45346e-15
    
15.         m[4]=0.0142225
    
16.         m[5]=0.920398
    
17.         m[6]=0.06538
    
18.         m[7]=1.95842e-13
    
19. h=243 w=33.658228
    
20.         m[3]=1.24972e-68
    
21.         m[4]=5.17132e-14
    
22.         m[5]=0.0173106
    
23.         m[6]=0.922215
    
24.         m[7]=0.0604749
    
25.         m[8]=9.90855e-13
    
26. h=729 w=66.881955
    
27.         m[3]=1.38857e-269
    
28.         m[4]=1.03378e-62
    
29.         m[5]=1.78279e-13
    
30.         m[6]=0.0200685
    
31.         m[7]=0.924867
    
32.         m[8]=0.0550648
    
33.         m[9]=9.18693e-13
    
34.         m[10]=1.65113e-88
    
35. h=2187 w=132.907027
    
36.         m[4]=2.537e-244
    
37.         m[5]=3.02525e-60
    
38.         m[6]=4.05497e-13
    
39.         m[7]=0.0228738
    
40.         m[8]=0.927353
    
41.         m[9]=0.049773
    
42.         m[10]=5.47905e-13
    
43.         m[11]=8.77733e-86
    
44. h=6561 w=264.153414
    
45.         m[5]=1.96039e-234
    
46.         m[6]=1.00005e-58
    
47.         m[7]=8.05039e-13
    
48.         m[8]=0.0258793
    
49.         m[9]=0.929332
    
50.         m[10]=0.0447889
    
51.         m[11]=2.80495e-13
    
52.         m[12]=1.92705e-86
    
53. h=19683 w=525.100205
    
54.         m[6]=1.17788e-228
    
55.         m[7]=1.6475e-57
    
56.         m[8]=1.51856e-12
    
57.         m[9]=0.0291525
    
58.         m[10]=0.930685
    
59.         m[11]=0.0401626
    
60.         m[12]=1.34692e-13
    
61.         m[13]=4.46662e-88
    
62. h=59049 w=1044.002650
    
63.         m[7]=3.59913e-224
    
64.         m[8]=2.11411e-56
    
65.         m[9]=2.79574e-12
    
66.         m[10]=0.0327297
    
67.         m[11]=0.931372
    
68.         m[12]=0.0358987
    
69.         m[13]=6.24828e-14
    
70.         m[14]=4.59493e-90
    

复制代码

mathe

根据前面Fans的观察，由于$\log_3 2\approx 0.63$, 所以结果可能会近似$N^{\log_3 2}$,
我们可以猜测$p_k(N)$会以$k-1-\log_3 N$为中心，而且对于充分大的N,这个图像会基本不发生变化。
我们猜测$ p_k(N) \to f(k-\log_3 N)$, 由此，可以对$33334\le N \le100000$范围内描出$ k-\log_3 N \to p_k(N)$的图像

如果我们再对p_k(N)取自然对数，可以描出图像$ k-\log_3 N \to \ln(p_k(N))$

KeyTo9_Fans

我前面的计算结果太粗糙了

我后来精确计算了一下，发现“石头数3倍、合并后2倍”是错觉

精确的结果应该是“石头数3.04383倍、合并后2倍”

或者是“石头数3倍、合并后1.98202倍”

总之N个石头合并后的体积应该是N的0.62271次方，不是N的0.63次方

因为当石头数增加至原来的3.04383倍时，合并后的体积分布与原来的分布是一样的，只是多了一个系数2

mathe你可以翻出你的原始数据看看，是不是这个结果

mathe

1. #include <math.h>
   
2. #include <time.h>
   
3. #include <stdio.h>
   
4. #define MAX(a,b) ((a)>(b))?(a):(b)
   
5. 
   
6. #define N 1000000
   
7. #define MD 3
   
8. #define LN (2*MD+1)
   
9. int base[N+1];
   
10. double matrix[N+1][LN];
    
11. double sum[N+1][LN];
    
12. 
    
13. int main()
    
14. {
    
15.         int h,s,a,b;
    
16.         time_t start;
    
17.         start=time(NULL);
    
18.         matrix[1][MD]=1.0;
    
19.         for(h=MD;h<LN;h++)sum[1][h]=1.0;
    
20.         base[1]=1;base[2]=2;
    
21.         matrix[2][MD]=1.0;
    
22.         for(h=MD;h<LN;h++)sum[2][h]=1.0;
    
23.         for(h=3;h<=N;h++){
    
24.                 base[h]=round(log(h)/log(3))+1;
    
25.                 for(a=0;a<LN;a++){
    
26.                         int hind = a+base[h]-MD;
    
27.                         double r=0.0;
    
28.                         for(s=1;s<=h-1;s++){
    
29.                                 int ua = hind+MD-base[s];
    
30.                                 if(ua<0||ua>=LN)continue;
    
31.                                 int ind2=h-hind;
    
32.                                 int b=hind-base[h-s]+MD-1;
    
33.                                 if(b<0)continue;
    
34.                                 double s2=0;
    
35.                                 if(b>=LN)s2=1.0;else s2=sum[h-s][b];
    
36.                                 r+=2.0*matrix[s][ua]*s2/(h-1);
    
37.                                 if(ua>0&&b<LN){
    
38.                                     r+=matrix[s][ua-1]*matrix[h-s][b]/(h-1);
    
39.                                 }
    
40.                         }
    
41.                         matrix[h][a]=r;
    
42.                         if(a==0)sum[h][a]=r;
    
43.                         else sum[h][a]=sum[h][a-1]+r;
    
44.                 }
    
45.                 if(h%128==0){
    
46.                         fprintf(stderr,"h=%d time %ds\n",h,(int)(time(NULL)-start));
    
47.                 }
    
48.         }
    
49.         for(h=2;h<=N;h++){
    
50.                 printf("h=%d\n",h);
    
51.                 double p=0;
    
52.                 double m=1.0;
    
53.                 for(s=0;s<LN;s++){
    
54.                         int a=base[h]+s-MD;
    
55.                         if(matrix[h][s]==0)continue;
    
56.                         printf("\tm[%d]=%g\n",a,matrix[h][s]);
    
57.                 }
    
58.         }
    
59.         return 0;
    
60. }
    

复制代码

这个比较乱的代码应该能够在一天以内运行到1000000，不知道误差有多大，你们有空可以试验一下

KeyTo9_Fans

8# mathe 点评：
我不知道你是怎么判断出来的，我觉得我的数据量不够评估，结果就是用3#的公式推算的

就是运行一下上面的代码，看一下前后两个峰值对应的N，差了多少倍啊

例如：

当n=1时，合并后的体积等于1的概率达到最大值，为1.00000000000000000
当n=2,3时，合并后的体积等于2的概率达到最大值，为1.00000000000000000
当n=8时，合并后的体积等于4的概率达到最大值，为0.97142857142857142
当n=24时，合并后的体积等于8的概率达到最大值，为0.94474148406846015
当n=74时，合并后的体积等于16的概率达到最大值，为0.93570797290841390
当n=224时，合并后的体积等于32的概率达到最大值，为0.93287352346257579
当n=683时，合并后的体积等于64的概率达到最大值，为0.93192170797609386
当n=2078时，合并后的体积等于128的概率达到最大值，为0.93160994038313749
当n=6326时，合并后的体积等于256的概率达到最大值，为0.93150737895583391
当n=19255时，合并后的体积等于512的概率达到最大值，为0.93147368686003040
当n=58610时，合并后的体积等于1024的概率达到最大值，为0.93146261826311938
当n=178399时，合并后的体积等于2048的概率达到最大值，为0.93145898182323894
当n=543016时，合并后的体积等于4096的概率达到最大值，为0.93145778712761851
……

由此可知当n=0.858558*3.043829^k时，合并后的体积等于2^k的概率达到最大值，为0.93145720258744+0.31741467/n

因此n个体积为1的石子合并后的体积大约是n的0.622709次方