怎么用C语言实现最小二乘法?

  最小二乘法常用于根据实测数据求线性方程的最近似解。根据如图（图片引用于百度百科）的描述，利用C语言求，使用最小二乘法算法求线性方程的解，程序如下：

#include <stdio.h>

#define N 4  //共有4个记录，根据需要增加记录

typedef struct Data{ //定义实验记录结构 

int w //实验次数 

double x  

double y

}DATA

//根据d中的n个DATA记录，计算出线性方程的a,b两值

void getcs(DATA *d,int n,double &a,double &b){

double fi11=0,fi12=0,fi21=0,fi22=0,f1=0,f2=0

int i

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

fi11+=d[i].w

fi12+=d[i].w*d[i].x

fi21=fi12

fi22+=d[i].w*d[i].x*d[i].x

f1+=d[i].w*d[i].y

f2+=d[i].w*d[i].x*d[i].y

}

//解一元一次方程

b=(f2*fi11/fi21-f1)/(fi22*fi11/fi21-fi12)

a=(f2*fi12/fi22-f1)/(fi21*fi12/fi22-fi11)

}

int main(){

DATA d[N]={  //定义时赋初值，共4个记录

{2,0.1,1.1},

{1,0.2,1.9},

{1,0.3,3.1},

{1,0.4,3.9}

}

double a,b

getcs(d,N,a,b)  //计算线性方程参数a,b

printf("线性方程是：Y=%.4lf+%.4lfX\n",a,b)

}

#include <stdio.h>

#include <conio.h>

#include <math.h>

#include <process.h>

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

