求C或C++语言编写的用最小二乘法进行曲线拟合

求C或C++语言编写的用最小二乘法进行曲线拟合,第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会有个常数项的。

简单的讲,所谓拟合是指已知某函数的若干离散函数值{f1,f2,…,fn},通过调整该函数中若干待定系数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=1i<=ni++)

{

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=1i<=mi++)

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=1i<=ni++)

z=z+x[i-1]/n/*z为各个x的平均值*/

b[0]=1

d1=n

p=0

c=0

for(i=1i<=ni++)

{

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=1i<=ni++)

{

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=3j<=mj++)

{

s[j-1]=t[j-2]

s[j-2]=-p*t[j-2]+t[j-3]

if(j>=4)

for(k=j-2k>=2k--)

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=1i<=ni++)

{

q=s[j-1]

for(k=j-1k>=1k--)

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-1k>=1k--)

{

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=1i<=ni++)

{

q=a[m-1]

for(k=m-1k>=1k--)

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)

}

一元线性回归的C语言程序是:利用最小二乘法来估计线性回归方程的参数,然后用这些参数来预测因变量的值1。例如,你可以参考下面的代码:

#include <stdio.h>#include <math.h>//定义一个函数,计算一元线性回归方程的参数a和bvoid linear_regression(double x[], double y[], int n, double *a, double *b){//定义变量

double sum_x = 0//x的和

double sum_y = 0//y的和

double sum_xy = 0//xy的和

double sum_x2 = 0//x平方的和

//遍历数组,计算各项和

for (int i = 0i <ni++)

{

sum_x += x[i]

sum_y += y[i]

sum_xy += x[i] * y[i]

sum_x2 += x[i] * x[i]

}//根据最小二乘法公式,计算a和b

*a = (n * sum_xy - sum_x * sum_y) / (n * sum_x2 - pow(sum_x, 2))

*b = (sum_y - (*a) * sum_x) / n

}//主函数int main(){//定义一个自变量数组x,存放观测值

double x[] = {1.0, 2.0, 3.0, 4.0} //定义一个因变量数组y,存放观测值

double y[] = {3.1, 4.9, 7.2, 8.9} //定义数组长度n

int n = sizeof(x) / sizeof(x[0]) //定义两个指针变量a和b,用来存放线性回归方程的参数

double a double b


欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/yw/12012180.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-05-20
下一篇 2023-05-20

发表评论

登录后才能评论

评论列表(0条)

保存