Question

matlab用最小二乘法求一形如Y=ct^m（c和m为待定系数）的多项式，使之与下列数据相拟合，

数据如下（19，19.0）（25,32.3）（31,49.0）（38,73.3）（44,97.8）

Answer 1

你的近似解析表达式为y=at+bt^2+ct^2
是不是想写成为y=at+bt^2+ct^3
但是实际拟合出来的表达式为y=a[3]+a[2]t+a[1]t^2+a[0]t^3会有个常数项的。
简单的讲，所谓拟合是指已知某函数的若干离散函数值，通过调整该函数中若干待定系数f(λ1, λ2,…,λ3), 使得该函数与已知点集的差别(最小二乘意义)最小。如果待定函数是线性，就叫线性拟合或者线性回归(主要在统计中)，否则叫作非线性拟合或者非线性回归。表达式也可以是分段函数，这种情况下叫作样条拟合。
曲线拟合：
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>

Smooth(double *x,double *y,double *a,int n,int m,double *dt1,double *dt2,double *dt3);
void main()
{
int i ,n ,m ;
double *x,*y,*a,dt1,dt2,dt3,b;
n = 12;// 12个样点
m = 4; //3次多项式拟合
b = 0; //x的初值为0
/*分别为x,y,a分配存贮空间*/
x = (double *)calloc(n,sizeof(double));
if(x == NULL)
{
printf("内存分配失败\n");
exit (0);
}
y = (double *)calloc(n,sizeof(double));
if(y == NULL)
{
printf("内存分配失败\n");
exit (0);
}
a = (double *)calloc(n,sizeof(double));
if(a == NULL)
{
printf("内存分配失败\n");
exit (0);
}
for(i=1;i<=n;i++)
{
x[i-1]=b+(i-1)*5;
/*每隔5取一个点，这样连续取12个点*/
}
y[0]=0;
y[1]=1.27;
y[2]=2.16;
y[3]=2.86;
y[4]=3.44;
y[5]=3.87;
y[6]=4.15;
y[7]=4.37;
y[8]=4.51;
y[9]=4.58;
y[10]=4.02;
y[11]=4.64;
/*x[i-1]点对应的y值是拟合已知值*/

Smooth(x,y,a,n,m,&dt1,&dt2,&dt3); /*调用拟合函数*/
for(i=1;i<=m;i++)
printf("a[%d] = %.10f\n",(i-1),a[i-1]);
printf("拟合多项式与数据点偏差的平方和为：\n");
printf("%.10e\n",dt1);
printf("拟合多项式与数据点偏差的绝对值之和为：\n");
printf("%.10e\n",dt2);
printf("拟合多项式与数据点偏差的绝对值最大值为：\n");
printf("%.10e\n",dt3);
free(x); /*释放存储空间*/
free(y); /*释放存储空间*/
free(a); /*释放存储空间*/
}