{

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

for(j=0j<T+1j++)

{

res=0

for(k=0k<Nk++)

{

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=0i<Ni++)

{

for(j=0j<T+1j++)

{

b[j][i]=a[i][j]

}

}

}

void init(float x_y[][2],int n)

{

int i

printf("请输入%d个已知点：\n",N)

for(i=0i<ni++)

{

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=0i<Ni++)

{

for(j=0j<T+1j++)

{

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

{

for(j=0j<T+1j++)

{

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

{

for(j=ij<nj++)

{

if(argu[i-1][i-1]==0)

{

for(p=ip<np++)

{

if(argu[p][i-1]!=0)

break

}

if(p==n)

{

printf("方程组无解!\n")

exit(0)

}

for(t=0t<n+1t++)

{

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-1k<n+1k++)

{

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

{

temp=argu[i][n]

for(j=n-1j>ij--)

{

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

{

temp=0

for(j=0j<Nj++)

{

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

{

for(j=0j<T+2j++)

{

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

{

printf("a[%d]=%g,",i+1,a[i])

}

printf("\n")

printf("拟合曲线方程为:\ny(x)=%g",a[0])

for(i=1i<ni++)

{

printf(" + %g",a[i])

for(j=0j<ij++)

{

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()

}

]]>

</Content>

<PostDateTime>2007-4-19 19:23:57</PostDateTime>

</Reply>

<Reply>

<PostUserNickName></PostUserNickName>

<rank>一级(初级)</rank>

<ranknum>user1</ranknum>

<credit>100</credit>

<ReplyID>40389872</ReplyID>

<TopicID>5478010</TopicID>

<PostUserId>1526752</PostUserId>

<PostUserName>jiangxc2004</PostUserName>

<Point>0</Point>

<Content>

<![CDATA[

你可以改一下

不从终端输入,直接在程序中给出参数

请输入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

]]>

</Content>

<PostDateTime>2007-4-19 19:26:11</PostDateTime>

</Reply>

<Reply>

<PostUserNickName></PostUserNickName>

<rank>一级(初级)</rank>

<ranknum>user1</ranknum>

<credit>100</credit>

<ReplyID>40390406</ReplyID>

<TopicID>5478010</TopicID>

<PostUserId>1526752</PostUserId>

<PostUserName>jiangxc2004</PostUserName>

<Point>0</Point>

<Content>

<![CDATA[

这样就可以直接调用process()函数了!

二次拟合的话就把宏 T 成2

拟合点的数目 N 也可以修改!

也可以去到注释的部分进行返回值的调用!

#include "stdafx.h"

#include <stdio.h>

#include <conio.h>

#include <stdlib.h>

#include <cmath>

#include <iostream>

using namespace std

void polyfit(int n, double x[], double y[], int poly_n, double a[])

{

int i, j

double *tempx, *tempy, *sumxx, *sumxy, *ata

void gauss_solve(int n, double A[], double x[], double b[])

tempx = new double[n]

sumxx = new double[poly_n * 2 + 1]

tempy = new double[n]

sumxy = new double[poly_n + 1]

ata = new double[(poly_n + 1)*(poly_n + 1)]

for (i = 0i<ni++)

{

tempx[i] = 1

tempy[i] = y[i]

}

for (i = 0i<2 * poly_n + 1i++)

for (sumxx[i] = 0, j = 0j<nj++)

{

sumxx[i] += tempx[j]

tempx[j] *= x[j]

}

for (i = 0i<poly_n + 1i++)

for (sumxy[i] = 0, j = 0j<nj++)

{

sumxy[i] += tempy[j]

tempy[j] *= x[j]

}

for (i = 0i<poly_n + 1i++)

for (j = 0j<poly_n + 1j++)

ata[i*(poly_n + 1) + j] = sumxx[i + j]

gauss_solve(poly_n + 1, ata, a, sumxy)

delete [] tempx

tempx = NULL

delete [] sumxx

sumxx = NULL

delete [] tempy

tempy = NULL

delete [] sumxy

sumxy = NULL

delete [] ata

ata = NULL

}

void gauss_solve(int n, double A[], double x[], double b[])

{

int i, j, k, r

double max

for (k = 0k<n - 1k++)

{

max = fabs(A[k*n + k])/*find maxmum*/

r = k

for (i = k + 1i<n - 1i++)

if (max<fabs(A[i*n + i]))

{

max = fabs(A[i*n + i])

r = i

}

if (r != k)

for (i = 0i<ni++)         /*change array:A[k]&A[r] */

{

max = A[k*n + i]

A[k*n + i] = A[r*n + i]

A[r*n + i] = max

}

max = b[k]                   /*change array:b[k]&b[r]     */

b[k] = b[r]

b[r] = max

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

{

for (j = k + 1j<nj++)

A[i*n + j] -= A[i*n + k] * A[k*n + j] / A[k*n + k]

b[i] -= A[i*n + k] * b[k] / A[k*n + k]

}

}

for (i = n - 1i >= 0x[i] /= A[i*n + i], i--)

for (j = i + 1, x[i] = b[i]j<nj++)

x[i] -= A[i*n + j] * x[j]

}

/*==================polyfit(n,x,y,poly_n,a)===================*/

/*=======拟合y=a0+a1*x+a2*x^2+……+apoly_n*x^poly_n========*/

/*=====n是数据个数 x y是数据值 poly_n是多项式的项数======*/

/*===返回a0,a1,a2,……a[poly_n]，系数比项数多一（常数项）=====*/

void main()

{

int  n = 9, poly_n = 2

//double x[20] = { 1 ,2   ,  3  ,   4   ,  7 ,8 ,9,11,12, 13,15,15,17 ,17, 19  ,18 ,20,19,20, 20 },

//y[20] = { 1,3,6 ,10 ,17,25,34,45,57,70,85,100,117,134,153,171,191,210 ,  230 ,  250 }

double x[9]={1,3,4,5,6,7,8,9,10},

y[9]={10,5,4,2,1,1,2,3,4}

double a[50]

polyfit(n, x, y, poly_n, a)

for (int i = 0i <poly_n + 1i++)/*这里是升序排列，Matlab是降序排列*/

cout <<"a[" <<i <<"]=" <<a[i] <<endl

}

运行结果，我是拟合的2次的，你可以拟合多次。

方程式：

0.267571*x*x-3.60531*x+13.4597=0

这个2次多项式的最低点还用我说吗，在那个区间上，自己代入；

---