预测一个数列的极限

wayne · 发表于 2010-8-3 09:07:11

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

wayne · 发表于 2010-8-3 09:19:14

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

wayne · 发表于 2010-8-3 09:24:36

上面的给出了拟合后的残差： {-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 · 发表于 2010-8-4 16:59:25

置信区间是怎么得出来的？这个数列比较特殊，越靠后的数据越有说服力。前面几项的值基本不需要理会，如何重点利用最后几项的值才是关键。不知道能不能加权处理。比如：后一项的权重是前一项的$2$倍。

wayne · 发表于 2010-8-4 17:34:27

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

wayne · 发表于 2010-8-4 17:37:51

其实，俺这方面也没经验

Buffalo · 发表于 2010-8-5 05:08:36

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

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

Buffalo · 发表于 2010-8-5 05:25:13

对这种单调趋于某个极限的数列，可以假设$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 · 发表于 2010-8-5 10:12:03

17# Buffalo 呵呵，一语中的啊

KeyTo9_Fans · 发表于 2010-9-16 15:19:13

18楼的方法看起来很不错，可是我怎么都学不会。结果我还是用了自己的那套土方法。我花了很大的成本，得到了$A_19$的值： $A_19=0.5905853419876694$ 拟合过程如下：

#include<cstdio>
#include<math.h>
double a[999999]={
0,
0.5,
0.54119610014619698439972320536639,
0.55929631600132353475576392078938,
0.56972413396802715322804051835303,
0.57581007321162765360503297431458,
0.57970275713244352141943997833056,
0.58235129508008298007347469183041,
0.58424146648984767386035113239817,
0.58564155686139651141699566635464,
0.58671003405340635980469047312441,
0.58754560137670707686509637674778,
0.58821247060644326317197307974183,
0.58875395365138276709785553207334,
0.58920017119372364434464005947868,
0.58957262470961612944817110769224,
0.5898870137232585,
0.5901550299615904,
0.5903855332606833,
0.5905853419876694
};
int i,n;
double b[999999],c[999999],d[999999],e[999999],f[999999],g[999999];
double f_(int i)
{
return -1.4976171+0.26715566*i;
}
int main()
{
n=20;
for(i=1;i<n;i++)
{
b[i]=exp(-log(a[i]-a[i-1])/2.75);
}
for(i=2;i<n;i++)
{
c[i]=b[i]-b[i-1];
}
for(i=3;i<n;i++)
{
d[i]=exp(-log(c[i-1]-c[i])/2.712);
}
for(i=4;i<n;i++)
{
e[i]=d[i]-d[i-1];
}
for(i=5;i<n;i++)
{
f[i]=exp(-log(e[i]-e[i-1])/3.853369);
printf("%d %.16lf %.16lf\n",i,f[i],f_(i));
}
for(i=n;i<999999;i++)
f[i]=f_(i);
for(i=n;i<999999;i++)
f[i]=f[i-1]+g[i];
for(i=n;i<999999;i++)
e[i]=e[i-1]+exp(-log(f[i])*3.853369);
for(i=n;i<999999;i++)
d[i]=d[i-1]+e[i];
for(i=n;i<999999;i++)
c[i]=c[i-1]-exp(-log(d[i])*2.712);
for(i=n;i<999999;i++)
b[i]=b[i-1]+c[i];
for(i=n;i<999999;i++)
a[i]=a[i-1]+exp(-log(b[i])*2.75);
for(i=0;i<999999;i+=(i>31?i:1))
printf("%6d: %.16lf\n",i,a[i]);
return 0;
}

复制代码

程序输出如下：

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

复制代码

结论：极限约为$0.59274596$，但是无法给出误差范围。我接下来会给出详细步骤。希望wayne等大牛仔细看，帮我估计一下最终结果的误差有多大。

账号		自动登录	找回密码
密码			欢迎注册

[讨论] 预测一个数列的极限

点评