Smooth(double *x,double *y,double *a,int n,int m,double *dt1,double *dt2,double *dt3)//(x,y,a,n,m,dt1,dt2,dt3 )
//double *x; /*实型一维数组，输入参数，存放节点的xi值*/
//double *y; /*实型一维数组，输入参数，存放节点的yi值*/
//double *a; /*双精度实型一维数组，长度为m。返回m一1次拟合多项式的m个系数*/
//int n; /*整型变量，输入参数，给定数据点的个数*/
//int m; /*整型变量，输入参数，拟合多项式的项数*/
//double *dt1; /*实型变量，输出参数，拟合多项式与数据点偏差的平方和*/
//double *dt2; /*实型变量，输出参数，拟合多项式与数据点偏差的绝对值之和*/
//double *dt3; /*实型变量，输出参数，拟合多项式与数据点偏差的绝对值最大值*/
{
int i ,j ,k ;
double *s,*t,*b,z,d1,p,c,d2,g,q,dt;
/*分别为s ,t ,b分配存贮空间*/
s = (double *)calloc(n,sizeof(double));
if(s == NULL)
{
printf("内存分配失败\n");
exit (0);
}
t = (double *)calloc(n,sizeof(double));
if(t == NULL)
{
printf("内存分配失败\n");
exit (0);
}
b = (double *)calloc(n,sizeof(double));
if(b == NULL)
{
printf("内存分配失败\n");
exit (0);
}
z = 0;
for(i=1;i<=n;i++)
z=z+x[i-1]/n; /*z为各个x的平均值*/
b[0]=1;
d1=n;
p=0;
c=0;
for(i=1;i<=n;i++)
{
p=p+x[i-1]-z;
c=c+y[i-1];
}
c=c/d1;
p=p/d1;
a[0]=c*b[0];
if(m>1)
{
t[1]=1;
t[0]=-p;
d2=0;
c=0;
g=0;
for(i=1;i<=n;i++)
{
q=x[i-1]-z-p;
d2=d2+q*q;
c=y[i-1]*q+c;
g=(x[i-1]-z)*q*q+g;
}
c=c/d2;
p=g/d2;
q=d2/d1;
d1=d2;
a[1]=c*t[1];
a[0]=c*t[0]+a[0];
}
for(j=3;j<=m;j++)
{
s[j-1]=t[j-2];
s[j-2]=-p*t[j-2]+t[j-3];
if(j>=4)
for(k=j-2;k>=2;k--)
s[k-1]=-p*t[k-1]+t[k-2]-q*b[k-1];
s[0]=-p*t[0]-q*b[0];
d2=0;
c=0;
g=0;
for(i=1;i<=n;i++)
{
q=s[j-1];
for(k=j-1;k>=1;k--)
q=q*(x[i-1]-z)+s[k-1];
d2=d2+q*q;
c=y[i-1]*q+c;
g=(x[i-1]-z)*q*q+g;
}
c=c/d2;
p=g/d2;
q=d2/d1;
d1=d2;
a[j-1]=c*s[j-1];
t[j-1]=s[j-1];
for(k=j-1;k>=1;k--)
{
a[k-1]=c*s[k-1]+a[k-1];
b[k-1]=t[k-1];
t[k-1]=s[k-1];
}
}
*dt1=0;
*dt2=0;
*dt3=0;
for(i=1;i<=n;i++)
{
q=a[m-1];
for(k=m-1;k>=1;k--)
q=q*(x[i-1]-z)+a[k-1];
dt=q-y[i-1];
if(fabs(dt)>*dt3)
*dt3=fabs(dt);
*dt1=*dt1+dt*dt;
*dt2=*dt2+fabs(dt);
}
/*释放存储空间*/
free(s);
free(t);
free(b);
return(1);
}

Answer 2

程序：
clear;
clc;
t = [19 25 31 38 44];
y = [19 32.3 49 73.8 97.8];
y = y';
t = t';
st_ = [0.5 0.5]; %%%% 起始点
ft_ = fittype('c*t^m' ,... %%%% 公式
'dependent',{'y'},'independent',{'t'},... %%% 变量和自变量
'coefficients',{'c', 'm'}); %%% 系数
[cf,good]= fit(t,y,ft_ ,'Startpoint',st_)
h_ = plot(cf,'fit',0.95);
legend off; % turn off legend from plot method call
set(h_(1),'Color',[1 0 0],...
'LineStyle','-', 'LineWidth',2,...
'Marker','none', 'MarkerSize',6);
hold on,plot(t,y,'*')

结果;
cf =

General model:
cf(t) = c*t^m
Coefficients (with 95% confidence bounds):
c = 0.0588 (0.05147, 0.06612)
m = 1.96 (1.926, 1.994)

good =

sse: 0.2267
rsquare: 0.9999
dfe: 3
adjrsquare: 0.9999
rmse: 0.2749
贴图：
这个就算了，避免被拉去审，效果很好，可以自己运行一下看看