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

wayne

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

KeyTo9_Fans

##### 9#的贴子是在哈尔滨的香格里拉大酒店里发的。 做个记号，留作纪念。 ##### 由于要参加ACM World Final的各项活动，所以睡前草草地编了一个C++程序跑起来就不管了。 我以为给她一个晚上的运行时间肯定足够了。 醒来后发现这个程序越跑越费劲，参数爬升了几级就爬不动了。 仔细一看才发现步长不合适了。 于是把步长改大几个数量级继续跑，然后我离开宾馆参加竞赛去了。 打完竞赛回到宾馆，发现参数又爬不动了。 于是再把步长改大，再跑一个晚上。 醒来后发现又不对劲了，参数爬升了几级又爬不动了。 然后才反应过来自己犯傻了，应该把步长设置为动态的！ 修改完毕后再次运行，看到参数刷的一声至上云霄，瞬间登顶！ 于是发贴汇报并自我检讨。 程序代码如下：

1. #include<cmath>
2. #include<cstdio>
3. #include<cstdlib>
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;
8. double f(double x) //把x代入，计算f(x)的值
9. {
10. return p1/(p2+exp(p3*log(x)))+p4*exp(p5*log(x));
11. }
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. }
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之间的函数关系了