预测一个数列的极限

wayne

如果是拟合的话，最好还要给出各个参数的 置信区间，才有意义。

wayne

俺也拟合了一下，选择的模型是
$y=a+b/{c+x^d}$
已经到了拟合的最高境界，几乎和插值的曲线重合了：

wayne

上面的给出了拟合后的残差：
{-0.0000433106, 0.000332, -0.000493819, -0.0000886767, 0.0000499126, 0.00010533, 0.000119389, 0.000113152, 0.0000965446, 0.0000747279, 0.0000505189, 0.0000255045, 5.90665*10^-7, -0.0000237077, -0.0000471074, -0.0000694644, -0.0000907193, -0.000110864}

c,a,b,t的置信区间是
{{-0.22857, -0.205862}, {1.19589, 1.45906}, {0.592995, 0.593745}, {1.44409, 1.53653}}

我把上面的图再放大了一点，注意有两条线，一条是粉红色的，是内插曲线，一条是拟合曲线

KeyTo9_Fans

置信区间是怎么得出来的？

这个数列比较特殊，越靠后的数据越有说服力。

前面几项的值基本不需要理会，如何重点利用最后几项的值才是关键。

不知道能不能加权处理。

比如：后一项的权重是前一项的$2$倍。

wayne

参数估计这块我很久没涉及了，脑袋里装的还是那个正态分布【u-3s,u+3s】的模型。
不同的概率分布有不同的说法，好像大致就是95%的范围

此题感觉用外插更合适

wayne

其实，俺这方面也没经验

Buffalo

8# Buffalo
呵呵，你给的方法我那本书里面大篇幅介绍了
wayne 发表于 2010-8-1 15:36

大篇幅介绍你也没看，至少是没学会，跟没有一样

Buffalo

对这种单调趋于某个极限的数列，可以假设$A_n=A+\sum_{k=1}^{\infty}\frac{C_k}{n^k}$，如果能够想办法找出系数$C_k$，立刻就可以求出极限$A=A_n-\sum_{k=1}^{\infty}\frac{C_k}{n^k}$。实际做的时候不是求出系数$C_k$，而是想办法直接消掉$C_k$，一种方法如下：
用$n^m$乘以方程$A_n=A+\sum_{k=1}^{\infty}\frac{C_k}{n^k}$两边，得$n^m A_n=An^m+\sum_{k=1}^{m}C_k n^{m-k}+\sum_{k=m+1}^{\infty}\frac{C_k}{n^k}$，两边对n作m阶差分得$\sum_{i=0}^m \frac{m!}{i!(m-i)!}(-1)^i (n+i)^m A_{n+i}=A m!+\sum_{k=m+1}^{\infty}\frac{D_k}{n^k}$，也就是说，新数列$\sum_{i=0}^m \frac{1}{i!(m-i)!}(-1)^i (n+i)^m A_{n+i}$和原数列$A_n$趋于同一个极限，但是收敛的速度更快，至少是$O(n^{-m-1})$量级。
对于本题做到m=7、8就可以发现新数列出现比较明显的长周期（周期约为12，振幅衰减的）振荡现象，不再是单调的，如果不想花时间对付这一项的话就只能到此为止。
新数列的最后三项为0.592781, 0.592776, 0.592772，可以估计出原数列的极限约为0.59276~0.59278。

wayne

17# Buffalo

呵呵，一语中的啊

KeyTo9_Fans

18楼的方法看起来很不错，可是我怎么都学不会。

结果我还是用了自己的那套土方法。

我花了很大的成本，得到了$A_19$的值：

$A_19=0.5905853419876694$

拟合过程如下：

1. #include<cstdio>
   
2. #include<math.h>
   
3. 
   
4. double a[999999]={
   
5.         0,
   
6.         0.5,
   
7.         0.54119610014619698439972320536639,
   
8.         0.55929631600132353475576392078938,
   
9.         0.56972413396802715322804051835303,
   
10.         0.57581007321162765360503297431458,
    
11.         0.57970275713244352141943997833056,
    
12.         0.58235129508008298007347469183041,
    
13.         0.58424146648984767386035113239817,
    
14.         0.58564155686139651141699566635464,
    
15.         0.58671003405340635980469047312441,
    
16.         0.58754560137670707686509637674778,
    
17.         0.58821247060644326317197307974183,
    
18.         0.58875395365138276709785553207334,
    
19.         0.58920017119372364434464005947868,
    
20.         0.58957262470961612944817110769224,
    
21.         0.5898870137232585,
    
22.         0.5901550299615904,
    
23.         0.5903855332606833,
    
24.         0.5905853419876694
    
25. };
    
26. 
    
27. int i,n;
    
28. double b[999999],c[999999],d[999999],e[999999],f[999999],g[999999];
    
29. 
    
30. double f_(int i)
    
31. {
    
32.         return -1.4976171+0.26715566*i;
    
33. }
    
34. 
    
35. int main()
    
