最小二乘法程序
最小二乘法程序
#include
#include
#include
#include
#define N 5//N个点
#define T 3 //T次拟合
#define W 1//权函数
#define PRECISION 0.00001
float pow_n(float a,int n)
{
int i;
if(n==0)
return(1);
float res=a;
for(i=1;i { res*=a; } return(res); } void mutiple(float a[][N],float b[][T+1],float c[][T+1]) { float res=0; int i,j,k; for(i=0;i for(j=0;j { res=0; for(k=0;k { res+=a[i][k]*b[k][j]; c[i][j]=res; } } } void matrix_trans(float a[][T+1],float b[][N]) { int i,j; for(i=0;i { for(j=0;j { b[j][i]=a[i][j]; } } } void init(float x_y[][2],int n) { int i; printf("请输入%d个已知点：\n",N); for(i=0;i { printf("(x%d y%d):",i,i); scanf("%f %f",&x_y[i][0],&x_y[i][1]); } } void get_A(float matrix_A[][T+1],float x_y[][2],int n) { int i,j; for(i=0;i { for(j=0;j { matrix_A[i][j]=W*pow_n(x_y[i][0],j); } } } void print_array(float array[][T+1],int n) { int i,j; for(i=0;i { for(j=0;j { printf("%-g",array[i][j]); } printf("\n"); } } void convert(float argu[][T+2],int n) { int i,j,k,p,t; float rate,temp; for(i=1;i { for(j=i;j { if(argu[i-1][i-1]==0) { for(p=i;p { if(argu[p][i-1]!=0) break; } if(p==n) { printf("方程组无解!\n"); exit(0); } for(t=0;t { temp=argu[i-1][t]; argu[i-1][t]=argu[p][t]; argu[p][t]=temp; } } rate=argu[j][i-1]/argu[i-1][i-1]; for(k=i-1;k { argu[j][k]-=argu[i-1][k]*rate; if(fabs(argu[j][k])<=PRECISION) argu[j][k]=0; } } } } void compute(float argu[][T+2],int n,float root[]) { int i,j; float temp; for(i=n-1;i>=0;i--) { temp=argu[i][n]; for(j=n-1;j>i;j--) { temp-=argu[i][j]*root[j]; } root[i]=temp/argu[i][i]; } } void get_y(float trans_A[][N],float x_y[][2],float y[],int n) { int i,j; float temp; for(i=0;i { temp=0; for(j=0;j { temp+=trans_A[i][j]*x_y[j][1]; } y[i]=temp; } } void cons_formula(float coef_A[][T+1],float y[],float coef_form[][T+2]) { int i,j; for(i=0;i { for(j=0;j { if(j==T+1) coef_form[i][j]=y[i]; else coef_form[i][j]=coef_A[i][j]; } } } void print_root(float a[],int n) { int i,j; printf("%d个点的%d次拟合的多项式系数为:\n",N,T); for(i=0;i { printf("a[%d]=%g,",i+1,a[i]); } printf("\n"); printf("拟合曲线方程为:\ny(x)=%g",a[0]); for(i=1;i { printf(" + %g",a[i]); for(j=0;j { printf("*X"); } } printf("\n"); } void process() { float x_y[N][2],matrix_A[N][T+1],trans_A[T+1][N],coef_A[T+1][T+1],coef_formu[T+1 ][T+2],y[T+1],a[T+1]; init(x_y,N); get_A(matrix_A,x_y,N); printf("矩阵A为：\n"); print_array(matrix_A,N); matrix_trans(matrix_A,trans_A); mutiple(trans_A,matrix_A,coef_A); printf("法矩阵为：\n"); print_array(coef_A,T+1); get_y(trans_A,x_y,y,T+1); cons_formula(coef_A,y,coef_formu); convert(coef_formu,T+1); compute(coef_formu,T+1,a); print_root(a,T+1); } void main() { process(); } ]]> 你可以改一下 不从终端输入,直接在程序中给出参数 请输入5个已知点： (x0 y0):-2 -0.1 (x1 y1):-1 0.1 (x2 y2):0 0.4 (x3 y3):1 0.9 (x4 y4):2 1.6 矩阵A为： 1 - 2 4 -8 1 -1 1 -1 1 0 0 0 1 1 1 1 1 2 4 8 法矩阵为： 5 0 10 0 0 10 0 34 10 0 34 0 0 34 0 130 5个点的3次拟合的多项式系数为: a[1]=0.408571, a[2]=0.391667, a[3]=0.0857143, a[4]=0.00833333, 拟合曲线方程为: y(x)=0.408571 + 0.391667*X + 0.0857143*X*X + 0.00833333*X*X*X ]]> #include void main () { int num,i; float x,y,l,m,n,p,a,b; i=1; l=0.0; m=0.0; n=0.0; p=0.0; printf ("请输入你想计算的x，y的个数："); scanf("%d",&num); if (num>=1) { while (i<=num); { printf("请输入x的值"); scanf ("%lf",&x); printf("请输入y的值"); scanf ("%lf",&y); l+=x; m+=y; n+=x*y; p+=x*x; i++; } a=(num*n-l*m)/(num*p-l*l); b=(p*m-n*l)/(num*p-l*l); printf("最小二乘法所算得的斜率和截距分别为%f和%f\n",a,b); } else printf("mun"输入有误！); } #include #include //////////////////////////////////////////////////////////////////// ////////////////////// //矩阵结构体 struct Matrix { int m,n;//m为行数,n为列数 double **pm;//指向矩阵二维数组的指针 }; //初始化矩阵mt，并置矩阵的行为m，列为n void InitMatrix(struct Matrix *mt,int m,int n) { int i; //指定矩阵的行和列 mt->m=m; mt->n=n; //为矩阵分配内存 mt->pm=new double *[m]; for(i=0;i { mt->pm[i]=new double[n]; } } //用0初始化矩阵mt，并置矩阵的行为m，列为n void InitMatrix0(struct Matrix *mt,int m,int n) { int i,j; InitMatrix(mt,m,n); for(i=0;i for(j=0;j } //销毁矩阵mt void DestroyMatrix(struct Matrix *mt) { int i; //释放矩阵内存 for(i=0;i { delete []mt->pm[i]; } delete []mt->pm; } //矩阵相乘mt3=mt1*mt2 //成功返回1，失败返回0 int MatrixMul(Matrix *mt1,Matrix *mt2,Matrix *mt3) { int i,j,k; double s; //检查行列号是否匹配 if(mt1->n!=mt2->m||mt1->m!=mt3->m||mt2->n!=mt3->n) return 0; for(i=0;i for(j=0;j { s=0.0; for(k=0;k mt3->pm[i][j]=s; } return 1; } //矩阵转置mt2=T(mt1) //成功返回1，失败返回0 int MatrixTran(Matrix *mt1,Matrix *mt2) { int i,j; //检查行列号是否匹配 if(mt1->m!=mt2->n||mt1->n!=mt2->m) return 0; for(i=0;i for(j=0;j return 1; } //控制台显示矩阵mt void ConsoleShowMatrix(Matrix *mt) { int i,j; for(i=0;i { printf("\n"); for(j=0;j } printf("\n"); } //////////////////////////////////////////////////////////////////// ////////////////////// //Guass列主元素消元法求解方程Ax=b,a=(A,b) int Guassl(double **a,double *x,int n) { int i,j,k,numl,*h,t; double *l,max; l=new double[n]; h=new int[n]; for(i=0;i for(i=1;i { max=fabs(a[h[i-1]][i-1]); numl=i-1; //列元的最大值 for(j=i;j { if(fabs(a[h[j]][i-1])>max) { numl=h[j]; max=fabs(a[h[j]][i-1]); } } if(max<0.00000000001) return 0; //交换行号 if(numl>i-1) { t=h[i]; h[i]=h[numl]; h[numl]=t; } for(j=i;j for(k=i;k a[h[j]][k]=a[h[j]][k]-l[j]*a[h[i-1]][k]; } for(i=n-1;i>=0;i--) { x[i]=a[h[i]][n]; for(j=i+1;j } //清除临时数组内存 delete []l; delete []h; return 1; } //////////////////////////////////////////////////////////////////// /////////////////////// //最小二乘法求解矩阵Ax=b int MinMul(Matrix A,Matrix b,double *x) { int i,j; if(b.n!=1) return 0; if(A.m!=b.m) return 0; Matrix TranA;//定义A的转置 InitMatrix0(&TranA,A.n,A.m); MatrixTran(&A,&TranA); Matrix TranA_A;//定义A的转置与A的乘积矩阵 InitMatrix0(&TranA_A,A.n,A.n); MatrixMul(&TranA,&A,&TranA_A);//A的转置与A的乘积 Matrix TranA_b;//定义A的转置与b的乘积矩阵 InitMatrix0(&TranA_b,A.n,1); MatrixMul(&TranA,&b,&TranA_b);//A的转置与b的乘积 DestroyMatrix(&TranA);//释放A的转置的内存 Matrix TranA_A_b;//定义增广矩阵 InitMatrix0(&TranA_A_b,TranA_A.m,TranA_A.m+1); //增广矩阵赋值 for(i=0;i { for(j=0;j } DestroyMatrix(&TranA_A); DestroyMatrix(&TranA_b); //Guass列主消元法求解 Guassl(TranA_A_b.pm,x,TranA_A_b.m); DestroyMatrix(&TranA_A_b); return 1; } int MinMul(double **A,double *b,int m,int n,double *x) { int r,i; Matrix Al,bl; Al.pm=new double *[m]; Al.m=m; Al.n=n; InitMatrix(&bl,m,1); for(i=0;i { Al.pm[i]=A[i]; bl.pm[i][0]=b[i]; } r=MinMul(Al,bl,x); delete []Al.pm; DestroyMatrix(&bl); return r; } 曲线拟合（curve-fitting ）：工程实践中，用测量到的一些离散的数据},...2,1,0),,{(m i y x i i =求一个近似的函数)(x ?来拟合这组数据，要求所得的拟合曲线能最好的反映数据的基本趋势（即使)(x ?最好地逼近()x f ，而不必满足插值原则。因此没必要取)(i x ?=i y ，只要使i i i y x -=)(?δ尽可能地小）。 原理： 给定数据点},...2,1,0),,{(m i y x i i =。求近似曲线)(x ?。并且使得近似曲线与()x f 的偏差最小。近似曲线在该点处的偏差i i i y x -=)(?δ，i=1,2,...,m 。 常见的曲线拟合方法: 1.使偏差绝对值之和最小 2.使偏差绝对值最大的最小 3.使偏差平方和最小 最小二乘法： 按偏差平方和最小的原则选取拟合曲线，并且采取二项式方程为拟合曲线的方法,称为最小二乘法。 推导过程： 1. 设拟合多项式为： 2. 各点到这条曲线的距离之和，即偏差平方和如下： 3. 问题转化为求待定系数0a ...k a 对等式右边求i a 偏导数，因而我们得到 了： ....... 4、 把这些等式化简并表示成矩阵的形式，就可以得到下面的矩阵： 5. 将这个范德蒙得矩阵化简后可得到: 6. 也就是说X*A=Y，那么A = (X'*X)-1*X'*Y，便得到了系数矩阵A，同时，我们也就得到了拟合曲线。 MATLAB实现： MATLAB提供了polyfit（）函数命令进行最小二乘曲线拟合。 调用格式：p=polyfit(x,y,n) [p,s]= polyfit(x,y,n) [p,s,mu]=polyfit(x,y,n) x,y为数据点，n为多项式阶数，返回p为幂次从高到低的多项式系数向量p。x必须是单调的。矩阵s包括R（对x进行QR分解的三角元素）、df(自由度）、normr(残差）用于生成预测值的误差估计。 [p,s,mu]=polyfit(x,y,n)在拟合过程中，首先对x进行数据标准化处理，以在拟合中消除量纲等影响，mu包含标准化处理过程中使用的x的均值和标准差。 polyval( )为多项式曲线求值函数，调用格式： y=polyval(p,x) [y,DELTA]=polyval(p,x,s) y=polyval(p,x)为返回对应自变量x在给定系数P的多项式的值。 [y,DELTA]=polyval(p,x,s) 使用polyfit函数的选项输出s得出误差估计Y DELTA。它假设polyfit函数数据输入的误差是独立正态的，并且方差为常数。则Y DELTA将至少包含50%的预测值。 如下给定数据的拟合曲线： x=[0.5,1.0,1.5,2.0,2.5,3.0]， y=[1.75,2.45,3.81,4.80,7.00,8.60]。 解：MATLAB程序如下： x=[0.5,1.0,1.5,2.0,2.5,3.0]; y=[1.75,2.45,3.81,4.80,7.00,8.60]; p=polyfit(x,y,2) x1=0.5:0.05:3.0; y1=polyval(p,x1); plot(x,y,'*r',x1,y1,'-b') 运行结果如图1 计算结果为： p =0.5614 0.8287 1.1560 即所得多项式为y=0.5614x^2+0.08287x+1.15560 图1 最小二乘法曲线拟合示例 对比检验拟合的有效性： 例：在[0,π]区间上对正弦函数进行拟合，然后在[0,2π]区间画出图形，比较拟合区间和非拟合区间的图形，考察拟合的有效性。 在MATLAB中输入如下代码： clear x=0:0.1:pi; y=sin(x); [p,mu]=polyfit(x,y,9) 最小二乘法及其应用 1． 引言 最小二乘法在19世纪初发明后,很快得到欧洲一些国家的天文学家和测地学家的广泛关注。据不完全统计,自1805年至1864年的60年间,有关最小二乘法的研究论文达256篇,一些百科全书包括1837年出版的大不列颠百科全书第7版,亦收入有关方法的介绍。同时,误差的分布是“正态”的,也立刻得到天文学家的关注及大量经验的支持。如贝塞尔( F. W. Bessel, 1784—1846)对几百颗星球作了三组观测,并比较了按照正态规律在给定范围内的理论误差值和实际值,对比表明它们非常接近一致。拉普拉斯在1810年也给出了正态规律的一个新的理论推导并写入其《分析概论》中。正态分布作为一种统计模型,在19世纪极为流行,一些学者甚至把19世纪的数理统计学称为正态分布的统治时代。在其影响下,最小二乘法也脱出测量数据意义之外而发展成为一个包罗极大,应用及其广泛的统计模型。到20世纪正态小样本理论充分发展后,高斯研究成果的影响更加显著。最小二乘法不仅是19世纪最重要的统计方法,而且还可以称为数理统计学之灵魂。相关回归分析、方差分析和线性模型理论等数理统计学的几大分支都以最小二乘法为理论基础。正如美国统计学家斯蒂格勒( S. M. Stigler)所说,“最小二乘法之于数理统计学犹如微积分之于数学”。最小二乘法是参数回归的最基本得方法所以研究最小二乘法原理及其应用对于统计的学习有很重要的意义。 2. 最小二乘法 所谓最小二乘法就是：选择参数10,b b ,使得全部观测的残差平方和最小. 用数学公式表示为： 21022)()(m in i i i i i x b b Y Y Y e --=-=∑∑∑∧ 为了说明这个方法，先解释一下最小二乘原理，以一元线性回归方程为例. i i i x B B Y μ++=10 （一元线性回归方程） "; for($j = 1;$j <= 9;$j++){ echo " } for($i = 1;$i <= $j;$i++){ echo " 用matlab实现最小二乘递推算法辨识系统参 数 自动化系统仿真实验室指导教师： 学生姓名班级计082-2 班学号撰写时间： 全文结束》》-3-1 成绩评定： 一．设计目的 1、学会用Matlab实现最小二乘法辨识系统参数。 2、进一步熟悉Matlab的界面及基本操作； 3、了解并掌握Matlab中一些函数的作用与使用；二．设计要求最小二乘递推算法辨识系统参数，利用matlab编程实现，设初始参数为零。z(k)-1、5*z(k-1)+0、7*z(k-2)=1*u(k-1)+0、5*u(k-2)+v(k); 选择如下形式的辨识模型：z(k)+a1*z(k- 1)+a2*z(k-2)=b1*u(k-1)+b2*u(k-2)+v(k);三．实验程序 m=3;N=100;uk=rand(1,N);for i=1:Nuk(i)=uk(i)*(-1)^(i-1);endyk=zeros(1,N); for k=3:N yk(k)=1、5*yk(k-1)-0、 7*yk(k-2)+uk(k-1)+0、5*uk(k-2); end%j=100;kn=0;%y=yk(m:j);%psi=[yk(m-1:j-1);yk(m-2:j-2);uk(m-1:j-1);uk(m-2:j- 2)];%pn=inv(psi*psi);%theta=(inv(psi*psi)*psi*y);theta=[0 ;0;0;0];pn=10^6*eye(4);for t=3:Nps=([yk(t-1);yk(t- 2);uk(t-1);uk(t-2)]);pn=pn- pn*ps*ps*pn*(inv(1+ps*pn*ps));theta=theta+pn*ps*(yk(t)-ps*theta);thet=theta;a1=thet(1);a2=thet(2);b1=thet(3);b2= thet(4); a1t(t)=a1;a2t(t)=a2;b1t(t)=b1;b2t(t)=b2;endt=1:N;plot(t,a 1t(t),t,a2t(t),t,b1t(t),t,b2t(t));text(20,1、 47,a1);text(20,-0、67,a2);text(20,0、97,b1);text(20,0、47,b2);四．设计实验结果及分析实验结果图：仿真结果表明，大约递推到第步时，参数辨识的结果基本到稳态状态，即a1=1、5999，b1=1，c1=0、5，d1=-0、7。五、设计感受这周的课程设计告一段落了，时间短暂，意义重大。通过这次次练习的机会，重新把matlab课本看了一遍，另外学习了系统辨识的有关内容，收获颇丰。对matlab的使用更加纯熟，也锻炼了自己在课本中搜索信息和知识的能力。在设计过程中虽然遇到了一些问题，但经过一次又一次的思考，一遍又一遍的检查终于找出了原因所在，也暴露出了前期我在这方面的知识欠缺和经验不足。同时我也进一步认识了matlab软件强大的功能。在以后的学习和工作中必定有很大的用处。 移动平均法对小儿科门诊人次的预测分析 目的运用移动平均法对小儿科门诊量进行预测。方法选择我院小儿科2010～2013年各季度门诊量，采用移动平均法建立一元线性回归模型，根据回归模型对2014年各季度数据进行评估预测并做95%置信区间检测。结果2014年各季度小儿科的门诊量均在置信区间内。结论移动平均法可以为我们提供较准确的预测数据，为科室的管理提供决策依据。 标签：移动平均；门诊人次；预测 小儿科门诊是医院里门诊量最大的科室之一，同时也是对季节波动较为敏感的科室，它随着季节变化而呈现有规律性的起伏波动[1]，现根据这种规律性的波动采用移动平均法对门诊人次进行预测，①可以对医院、科室的发展规划提供参考，②可据此来调整门诊医师的做诊班次，在患者高峰来临之前做好准备，减少患者的等待时间，提高患者满意度，因而具有十分重要的意义[2]。 1 资料与方法 1.1一般资料根据我院2010～2013年每季度小儿科门诊人次数据为资料进行统计分析。数据来源我院信息科统计报表，真实可靠。 1.2统计学方法数据采用spss 19.0进行分析，以时间序列的实际值为依据，计算出移动平均的季节比率、校正系数和预测误差，求出长期趋势线，以95%的置信区间对长期趋势值进行预测。 2 结果 2.1计算移动平均比率依据门诊人次（Yt）计算4个季度的移动平均数，再对4项移动平均进行2项移动平均，以剔除时间数列中季节的变动和不规则变动的影响，计算移动平均比率，并计算平均指数的均值，本例为1.0010541；按平均1来校正移动平均比率，本例中校正系数=1/1.0010541。季节修正指数=各季节平均指数×校正系数，并求得调和时间序列（Tt），见表1、表2。 2.2计算调和时间序列并进行季节比率预测调和时间序列等于每一实际门诊数值/季节修正指数。以时间序列号为自变量，以调和时间序列值为应变量建立长期时间趋势序列直线：T=a+b×t，根據最小二乘法原则，采用spss 19.0软件计算：Tt=14984.422+524.141×t。 2.3计算误差及置信区间根据所得直线方程求出各期预测值：Yt=（14984.422+524.141×t）×季节修正指数。并再次利用最小二乘原则计算误差平方e2=（Yt-Yt）2，并对其进行合理性分析得到95%相应置信区间，得到其合理性检验并据此对14年各季度门诊人次进行预测和各季度人次分析，见表3。 ''' 最小二乘法拟合函数曲线f(x) 1、拟合多项式为：y=a0+a1*x+a2*x^2+...+ak*x^k 2、求每个点到曲线的距离之和：Loss=∑(yi-(a0+a1*x+a2*x^2+...+ak*x^k))^2 3、最优化Loss函数，即求Loss对a0,a1,...ak的偏导数为0 3.1、数学解法——求解线性方程组： 整理最优化的偏导数矩阵为：X：含有xi的系数矩阵，A：含有ai的系数矩阵，Y：含有yi的系数矩阵 求解：XA=Y中的A矩阵 3.2、迭代解法——梯度下降法： 计算每个系数矩阵A[k]的梯度，并迭代更新A[k]的梯度 A[k]=A[k]-(learn_rate*gradient) ''' import numpy as np import matplotlib.pyplot as plt plt.rcParams['font.sans-serif']=['SimHei'] plt.rcParams['axes.unicode_minus']=False ''' 高斯列主消元算法 ''' #得到增广矩阵 def get_augmented_matrix(matrix,b): row,col=np.shape(matrix) matrix=np.insert(matrix,col,values=b,axis=1) return matrix #取出增广矩阵的系数矩阵（第一列到倒数第二列） def get_matrix(a_matrix): return a_matrix[:,:a_matrix.shape[1]-1] #选列主元，在第k行后的矩阵里，找出最大值和其对应的行号和列号 def get_pos_j_max(matrix,k): max_v=np.max(matrix[k:,:]) pos=np.argwhere(matrix==max_v) i,_=pos[0] return i,max_v #矩阵的第k行后，行交换 def exchange_row(matrix,r1,r2,k): matrix[[r1,r2],k:]=matrix[[r2,r1],k:] return matrix #消元计算(初等变化) def elimination(matrix,k): row,col=np.shape(matrix) for i in range(k+1,row): m_ik=matrix[i][k]/matrix[k][k] matrix[i]=-m_ik*matrix[k]+matrix[i] return matrix #回代求解 def backToSolve(a_matrix): matrix=a_matrix[:,:a_matrix.shape[1]-1]#得到系数矩阵 b_matrix=a_matrix[:,-1]#得到值矩阵 row,col=np.shape(matrix) x=[None]*col#待求解空间X 1、最小二乘法： 1）（用1 T A A 方法计算逆矩阵） #include double mean,temp=0,variance=0; int i; for(i=0;i matlab 基本语句 1.循环语句for for i=s1:s3：s2 循环语句组 end 解释：首先给i赋值s1；然后，判断i是否介于s1与s2之间；如果是，则执行循环语句组,i=i+s3（否则，退出循环.）；执行完毕后，继续下一次循环。 例：求1到100的和，可以编程如下： sum=0 for i=1:1:100 sum=sum+i end 这个程序也可以用while语句编程。 注：for循环可以通过break语句结束整个for循环. 2.循环语句while 例：sum=0;i=1; while(i<=100) sum=sum+i;i=i+1; end 3.if语句 if(条件) 语句 end if(条件) 语句 else 语句 end if(条件) 语句 elseif 语句 end 4.关系表达式： =,>,<,>=,<=,==(精确等于) 5.逻辑表达式：|(或),&(且) 6.[n,m]=size(A)(A为矩阵) 这样可以得到矩阵A的行和列数 n=length(A)，可以得到向量A的分量个数；如果是矩阵，则得到矩阵A的行与列数这两个数字中的最大值。 7.！后面接Dos命令可以调用运行一个dos程序。 8.常见函数： poly():为求矩阵的特征多项式的函数，得到的为特征多项式的各个系数。如 a=[1,0,0;0,2,0;0,0,3]，则poly(a)＝1 -6 11 -6。相当于poly(a)＝1入^3+(-6)入^2+11入+(-6)。 compan():可以求矩阵的伴随矩阵. sin()等三角函数。 MATLAB在数学建模中的应用(3) 一、程序设计概述 MATLAB所提供的程序设计语言是一种被称为第四代编程语言的高级程序设计语言，其程序简洁，可读性很强，容易调试。同时，MATLAB的编程效率比C/C++语言要高得多。 MATLAB编程环境有很多。常用的有： 1.命令窗口 2.word窗口 3.M-文件编辑器，这是最好的编程环境。 M-文件的扩展名为“.m”。M-文件的格式分为两种： ①λ M-脚本文件，也可称为“命令文件”。 ② M-函数文件。这是matlab程序设计的主流。λ 保存后的文件可以随时调用。 二、MATLAB程序结构 按照现代程序设计的观点，任何算法功能都可以通过三种基本程序结构来实现，这三种结构是：顺序结构、选择结构和循环结构。其中顺序结构是最基本的结构，它依照语句的自然顺序逐条地执行程序的各条语句。如果要根据输入数据的实际情况进行逻辑判断，对不同的结果进行不同的处理，可以使用选择结构。如果需要反复执行某些程序段落，可以使用循环结构。 1 顺序结构 顺序结构是由两个程序模块串接构成。一个程序模块是完成一项独立功能的逻辑单元，它可以是一段程序、一个函数，或者是一条语句。 看图可知，在顺序结构中，这两个程序模块是顺序执行的，即先执行<程序 科研探索 知识创新 与109 ——科协论坛?2013年第01期（下）—— 基于移动最小二乘法的列车牵引特性曲线拟合 □ 赵笑龙 徐中伟 朱龙 (同济大学 上海 201804) 摘 要：为了精确地得到列车运行时产生的牵引力，需对列车牵引特性曲线进行拟合。详细阐述移动最小二乘法的基本原理，并以6K 型电力机车为例，分别使用分段最小二乘法和移动最小二乘法进行列车牵引特性曲线的拟合，并将拟合的结果进行对比，通过分析误差的大小及曲线的平滑性，表明该方法的优越性和有效性。关键词：移动最小二乘法牵引特性曲线曲线拟合中图分类号：O241.5 文献标识码：A 文章编号：1007-3973（2013）001-109-03 1引言 曲线拟合在工程设计、计算机图形等方面有着广泛的应用。城市轨道列车牵引力的计算一般分为图解法和分析法，分析法是直接根据物理公式进行计算，求出牵引力，而图解法是通过对牵引特性曲线的拟合，来得到牵引力。由于关于列车牵引力计算的大量公式都是通过经验得来的，因此用公式计算法并不比用图解法得到的结果精确度高。所以我们经常通过图解法来得到列车的牵引力。可以通过测量得到一组关 于速度与牵引力的数据点(),i =1,2…n ，由于不知道牵引特性曲线的原型函数f (x )，而且测量数据点也必然会有误差，因此可以根据离散数据点作出拟合曲线。曲线拟合方法中多使用传统的最小二乘法，即求解一个方程组，使得曲线的误差平方和最小。通过求解该方程组，即可得到期望的拟合曲线。但是传统最小二乘法在拟合的过程中有着较大的计算量，并且在数据量比较大的时候，由于形状比较复杂，可能需要分段拟合，这样会造成拟合曲线的不连续。为了克服上述缺点，本文将使用移动最小二乘法(MLS )来拟合列车牵引曲线。2移动最小二乘法(MLS )的逼近原理 本节主要介绍移动最小二乘法(MLS )的基本原理及其数学表达，以任一计算点处拟合系数的计算为例，来解释移动最小二乘法(MLS )的处理方式及思想。 设全局近似函数为(x )，则它在点x 的局部近似为 (1) 其中：(x )(i =1,2,…,m )是m 次完备单项式基函数，是相应的系数，这些系数是空间坐标x 的函数。基函数(x )常取单项式，在一维问题中： 线性基： ()=[1]，=2(2)二次基： ()=[12]，=3 (3) 为了体现不同节点对计算点的影响，引入了权函数，它是一个紧支撑函数，只和以i 为中心的有效区域内的结点值相关，并与节点到计算点的距离有关。其中，高斯型、紧支撑径向基函数、样条函数和指数型是常用的权函数。 在移动最小二乘法（MLS ）的中，的选取使近似函数在计算点x 领域内对待求函数f (x )的局部逼近误差最小。定义 (4) 令J 取最小值，则 (5) 由此可得 (6) 将式(4)写为矩阵形式，可以得到()=1()()( 7) 其中 (8)()=[(-1)(1)，(-2)(2)，…，(-)()(9) =[(1)，(2)， …， ()](10)()=1()，2()，…，()(11)可见，点x 处的系数有两方面决定：一是该节点值， 二是它的权函数。对不同的节点，随着权函数的变化，系数也跟着不停地变化， 相应的近似函数表示为() =( )()=()1()()(o ) (12) 将移动最小二乘法(MLS )和分段最小二乘法的拟合进行对比，可以发现：分段拟合是将节点域分成了若干段，这就等于对分段的节点分别加权；并且在每个分段内，各节点对拟合系数有着相同的影响，即加权函数取为矩形。而在移动最小二乘法(MLS )中，通过引入紧支撑权函数，对每个计算点附近的节点进行加权，计算点即为中心点，权函数的支撑域即为对中心点有影响的点的范围，随着计算点的改变，权函数的支撑域不断地移动；同时，在每个权函数内，各节点对于计算点的影响与权函数取值有关，因此权函数可以发挥一个软分段的作用。 3移动最小二乘法(MLS )的实施方案 在进行移动最小二乘法(MLS )拟合时,根据具体情况的不同，主要要考虑两方面的因素： (1)基函数的选择：基函数可以分为线性基、二次基和高次基。通常，拟合曲线的平滑性会随着基函数次数的增高越来越好，但是，它会造成就计算量的急剧增大，并且还有可能会出现拟合方程组病态的情况。因此，在拟合过程中，最多取三次基。本文将选择二次基进行实验仿真。(2)权函数的选择：在移动二乘法中，权函数起着至关重要的作用。移动二乘法中的权函数需要满足以下条件：1)在支撑域内，( - ) ＞ 0； 最小二乘法的多项式拟 合m a t l a b实现 Document serial number【NL89WT-NY98YT-NC8CB-NNUUT-NUT108】 用最小二乘法进行多项式拟合（matlab 实现） 西安交通大学 徐彬华 算法分析： 对给定数据 (i=0 ,1,2,3,..,m),一共m+1个数据点，取多项式P(x),使 函数P(x)称为拟合函数或最小二乘解，令似的 使得 其中，a0，a1,a2,…,an 为待求未知数，n 为多项式的最高次幂，由此，该问题化为求 的极值问题。由多元函数求极值的必要条件： j=0,1,…,n 得到： j=0,1,…,n 这是一个关于a0，a1,a2,…,an 的线性方程组，用矩阵表示如下： 因此，只要给出数据点 及其个数m ，再给出所要拟合的参数n ，则即可求出未知数矩阵（a0，a1,a2,…,an ） 试验题1 编制以函数 为基的多项式最小二乘拟合程序，并用于对下列数据作三次多项式最小二乘拟合(取权函数wi ≡1) x i y i 总共有7个数据点，令m=6 第一步：画出已知数据的的散点图，确定拟合参数n; x=::;y=[,,,,,,]; plot(x,y,'*') xlabel 'x 轴' ylabel 'y 轴' title '散点图' hold on {} n k k x 0= 因此将拟合参数n设为3. 第二步：计算矩阵 A= 注意到该矩阵为（n+1）*(n+1)矩阵， 多项式的幂跟行、列坐标（i,j）的关系为i+j-2,由此可建立循环来求矩阵的各个元素，程序如下： m=6;n=3; A=zeros(n+1); for j=1:n+1 for i=1:n+1 for k=1:m+1 A(j,i)=A(j,i)+x(k)^(j+i-2) end end 九九乘法表的C语言代码，黄路平编写与2012.3.6 代码一：#include case 6:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 7:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 8:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; case 9:for (i=1;i<=j;i++) printf ("%d*%d=%d\t",i,j,i*j); printf("\n"); break; } } } MATLAB机械工程 最小二乘法曲线拟合的应用实例 班级: 姓名: 学号: 指导教师: 一，实验目的 通过Matlab上机编程，掌握利用Matlab软件进行数据拟合分析及数据可视化方法 二，实验内容 1.有一组风机叶片的耐磨实验数据，如下表所示，其中X为使用时间，单位为小时h，Y为磨失质量，单位为克g。要求： 对该数据进行合理的最小二乘法数据拟合得下列数据。 x=[10000 11000 12000 13000 14000 15000 16000 17000 18000 19000 2 0000 21000 22000 23000]; y=[24.0 26.5 29.8 32.4 34.7 37.7 41.1 42.8 44.6 47.3 65.8 87.5 137.8 174. 2] 三，程序如下 X=10000:1000:23000; Y=[24.0,26.5,29.8,32.4,34.7,37.7,41.1,42.8,44.6,47.3,65.8,87.5,137.8,17 4.2] dy=1.5; %拟合数据y的步长for n=1:6 [a,S]=polyfit(x,y,n); A{n}=a; da=dy*sqrt(diag(inv(S.R′*S.R))); Da{n}=da′; freedom(n)=S.df; [ye,delta]=polyval(a,x,S); YE{n}=ye; D{n}=delta; chi2(n)=sum((y-ye).^2)/dy/dy; end Q=1-chi2cdf(chi2,freedom); %判断拟合良好度 clf,shg subplot(1,2,1),plot(1：6,abs(chi2-freedom),‘b’) xlabel(‘阶次’)，title(‘chi2与自由度’) subplot(1,2,2),plot(1：6,Q,‘r’,1：6,ones(1,6)*0.5) xlabel(‘阶次’),title(‘Q与0.5线’) nod=input(‘根据图形选择适当的阶次（请输入数值）’)； elf,shg, plot(x,y,‘kx’)；xlabel(‘x’),ylabel(‘y’)； axis([8000,23000,20.0,174.2])；hold on errorbar(x,YE{nod},D{nod},‘r’)；hold off title(‘较适当阶次的拟合’) text(10000,150.0,[‘chi2=’num2str(chi2(nod))‘~’int2str(freedom(nod))]) 为明确解释变量和随机误差各产生的效应是多少，统计上把数据点与它在回归直线上相应位置的差异称为残差，把每个残差的平方和称为残差平方和，它表示随机误差的效应。20()n i i S y y == -∑ 设所示直线议程为y ax b =+，最小二乘法就是示使得残差平方和2 1[()]n i i i M y ax b ==-+∑最 小时a 和b 的值。把M 看作a 和b 的函数，通过求多元函数偏导求最小值时的a 和b 。 1(,)2[()]0n a i i i M M a b y ax b a =?==--+=?∑ 1 (,)2[()]0n b i i i M M a b y a x b b =?==--+=?∑ 即：2111n n n i i i i i i i a x b x x y ===+=∑∑∑ 11 n n i i i i a x nb y ==+=∑∑ 化简得：11122 11()n n n i i i i i i i n n i i i i n x y x y a n x x =====-=-∑∑∑∑∑ b y a x =- 求出其取极值时的a 和b 值取可得直线拟合的方程。其C 语言程序如下： #include 曲线拟合的最小二乘法 学院：光电信息学院 姓名：赵海峰 学号： 200820501001 一、曲线拟合的最小二乘法原理： 由已知的离散数据点选择与实验点误差最小的曲线 S( x) a 0 0 ( x) a 1 1(x) ... a n n ( x) 称为曲线拟合的最小二乘法。 若记 m ( j , k ) i (x i ) j (x i ) k (x i ), 0 m (f , k ) i0 (x i )f (x i ) k (x i ) d k n 上式可改写为 ( k , jo j )a j d k ; (k 0,1,..., n) 这个方程成为法方程，可写成距阵 形式 Ga d 其中 a (a 0,a 1,...,a n )T ,d (d 0,d 1,...,d n )T , 、 数值实例： 下面给定的是乌鲁木齐最近 1个月早晨 7:00左右(新疆时间 )的天气预报所得 到的温度数据表，按照数据找出任意次曲线拟合方程和它的图像。 它的平方误差为： || 2 | 2 ] x ( f (2008 年 10 月 26~11 月 26) F 面应用Matlab 编程对上述数据进行最小二乘拟合 三、Matlab 程序代码: x=[1:1:30]; y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1]; %三次多项式拟合% %九次多项式拟合% %十五次多项式拟合% %三次多项式误差平方和 % %九次次多项式误差平方和 % %十五次多项式误差平方和 % %用*画出x,y 图像% %用红色线画出x,b1图像% %用绿色线画出x,b2图像% %用蓝色o 线画出x,b3图像% 四、数值结果： 不同次数多项式拟和误差平方和为: r1 = 67.6659 r2 = 20.1060 r3 = 3.7952 r1、r2、r3分别表示三次、九次、十五次多项式误差平方和 拟和曲线如下图: a 仁polyfit(x,y,3) a2= polyfit(x,y,9) a3= polyfit(x,y,15) b1= polyval(a1,x) b2= polyval(a2,x) b3= polyval(a3,x) r1= sum((y-b1).A 2) r2= sum((y-b2).A2) r3= sum((y-b3).A2) plot(x,y,'*') hold on plot(x,b1, 'r') hold on plot(x,b2, 'g') hold on plot(x,b3, 'b:o') Private Sub Command1_Click() For i = 1 To 4 For j = 1 To 6 Print "*"; Next j Print Next i End Sub Private Sub Command2_Click() For i = 1 To 4 Print Tab(i); For j = 1 To 6 Print "*"; Next j Print Next i End Sub Private Sub Command3_Click() For i = 1 To 4 Print Tab(5 - i); For j = 1 To 6 Print "*"; Next j Print Next i End Sub Private Sub Command4_Click() For i = 1 To 9 For j = 1 To 9 Print i; "*"; j; "="; i * j; Next j Print Next i End Sub Private Sub Command5_Click() Dim se As String Print Tab(35); "乘法表" For i = 1 To 9 For j = 1 To 9 se = i & "*" & j & "=" & i * j Print Tab((j - 1) * 9); se; Next j Picture1.Print Next i End Sub Private Sub Command6_Click() End End Sub Private Sub Command7_Click() Print Tab(35); "乘法表" For i = 1 To 9 For j = i To 9 se = i & "*" & j & "=" & i * j Print Tab((j - 1) * 9); se; Next j Print Next i End Sub Private Sub Command8_Click() Print Tab(35); "乘法表" For i = 1 To 9 For j = 1 To i se = i & "*" & j & "=" & i * j Print Tab((j - 1) * 9); se; Next j Print Next i End Sub Private Sub Command9_Click() Cls End Sub Private Sub Picture1_Click() Dim se As String Picture1.Print Tab(35); "乘法表" For i = 1 To 9 For j = 1 To 9 se = i & "*" & j & "=" & i * j Picture1.Print Tab((j - 1) * 9); se; Next j 移动最小二乘法 2.1 移动最小二乘曲线拟合 将拟合函数表述为如下形式： 1 ()()()()()m T i i i f x p x a x p x a x ===∑, (3) 其中a (x )=(a 1(x ), a 2(x ),…, a m (x ))T 为待定系数，p (x )=(p 1(x ), p 2(x ),…, p m (x ))T 为基函数向量，通常需要选择完备多项式基，例如二维情况 线性基 p (x ) = (1, x , y )T (m =3) 二次基 p (x ) = (1, x , y , x 2, xy , y 2)T (m=6) 为了得到较为精确的局部近似值，需使局部近似值f (x i )和节点值y i 之差平方带权最小，因此残差的离散加权L 2范式为： 2 21 1 ()[()]()[()()]n n T i i i i i i i J w x x f x y w x x p x a x y ===--=--∑∑, (4) 其中n 是求解区域内的节点数，f (x )是拟合函数，w (x －x i )是节点x i 的权函数。 权函数应该是非负的，且随着2 i x x -的增加单调递减，权函数还应该具有紧支性，即 在支持域（x 的影响区域）内不等于0，在支持域之外全为0，一般选用圆形作为权函数的 支持域，半径记为r 。常用的权函数是样条函数，记i s x x '=-，s s r ' =，则三次样条函数形式如下： 23 23 214432 4 41 ()4413 32 0 1. s s s s s s s s s ω?-+≤ ?? ?=-+-<≤??>? ?? (5) 要求出待定系数a (x ),先要使J 取得最小值，先将（4）式写成矩阵形式： J = (Pa (x )－Y )T W (x ) (Pa (x )－Y ) 其中 Y = (y 1, y 2,…, y n )T , W (x ) = diag (w 1(x ), w 2(x ),…, w n (x ))，w i (x ) =()i w x x -. 1121112 22212()()()()()()() ()()m m n n m n p x p x p x p x p x p x P p x p x p x ??????=?? ? ??? . 根据最小二乘原理求得待定系数为： a (x ) = A -1(x )B (x )Y 其中A (x ) = P T W (x ) P , B (x ) = P T W (x )。最小二乘法曲线拟合 原理及matlab实现
最小二乘法及其应用..
php九九乘法表
"; for($z = 0;$z <9 - $j;$z++){ echo " "; } echo "";*/ //第三种 /*echo " "; } for($i = $j;$i >= 1;$i--){ echo "$i*$j=".$i*$j." "; } echo ""; for($j = 9;$j >= 1;$j--){ echo "
";*/ /*第二种 for($j = 9;$j >= 1;$j--){ for($i = 1; $i <= $j; $i++){ //echo "$i x $j = ".$i * $j; //echo " "; echo "$i x $j =".$i*$j." "; } echo ""; for($z = 1;$z <= 9 - $j;$z++ ){ echo " "; } echo " "; //echo "$i*$j=".$i*$j."";$i*$j=".$i*$j." "; } echo "
"; } */ /*第一种 for($j = 1;$j <= 9;$j++){用matlab实现最小二乘递推算法辨识系统参数
移动平均法对小儿科门诊人次的预测分析
最小二乘算法程序Python版
最小二乘法的编程实现
matlab循环语句
基于移动最小二乘法的列车牵引特性曲线拟合
最小二乘法的多项式拟合matlab实现
九九乘法表的C语言代码
Matlab最小二乘法曲线拟合的应用实例
最小二乘法C语言程序
曲线拟合的最小二乘法matlab举例
九九乘法表源代码(vb)
移动最小二乘法
- 最小二乘算法程序Python版
- 最小二乘法程序说明及流程图
- 最小二乘法matlab程序
- 移动最小二乘法的程序
- 最小二乘法程序
- 基于移动最小二乘法的曲线曲面拟合_曾清红
- 最小二乘法算法程序
- 各种最小二乘法汇总(算例及MATLAB程序)
- 各种最小二乘法汇总(算例及MATLAB程序)00177
- 移动最小二乘法
- ESPRIT算法(最小二乘法)matlab程序
- 最小二乘法递推最小二乘法C语言源程序
- 移动最小二乘法MLS程序
- 移动最小二乘法
- 修正的移动最小二乘方法
- 基于移动最小二乘法的列车牵引特性曲线拟合
- 递推最小二乘法与广义最小二乘法matlab源程序
- 最小二乘法程序编写.doc
- 最小二乘法C语言的实现
- 最小二乘法及matlab程序