是用gnuplot画的吗 #####
9#的贴子是在哈尔滨的香格里拉大酒店里发的。
做个记号，留作纪念。
#####
由于要参加ACM World Final的各项活动，所以睡前草草地编了一个C++程序跑起来就不管了。
我以为给她一个晚上的运行时间肯定足够了。
醒来后发现这个程序越跑越费劲，参数爬升了几级就爬不动了。
仔细一看才发现步长不合适了。
于是把步长改大几个数量级继续跑，然后我离开宾馆参加竞赛去了。
打完竞赛回到宾馆，发现参数又爬不动了。
于是再把步长改大，再跑一个晚上。
醒来后发现又不对劲了，参数爬升了几级又爬不动了。
然后才反应过来自己犯傻了，应该把步长设置为动态的！
修改完毕后再次运行，看到参数刷的一声至上云霄，瞬间登顶！
于是发贴汇报并自我检讨。
程序代码如下：#include<cmath>
#include<cstdio>
#include<cstdlib>
double x,y,a,p1,p2,p3,p4,p5,d1,d2,d3,d4,d5,min,tmp,b1,b2,b3,b4,b5;
int i,n,t;
double f(double x) //把x代入，计算f(x)的值
{
return p1/(p2+exp(p3*log(x)))+p4*exp(p5*log(x));
}
double ss() //把所有的x逐个代入，计算总偏差
{
int i;
double s,s2;
s2=0;
for(i=0;i<n;i++)
{
s=f(x)-y; //计算一个数据点的偏差
s2+=s*s; //把偏差平方并累计
}
return s2;
}
int main()
{
scanf("%d",&n);
for(i=0;i<n;i++)
scanf("%lf%lf",&x,&y); //读入数据
p1=1.3738626629060553e+032; //设置初值
p2=2.4369393036751976e+028; //设置初值
p3=16.719431561664958; //设置初值
p4=0.0012054475996989941; //设置初值
p5=3.0453779404292880; //设置初值
min=ss(); //计算偏差
while(1) //死循环（在适当的时候人工终止运行，修改代码和参数，然后重新运行）
{
a=p1/1e12; //设置相对步长
d1=(rand()-16384)*a; //设置绝对步长
d2=(rand()-16384)*a/1e4; //设置绝对步长
d3=rand()/1e13-1.6e-9; //设置绝对步长
d4=rand()/1e16-1.6e-12; //设置绝对步长
d5=rand()/1e14-1.6e-10; //设置绝对步长
p1+=d1; //根据步长调整参数
p2+=d2; //根据步长调整参数
p3+=d3; //根据步长调整参数
p4+=d4; //根据步长调整参数
p5+=d5; //根据步长调整参数
tmp=ss(); //计算调整参数后的偏差
if(tmp<min) //如果偏差比原来小
{
min=tmp; //更新当前偏差，保留调整后的参数
// t++; //累计更新次数
// if(t>100) //每100次更新输出一次
// {
// t=0;
printf("min=%lf\np1=%.0lf;\np2=%.0lf;\np3=%lf;\np4=%.9lf;\np5=%lf;\n\n",min,p1,p2,p3,p4,p5); //输出参数
// }
}
else //如果偏差比原来大
{
p1-=d1; //还原参数
p2-=d2; //还原参数
p3-=d3; //还原参数
p4-=d4; //还原参数
p5-=d5; //还原参数
}
}
return 0;
}第一次做高难度的拟合作业，完全没有经验。
所以难以写出一气呵成的代码。
上述代码是跑跑停停，停停改改，改改又跑跑，才得到最终结果的。
其中$p_1$和$p_2$是分别从几百万和几千开始跑的。
然后眼看着它们一步一步地往上爬，一直没有停过。
最后才反应过来它们可能会直上云霄。
所以要把步长设置为相对的，而不是绝对的。
果然登顶的时候它们出奇的大，得用科学计数法来表示。
最后把步上逐渐减小，最终达到最佳拟合效果。
图是用curxpt软件画的。
curxpt软件自带拟合功能，但你这个数据它拟合不出来。
所以我才编了一个C++程序自己来拟合的。 程序很清晰，我还以为你用到了某相关算法，有时间了我也试试用C++来做 原来高手们早就求解过1stOpt的优化题。
一方面，高手就是高手；另一方面，1stOpt的算法的确精悍。
1stOpt的9道优化题，Forcal目前只解出了7道（使用随机初值），包括这个题。 14# forcal
很有意思，有时间也来尝试一下。 14# forcal
很有意思，有时间也来尝试一下。
G-Spider 发表于 2011-2-24 22:37 http://bbs.emath.ac.cn/images/common/back.gif
恩，作为练手，或者爱好，都值得一试。 6# wayne
这个是什么原理？为什么可以除以40和1000，那原来的公式还能满足吗？ 17# gogdizzy
满足啊，那就变成了y/1000与x/40之间的函数关系了
页:
1
[2]