36. {
    
37.         n=20;
    
38.         for(i=1;i<n;i++)
    
39.         {
    
40.                 b[i]=exp(-log(a[i]-a[i-1])/2.75);
    
41.         }
    
42.         for(i=2;i<n;i++)
    
43.         {
    
44.                 c[i]=b[i]-b[i-1];
    
45.         }
    
46.         for(i=3;i<n;i++)
    
47.         {
    
48.                 d[i]=exp(-log(c[i-1]-c[i])/2.712);
    
49.         }
    
50.         for(i=4;i<n;i++)
    
51.         {
    
52.                 e[i]=d[i]-d[i-1];
    
53.         }
    
54.         for(i=5;i<n;i++)
    
55.         {
    
56.                 f[i]=exp(-log(e[i]-e[i-1])/3.853369);
    
57.                 printf("%d %.16lf %.16lf\n",i,f[i],f_(i));
    
58.         }
    
59.         for(i=n;i<999999;i++)
    
60.                 f[i]=f_(i);
    
61.         for(i=n;i<999999;i++)
    
62.                 f[i]=f[i-1]+g[i];
    
63.         for(i=n;i<999999;i++)
    
64.                 e[i]=e[i-1]+exp(-log(f[i])*3.853369);
    
65.         for(i=n;i<999999;i++)
    
66.                 d[i]=d[i-1]+e[i];
    
67.         for(i=n;i<999999;i++)
    
68.                 c[i]=c[i-1]-exp(-log(d[i])*2.712);
    
69.         for(i=n;i<999999;i++)
    
70.                 b[i]=b[i-1]+c[i];
    
71.         for(i=n;i<999999;i++)
    
72.                 a[i]=a[i-1]+exp(-log(b[i])*2.75);
    
73. 
    
74.         for(i=0;i<999999;i+=(i>31?i:1))
    
75.                 printf("%6d: %.16lf\n",i,a[i]);
    
76.         return 0;
    
77. }

复制代码

程序输出如下：

1. 5 -1.#IND000000000000 -0.1618388000000000
   
2. 6 -1.#IND000000000000 0.1053168600000001
   
3. 7 -1.#IND000000000000 0.3724725200000001
   
4. 8 -1.#IND000000000000 0.6396281800000001
   
5. 9 -1.#IND000000000000 0.9067838400000003
   
6. 10 1.#QNAN00000000000 1.1739395000000001
   
7. 11 1.2400299956402463 1.4410951599999999
   
8. 12 1.7184718045516196 1.7082508200000002
   
9. 13 1.9384902061929288 1.9754064800000004
   
10. 14 2.2365949301518984 2.2425621400000004
    
11. 15 2.5061186902492905 2.5097177999999998
    
12. 16 2.7768734284468222 2.7768734600000000
    
13. 17 3.0440290793828071 3.0440291200000003
    
14. 18 3.3111847500498381 3.3111847800000005
    
15. 19 3.5783404007816051 3.5783404399999998
    
16.      0: 0.0000000000000000
    
17.      1: 0.5000000000000000
    
18.      2: 0.5411961001461970
    
19.      3: 0.5592963160013236
    
20.      4: 0.5697241339680271
    
21.      5: 0.5758100732116277
    
22.      6: 0.5797027571324436
    
23.      7: 0.5823512950800830
    
24.      8: 0.5842414664898477
    
25.      9: 0.5856415568613965
    
26.     10: 0.5867100340534064
    
27.     11: 0.5875456013767071
    
28.     12: 0.5882124706064432
    
29.     13: 0.5887539536513827
    
30.     14: 0.5892001711937237
    
31.     15: 0.5895726247096161
    
32.     16: 0.5898870137232585
    
33.     17: 0.5901550299615904
    
34.     18: 0.5903855332606833
    
35.     19: 0.5905853419876694
    
36.     20: 0.5907597763774112
    
37.     21: 0.5909130395624155
    
38.     22: 0.5910484896187661
    
39.     23: 0.5911688369727856
    
40.     24: 0.5912762897683751
    
41.     25: 0.5913726623579062
    
42.     26: 0.5914594572711018
    
43.     27: 0.5915379278478790
    
44.     28: 0.5916091265966935
    
45.     29: 0.5916739428925802
    
46.     30: 0.5917331326283863
    
47.     31: 0.5917873417313245
    
48.     32: 0.5918371249591274
    
49.     64: 0.5924660692652847
    
50.    128: 0.5926612272674896
    
51.    256: 0.5927205350417009
    
52.    512: 0.5927383640078597
    
53.   1024: 0.5927436941980384
    
54.   2048: 0.5927452833002713
    
55.   4096: 0.5927457564031907
    
56.   8192: 0.5927458971557772
    
57. 16384: 0.5927459390163858
    
58. 32768: 0.5927459514637898
    
59. 65536: 0.5927459551647473
    
60. 131072: 0.5927459562650940
    
61. 262144: 0.5927459565922527
    
62. 524288: 0.5927459566891469

复制代码

结论：极限约为$0.59274596$，但是无法给出误差范围。

我接下来会给出详细步骤。

希望wayne等大牛仔细看，帮我估计一下最终结果的误差有多大。