一道很棘手的函数拟合问题

wayne

你的图很特别
是用gnuplot画的吗

KeyTo9_Fans

#####
9#的贴子是在哈尔滨的香格里拉大酒店里发的。

做个记号，留作纪念。
#####

由于要参加ACM World Final的各项活动，所以睡前草草地编了一个C++程序跑起来就不管了。

我以为给她一个晚上的运行时间肯定足够了。

醒来后发现这个程序越跑越费劲，参数爬升了几级就爬不动了。

仔细一看才发现步长不合适了。

于是把步长改大几个数量级继续跑，然后我离开宾馆参加竞赛去了。

打完竞赛回到宾馆，发现参数又爬不动了。

于是再把步长改大，再跑一个晚上。

醒来后发现又不对劲了，参数爬升了几级又爬不动了。

然后才反应过来自己犯傻了，应该把步长设置为动态的！

修改完毕后再次运行，看到参数刷的一声至上云霄，瞬间登顶！

于是发贴汇报并自我检讨。

程序代码如下：

1. #include<cmath>
   
2. #include<cstdio>
   
3. #include<cstdlib>
   
4. 
   
5. double x[180],y[180],a,p1,p2,p3,p4,p5,d1,d2,d3,d4,d5,min,tmp,b1,b2,b3,b4,b5;
   
6. int i,n,t;
   
7. 
   
8. double f(double x)        //把x代入，计算f(x)的值
   
9. {
   
10.         return p1/(p2+exp(p3*log(x)))+p4*exp(p5*log(x));
    
11. }
    
12. 
    
13. double ss()        //把所有的x逐个代入，计算总偏差
    
14. {
    
15.         int i;
    
16.         double s,s2;
    
17.         s2=0;
    
18.         for(i=0;i<n;i++)
    
19.         {
    
20.                 s=f(x[i])-y[i];        //计算一个数据点的偏差
    
21.                 s2+=s*s;        //把偏差平方并累计
    
22.         }
    
23.         return s2;
    
24. }
    
25. 
    
26. int main()
    
27. {
    
28.         scanf("%d",&n);
    
29.         for(i=0;i<n;i++)
    
30.                 scanf("%lf%lf",&x[i],&y[i]);        //读入数据
    
31.         p1=1.3738626629060553e+032;        //设置初值
    
32.         p2=2.4369393036751976e+028;        //设置初值
    
33.         p3=16.719431561664958;        //设置初值
    
34.         p4=0.0012054475996989941;        //设置初值
    
35.         p5=3.0453779404292880;        //设置初值
    
36.         min=ss();        //计算偏差
    
37.         while(1)        //死循环（在适当的时候人工终止运行，修改代码和参数，然后重新运行）
    
38.         {
    
39.                 a=p1/1e12;        //设置相对步长
    
40.                 d1=(rand()-16384)*a;        //设置绝对步长
    
41.                 d2=(rand()-16384)*a/1e4;        //设置绝对步长
    
42.                 d3=rand()/1e13-1.6e-9;        //设置绝对步长
    
43.                 d4=rand()/1e16-1.6e-12;        //设置绝对步长
    
44.                 d5=rand()/1e14-1.6e-10;        //设置绝对步长
    
45.                 p1+=d1;        //根据步长调整参数
    
46.                 p2+=d2;        //根据步长调整参数
    
47.                 p3+=d3;        //根据步长调整参数
    
48.                 p4+=d4;        //根据步长调整参数
    
49.                 p5+=d5;        //根据步长调整参数
    
50.                 tmp=ss();        //计算调整参数后的偏差
    
51.                 if(tmp<min)        //如果偏差比原来小
    
52.                 {
    
53.                         min=tmp;        //更新当前偏差，保留调整后的参数
    
54. //                        t++;        //累计更新次数
    
55. //                        if(t>100)        //每100次更新输出一次
    
56. //                        {
    
57. //                                t=0;
    
58.                                 printf("min=%lf\np1=%.0lf;\np2=%.0lf;\np3=%lf;\np4=%.9lf;\np5=%lf;\n\n",min,p1,p2,p3,p4,p5);        //输出参数
    
59. //                        }
    
60.                 }
    
61.                 else        //如果偏差比原来大
    
62.                 {
    
63.                         p1-=d1;        //还原参数
    
64.                         p2-=d2;        //还原参数
    
65.                         p3-=d3;        //还原参数
    
66.                         p4-=d4;        //还原参数
    
67.                         p5-=d5;        //还原参数
    
68.                 }
    
69.         }
    
70.         return 0;
    
71. }

复制代码

第一次做高难度的拟合作业，完全没有经验。

所以难以写出一气呵成的代码。

上述代码是跑跑停停，停停改改，改改又跑跑，才得到最终结果的。

其中$p_1$和$p_2$是分别从几百万和几千开始跑的。

然后眼看着它们一步一步地往上爬，一直没有停过。

最后才反应过来它们可能会直上云霄。

所以要把步长设置为相对的，而不是绝对的。

果然登顶的时候它们出奇的大，得用科学计数法来表示。

最后把步上逐渐减小，最终达到最佳拟合效果。

图是用curxpt软件画的。

curxpt软件自带拟合功能，但你这个数据它拟合不出来。

所以我才编了一个C++程序自己来拟合的。

wayne

程序很清晰，我还以为你用到了某相关算法，有时间了我也试试用C++来做

forcal

原来高手们早就求解过1stOpt的优化题。
一方面，高手就是高手；另一方面，1stOpt的算法的确精悍。
1stOpt的9道优化题，Forcal目前只解出了7道（使用随机初值），包括这个题。

G-Spider

14# forcal
很有意思，有时间也来尝试一下。

forcal

14# forcal
很有意思，有时间也来尝试一下。
G-Spider 发表于 2011-2-24 22:37

恩，作为练手，或者爱好，都值得一试。

gogdizzy

6# wayne

这个是什么原理？为什么可以除以40和1000，那原来的公式还能满足吗？

wayne

17# gogdizzy
满足啊，那就变成了y/1000与x/40之间的函数关系了