NIST非线性拟合-----接力赛
原始数据以及 待拟合的 目标公式都在该链接里,http://www.itl.nist.gov/div898/strd/nls/nls_main.shtml
我们不妨来 接力攻破一下~~ 给楼主一个建议:能否一个一个地给出题目的数据及中文说明?即给出一个,解决一个,再给出一个。
像我这样英文水平的人就也可以进行尝试了。也能方便更多像我这样的人。
低难度组第1题: Misra1a
本帖最后由 wayne 于 2011-2-21 21:30 编辑参考链接:http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml
题目背景意义: 关于牙医研究中单分子的吸收率。
=================================
待拟合的公式: y = b1*(1-exp[-b2*x])
其中b1,b2是参数,自变量x是体积,因变量y是压强。共14组观测数据:
Data: y x
10.07E0 77.6E0
14.73E0 114.9E0
17.94E0 141.1E0
23.93E0 190.8E0
29.61E0 239.9E0
35.18E0 289.0E0
40.02E0 332.8E0
44.82E0 378.4E0
50.76E0 434.8E0
55.05E0 477.3E0
61.01E0 536.8E0
66.40E0 593.1E0
75.47E0 689.1E0
81.78E0 760.0E0
第一题Forcal求解:!using["fcopt","math"];
init(::Array,max)=
{
max=14,
Array=arrayinitns{max,2 :
" 10.07E0 77.6E0
14.73E0 114.9E0
17.94E0 141.1E0
23.93E0 190.8E0
29.61E0 239.9E0
35.18E0 289.0E0
40.02E0 332.8E0
44.82E0 378.4E0
50.76E0 434.8E0
55.05E0 477.3E0
61.01E0 536.8E0
66.40E0 593.1E0
75.47E0 689.1E0
81.78E0 760.0E0"
}.free()
};
f(b1, b2 : i,s,x,y,z : Array,max)=
{
s=0,i=0,(i<max).while{
y=Array, x=Array,
s=s+[ b1*(1-exp[-b2*x])-y]^2,
i++
},
s
};
Opt;结果(前面是最优参数,最后一个数是误差,Forcal结果基本都是这样的):
238.9421251876012 5.501564435330968e-004 0.1245513889445084
从Forcal的角度,难度比较小。 3# wayne
用Mathematica解:
1、直接从网络抓取数据,不用拷贝:data = Import["http://www.itl.nist.gov/div898/strd/nls/data/LINKS/DATA/Misra1a.dat", "Data"];aa = Cases
2、数据有理化xy = Cases, {x_, y_} -> {y/100, x/10}]3、拟合参数值:{10 b1, b2/100} /. FindFit), {b1, b2}, x, WorkingPrecision -> 20]===================================================
{238.94212917913832416, 0.00055015643180517008397} 即便没有用过Mathematica,也开眼界,长见识。Mathematica的优化能力在软件四大家中恐怕要居首位吧?个人浅陋的理解,虽然都没有用过。
该推出第二题了吧?速度稍快些也无妨,呵呵。 英文看起来还是比较费劲的,我来第二题,错了大家给指出哦。
低难度组第2题: Chwirut2
参考链接:http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml
=================================
待拟合的公式: y = exp(-b1*x)/(b2+b3*x)
其中b1,b2,b3是拟合参数。共54组观测数据:Data:y x
92.9000E0 0.500E0
57.1000E0 1.000E0
31.0500E0 1.750E0
11.5875E0 3.750E0
8.0250E0 5.750E0
63.6000E0 0.875E0
21.4000E0 2.250E0
14.2500E0 3.250E0
8.4750E0 5.250E0
63.8000E0 0.750E0
26.8000E0 1.750E0
16.4625E0 2.750E0
7.1250E0 4.750E0
67.3000E0 0.625E0
41.0000E0 1.250E0
21.1500E0 2.250E0
8.1750E0 4.250E0
81.5000E0 .500E0
13.1200E0 3.000E0
59.9000E0 .750E0
14.6200E0 3.000E0
32.9000E0 1.500E0
5.4400E0 6.000E0
12.5600E0 3.000E0
5.4400E0 6.000E0
32.0000E0 1.500E0
13.9500E0 3.000E0
75.8000E0 .500E0
20.0000E0 2.000E0
10.4200E0 4.000E0
59.5000E0 .750E0
21.6700E0 2.000E0
8.5500E0 5.000E0
62.0000E0 .750E0
20.2000E0 2.250E0
7.7600E0 3.750E0
3.7500E0 5.750E0
11.8100E0 3.000E0
54.7000E0 .750E0
23.7000E0 2.500E0
11.5500E0 4.000E0
61.3000E0 .750E0
17.7000E0 2.500E0
8.7400E0 4.000E0
59.2000E0 .750E0
16.3000E0 2.500E0
8.6200E0 4.000E0
81.0000E0 .500E0
4.8700E0 6.000E0
14.6200E0 3.000E0
81.7000E0 .500E0
17.1700E0 2.750E0
81.3000E0 .500E0
28.9000E0 1.750E0 第二题Forcal求解:!using["fcopt","math"];
init(::Array,max)=
{
max=54,
Array=arrayinitns{max,2 :
" 92.9000E0 0.500E0
57.1000E0 1.000E0
31.0500E0 1.750E0
11.5875E0 3.750E0
8.0250E0 5.750E0
63.6000E0 0.875E0
21.4000E0 2.250E0
14.2500E0 3.250E0
8.4750E0 5.250E0
63.8000E0 0.750E0
26.8000E0 1.750E0
16.4625E0 2.750E0
7.1250E0 4.750E0
67.3000E0 0.625E0
41.0000E0 1.250E0
21.1500E0 2.250E0
8.1750E0 4.250E0
81.5000E0 .500E0
13.1200E0 3.000E0
59.9000E0 .750E0
14.6200E0 3.000E0
32.9000E0 1.500E0
5.4400E0 6.000E0
12.5600E0 3.000E0
5.4400E0 6.000E0
32.0000E0 1.500E0
13.9500E0 3.000E0
75.8000E0 .500E0
20.0000E0 2.000E0
10.4200E0 4.000E0
59.5000E0 .750E0
21.6700E0 2.000E0
8.5500E0 5.000E0
62.0000E0 .750E0
20.2000E0 2.250E0
7.7600E0 3.750E0
3.7500E0 5.750E0
11.8100E0 3.000E0
54.7000E0 .750E0
23.7000E0 2.500E0
11.5500E0 4.000E0
61.3000E0 .750E0
17.7000E0 2.500E0
8.7400E0 4.000E0
59.2000E0 .750E0
16.3000E0 2.500E0
8.6200E0 4.000E0
81.0000E0 .500E0
4.8700E0 6.000E0
14.6200E0 3.000E0
81.7000E0 .500E0
17.1700E0 2.750E0
81.3000E0 .500E0
28.9000E0 1.750E0"
}.free()
};
f(b1, b2, b3 : i,s,x,y : Array,max)=
{
s=0,i=0,(i<max).while{
y=Array, x=Array,
s=s+[ exp(-b1*x)/(b2+b3*x)-y]^2,
i++
},
s
};
Opt;结果:
0.1665774524075732 5.165354127091879e-003 1.214996964994243e-002 513.0480294504323 低难度组第3题: Chwirut1
参考链接:http://www.itl.nist.gov/div898/strd/nls/data/chwirut1.shtml
=================================
待拟合的公式: exp[-b1*x]/(b2+b3*x)
其中b1,b2,b3是拟合参数。共214组观测数据:Data:y x
92.9000E0 0.5000E0
78.7000E0 0.6250E0
64.2000E0 0.7500E0
64.9000E0 0.8750E0
57.1000E0 1.0000E0
43.3000E0 1.2500E0
31.1000E0 1.7500E0
23.6000E0 2.2500E0
31.0500E0 1.7500E0
23.7750E0 2.2500E0
17.7375E0 2.7500E0
13.8000E0 3.2500E0
11.5875E0 3.7500E0
9.4125E0 4.2500E0
7.7250E0 4.7500E0
7.3500E0 5.2500E0
8.0250E0 5.7500E0
90.6000E0 0.5000E0
76.9000E0 0.6250E0
71.6000E0 0.7500E0
63.6000E0 0.8750E0
54.0000E0 1.0000E0
39.2000E0 1.2500E0
29.3000E0 1.7500E0
21.4000E0 2.2500E0
29.1750E0 1.7500E0
22.1250E0 2.2500E0
17.5125E0 2.7500E0
14.2500E0 3.2500E0
9.4500E0 3.7500E0
9.1500E0 4.2500E0
7.9125E0 4.7500E0
8.4750E0 5.2500E0
6.1125E0 5.7500E0
80.0000E0 0.5000E0
79.0000E0 0.6250E0
63.8000E0 0.7500E0
57.2000E0 0.8750E0
53.2000E0 1.0000E0
42.5000E0 1.2500E0
26.8000E0 1.7500E0
20.4000E0 2.2500E0
26.8500E0 1.7500E0
21.0000E0 2.2500E0
16.4625E0 2.7500E0
12.5250E0 3.2500E0
10.5375E0 3.7500E0
8.5875E0 4.2500E0
7.1250E0 4.7500E0
6.1125E0 5.2500E0
5.9625E0 5.7500E0
74.1000E0 0.5000E0
67.3000E0 0.6250E0
60.8000E0 0.7500E0
55.5000E0 0.8750E0
50.3000E0 1.0000E0
41.0000E0 1.2500E0
29.4000E0 1.7500E0
20.4000E0 2.2500E0
29.3625E0 1.7500E0
21.1500E0 2.2500E0
16.7625E0 2.7500E0
13.2000E0 3.2500E0
10.8750E0 3.7500E0
8.1750E0 4.2500E0
7.3500E0 4.7500E0
5.9625E0 5.2500E0
5.6250E0 5.7500E0
81.5000E0 .5000E0
62.4000E0 .7500E0
32.5000E0 1.5000E0
12.4100E0 3.0000E0
13.1200E0 3.0000E0
15.5600E0 3.0000E0
5.6300E0 6.0000E0
78.0000E0 .5000E0
59.9000E0 .7500E0
33.2000E0 1.5000E0
13.8400E0 3.0000E0
12.7500E0 3.0000E0
14.6200E0 3.0000E0
3.9400E0 6.0000E0
76.8000E0 .5000E0
61.0000E0 .7500E0
32.9000E0 1.5000E0
13.8700E0 3.0000E0
11.8100E0 3.0000E0
13.3100E0 3.0000E0
5.4400E0 6.0000E0
78.0000E0 .5000E0
63.5000E0 .7500E0
33.8000E0 1.5000E0
12.5600E0 3.0000E0
5.6300E0 6.0000E0
12.7500E0 3.0000E0
13.1200E0 3.0000E0
5.4400E0 6.0000E0
76.8000E0 .5000E0
60.0000E0 .7500E0
47.8000E0 1.0000E0
32.0000E0 1.5000E0
22.2000E0 2.0000E0
22.5700E0 2.0000E0
18.8200E0 2.5000E0
13.9500E0 3.0000E0
11.2500E0 4.0000E0
9.0000E0 5.0000E0
6.6700E0 6.0000E0
75.8000E0 .5000E0
62.0000E0 .7500E0
48.8000E0 1.0000E0
35.2000E0 1.5000E0
20.0000E0 2.0000E0
20.3200E0 2.0000E0
19.3100E0 2.5000E0
12.7500E0 3.0000E0
10.4200E0 4.0000E0
7.3100E0 5.0000E0
7.4200E0 6.0000E0
70.5000E0 .5000E0
59.5000E0 .7500E0
48.5000E0 1.0000E0
35.8000E0 1.5000E0
21.0000E0 2.0000E0
21.6700E0 2.0000E0
21.0000E0 2.5000E0
15.6400E0 3.0000E0
8.1700E0 4.0000E0
8.5500E0 5.0000E0
10.1200E0 6.0000E0
78.0000E0 .5000E0
66.0000E0 .6250E0
62.0000E0 .7500E0
58.0000E0 .8750E0
47.7000E0 1.0000E0
37.8000E0 1.2500E0
20.2000E0 2.2500E0
21.0700E0 2.2500E0
13.8700E0 2.7500E0
9.6700E0 3.2500E0
7.7600E0 3.7500E0
5.4400E0 4.2500E0
4.8700E0 4.7500E0
4.0100E0 5.2500E0
3.7500E0 5.7500E0
24.1900E0 3.0000E0
25.7600E0 3.0000E0
18.0700E0 3.0000E0
11.8100E0 3.0000E0
12.0700E0 3.0000E0
16.1200E0 3.0000E0
70.8000E0 .5000E0
54.7000E0 .7500E0
48.0000E0 1.0000E0
39.8000E0 1.5000E0
29.8000E0 2.0000E0
23.7000E0 2.5000E0
29.6200E0 2.0000E0
23.8100E0 2.5000E0
17.7000E0 3.0000E0
11.5500E0 4.0000E0
12.0700E0 5.0000E0
8.7400E0 6.0000E0
80.7000E0 .5000E0
61.3000E0 .7500E0
47.5000E0 1.0000E0
29.0000E0 1.5000E0
24.0000E0 2.0000E0
17.7000E0 2.5000E0
24.5600E0 2.0000E0
18.6700E0 2.5000E0
16.2400E0 3.0000E0
8.7400E0 4.0000E0
7.8700E0 5.0000E0
8.5100E0 6.0000E0
66.7000E0 .5000E0
59.2000E0 .7500E0
40.8000E0 1.0000E0
30.7000E0 1.5000E0
25.7000E0 2.0000E0
16.3000E0 2.5000E0
25.9900E0 2.0000E0
16.9500E0 2.5000E0
13.3500E0 3.0000E0
8.6200E0 4.0000E0
7.2000E0 5.0000E0
6.6400E0 6.0000E0
13.6900E0 3.0000E0
81.0000E0 .5000E0
64.5000E0 .7500E0
35.5000E0 1.5000E0
13.3100E0 3.0000E0
4.8700E0 6.0000E0
12.9400E0 3.0000E0
5.0600E0 6.0000E0
15.1900E0 3.0000E0
14.6200E0 3.0000E0
15.6400E0 3.0000E0
25.5000E0 1.7500E0
25.9500E0 1.7500E0
81.7000E0 .5000E0
61.6000E0 .7500E0
29.8000E0 1.7500E0
29.8100E0 1.7500E0
17.1700E0 2.7500E0
10.3900E0 3.7500E0
28.4000E0 1.7500E0
28.6900E0 1.7500E0
81.3000E0 .5000E0
60.9000E0 .7500E0
16.6500E0 2.7500E0
10.0500E0 3.7500E0
28.9000E0 1.7500E0
28.9500E0 1.7500E0 第三题Forcal求解:!using["fcopt","math"];
init(::Array,max)=
{
max=214,
Array=arrayinitns{max,2 :
" 92.9000E0 0.5000E0
//(此处略去10000字)
28.9500E0 1.7500E0"
}.free()
};
f(b1, b2, b3 : i,s,x,y : Array,max)=
{
s=0,i=0,(i<max).while{
y=Array, x=Array,
s=s+[ exp[-b1*x]/(b2+b3*x)-y]^2,
i++
},
s
};
Opt;结果:
0.1902757520340506 6.131347827474484e-003 1.053102614699184e-002 2384.477139692282