实验一 用线性回归进行一元非线性模型的参数估计
【实验目的】
通过编写具有一定通用型的可转化为现行的医院非线性模型程序,学习线性回归这种参数估计方法,体会作为经验模型主要形式之一的拟合模型的建立方法。 【实验内容】
试用线性回归分析建立某无烟煤的粒度模型。其粒度试验结果见下表:
某无烟煤的粒度特性 25 13 6 3 0.5 筛孔/mm 50 正累积产率/R% 9.91 18.82 32.32 47.38 63.02 91.44 【实验要求】绘制详细流程图,编写程序并上机验证,最后撰写实验报告。 【实验流程】 g=g-n*e*e; 开始 l=l-n*e*f; m=m-n*f*f; e=e/n; 实验粒 输入样本容量n、f=f/n; 选择模型 度数据x[i]y[i]、h=sqrt(g*m); h=l/h; i=0 c=l/g; d=f-c*e; N i 【实验程序】 #inclede Int n, j,k; float x[],y[],z[],u[],v[],p[],e=0,f=0,g=0,l=0,m=0,a,b,c,d,q=0,sm; printf(“请输入样本容量n=”); scanf(“%f”,&n); printf(“请输入粒度试验结果\\n”); for(j=0;j scanf(“%f%f”,&x[j],&y[j]); } printf(“请选择模型: \\n1.y=a+bx\\n2.y=100e^(-ax^b)\\n3.y=ax^b\\n4.y=ae^(bx)\\n5.y=ae^(b/x)\\n6.y=a+lgx\\n7.y=1/(a+be^(-x))\\n8.y=ab^x\\n”); scanf(“%d”,&k); switch(k) /*选择模型*/ { case 1:for(i=0;i u[i]=x[i]; v[i]=y[i]; } break; case 2:for(i=0;i u[i]=x[i]; v[i]=log(log(100/y[i])); } break; case 3:for(i=0;i u[i]=log10(x[i]); v[i]=log10(y[i]); } break; case 4:for(i=0;i u[i]=x[i]; v[i]=log(y[i]); } break; case 5:for(i=0;i u[i]=1/x[i]; v[i]=log(y[i]); } break; case 6:for(i=0;i u[i]=log10(x[i]); v[i]=1/y[i]; } break; case 7:for(i=0;i u[i]=exp(-x[i]); v[i]=1/y[i]; } break; case 8:for(i=0;i u[i]=x[i]; v[i]=log(y[i]); } break; default:printf(“不合法的模型号!\\n”);break; } for(j=0;j e=e+u[j]; /*∑x*/ f=f+v[j]; /*∑y*/ g=g+u[i]*u[i]; l=l+u[i]*v[i]; m=m+v[i]*v[i]; } g=g-n*e*e; /*∑x^2*/ l=l-n*e*f; /*∑xy*/ m=m-n*f*f; /*∑y^2*/ e=e/n; f=f/n; h=sqrt(g*m); h=l/h; c=l/g; d=f-c*e; switch(k) { case 1: case 6: case 7:a=d; b=c;break; case 2: case 3: case 4: case 5:a=exp(d); b=c;break; case 8:a=exp(d);b=exp?;break; default:printf(“非法模型号!\\n”);break; } for(i=0;i switch(k) { case 1:p[i]=a+b*x[i];printf(“拟合结果为:y=%f+%fx\\n”,&a,&b);break; case 2:p[i]=100*exp(-a*pow(x[i],b)); printf(“拟合结果为:y=100e^(-%fx^%f)\\n”,a,b);break; case 3:p[i]=a*pow(x[i],b); printf(“拟合结果为:y=%fx^%f\\n”,a,b);break; case 4:p[i]=a*exp(b*x[i]); printf(“拟合结果为:y=t^(%f x)\\n”,a,b);break; case 5:p[i]=a*exp(b/x[i]); printf(“拟合结果为:.y=%f e^(%f/x)\\n”,a,b);break; case 6:p[i]=a+b*log10(x[i]); printf(“拟合结果为:y=%f+%flgx \\n”,a,b);break; case 7:p[i]=1/(a+b*exp(-x[i])); printf(“拟合结果为:y=1/(%f+t^(-x))\\n”,a,b);break case 8:p[i]=a*pow(b,x[i]); printf(“拟合结果为:y=%f%f^x \\n”,a,b);break; default:printf(“非法模型号!\\n”);break; } } for(i=0;i q=q+(y[i]-p[i])*(y[i]-p[i]); } sm=sqrt(q/(n-2)); } printf{“a=%f b=%f 标准差=%f\\n”,a,b,sm}; printf(“***********************************\\n); printf(“ x y y的估计 y-(y的估计”); for(j=0;j printf(“ %f %f %f %f”,x[j],y[j],p[j],y[j]-p[j]); } } 【实验结果】