『壹』 求一個最小二乘法c語言程序
#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=1;i<n;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<T+1;i++)
for(j=0;j<T+1;j++)
{
res=0;
for(k=0;k<N;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<N;i++)
{
for(j=0;j<T+1;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<n;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<N;i++)
{
for(j=0;j<T+1;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<n;i++)
{
for(j=0;j<T+1;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<n;i++)
{
for(j=i;j<n;j++)
{
if(argu[i-1][i-1]==0)
{
for(p=i;p<n;p++)
{
if(argu[p][i-1]!=0)
break;
}
if(p==n)
{
printf("方程組無解!\n");
exit(0);
}
for(t=0;t<n+1;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<n+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<n;i++)
{
temp=0;
for(j=0;j<N;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<T+1;i++)
{
for(j=0;j<T+2;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<n;i++)
{
printf("a[%d]=%g,",i+1,a[i]);
}
printf("\n");
printf("擬合曲線方程為:\ny(x)=%g",a[0]);
for(i=1;i<n;i++)
{
printf(" + %g",a[i]);
for(j=0;j<i;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();
}
]]>
</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 也可以修改!
也可以去到注釋的部分進行返回值的調用!
『貳』 系統辨識中遞推最小二乘法的C語言編程問題
'/0'是多字元字元常量，有的編譯器不支持，即便支持出來的效果也前插萬別。不過就這上下文的話，因為是列印%c所以只考慮低位元組'/0'和'0'就沒區別了（多數編譯器低位元組在後面，和字元串相反）。'0'是0x30，也就是48，68+48時116也就是't'
以上作為你沒寫錯考慮的
如果是'\0'的話那就是0，68+0還是0。68對應的是'D'
我可以幫助你，你先設置我最佳答案後，我網路Hii教你。
『叄』 計算方法中最小二乘法如何用C語言編程
#include <stdio.h>
#include <math.h>
#define epsilon 1e-6
void nihe1(int n,int m,float sum_x,float sum_y,float sum_xy,float x2);
void nihe2(int n,int m,float sum_x,float sum_y,float sum_xy,float x2,float x2y,float x3,float x4);
int main(){
float x[100]={0.0};
float y[100]={0.0};
int n,i,flag=1;
float sum_y=0.0,sum_x=0,x2=0,sum_xy=0.0,x3=0,x4=0,x2y=0.0;
printf("請你輸入需要測試的數據(先輸入x[]，後輸入y[])的個數:");
scanf("%d",&n);
for(i = 0; i < n; i++){
scanf("%f",&x[i]);}
for(i = 0; i < n; i++){
scanf("%f",&y[i]);}
for(i = 0; i < n; i++){
sum_x += x[i];
sum_y += y[i];
sum_xy += x[i]*y[i];
x2 += x[i]*x[i];
x2y += x[i]*x[i]*y[i];
x3 += x[i]*x[i]*x[i];
x4 += x[i]*x[i]*x[i]*x[i];}
printf("---------------請你輸入的要擬合的函數------------------\n");
printf(" 1、擬合一次函數\n");
printf(" 2、擬合二次函數\n");
scanf("%d",&flag);
switch(flag){
case 1:
nihe1(n,flag+1,sum_x,sum_y,sum_xy,x2); break;
case 2:
nihe2(n,flag+1,sum_x,sum_y,sum_xy,x2,x2y,x3,x4); break;
default:
printf("ERROR\n");}
return 0;}
void nihe1(int n,int m,float sum_x,float sum_y,float sum_xy,float x2){
int i,k,j;
float t,s=0;
float a[2][3] = {{(float)n,sum_x,sum_y},{sum_x,x2,sum_xy}};
n=m;
//if(m == 3)
// a[3][4] = {{n,sum_x,sum_y},{sum_x,x2,x3,sum_xy},{x2,x3,x4,x2y}};
for(k=0;k<n-1;k++) {
for(i=k+1;i<n;i++)
if( abs((int)a[i][k]) > abs((int)a[k][k]) )
for(j=k;j<n+1;j++) {
t=a[k][j];
a[k][j]=a[i][j];
a[i][j]=t; }
if( abs((int)a[k][k]) < epsilon) {
printf("\nError,主元消去法 cann't be rable,break at %d!\n",k+1);
return; }
for(i=k+1;i<n;i++){
a[i][k]=a[i][k] / a[k][k];
for(j=k+1;j<n+1;j++)
a[i][j]=a[i][j]-a[i][k] * a[k][j]; }}
a[n-1][n]=a[n-1][n] / a[n-1][n-1];
for(k=n-2;k>=0;k--) {
s=0;
for(j=k+1;j<n;j++)
s+=a[k][j]*a[j][n];
a[k][n]=( a[k][n]-s ) / a[k][k]; }
printf("\n*****The Result*****\n");
for(i=0;i<n;i++)
printf(" x[%d]=%.4f\n",i+1,a[i][n]);
printf("函數為：p(x) = %.4f + (%.4f)*x\n",a[0][n],a[1][n]);
getchar();}
void nihe2(int n,int m,float sum_x,float sum_y,float sum_xy,float x2,float x2y,float x3,float x4){
int i,k,j;
float t,s=0;
float a[3][4]=
{{(float)n,sum_x,x2,sum_y},{sum_x,x2,x3,sum_xy},{x2,x3,x4,x2y}};
n=m;
for(k=0;k<n-1;k++) {
for(i=k+1;i<n;i++)
if( abs((int)a[i][k]) > abs((int)a[k][k]) )
for(j=k;j<n+1;j++) {
t=a[k][j];
a[k][j]=a[i][j];
a[i][j]=t; }
if( abs((int)a[k][k]) < epsilon) {
printf("\nError,主元消去法 cann't be rable,break at %d!\n",k+1);
return; }
for(i=k+1;i<n;i++){
a[i][k]=a[i][k] / a[k][k];
for(j=k+1;j<n+1;j++)
a[i][j]=a[i][j]-a[i][k] * a[k][j]; } }
a[n-1][n]=a[n-1][n] / a[n-1][n-1];
for(k=n-2;k>=0;k--) {
s=0;
for(j=k+1;j<n;j++)
s+=a[k][j]*a[j][n];
a[k][n]=( a[k][n]-s ) / a[k][k]; }
printf("\n*****The Result*****\n");
for(i=0;i<n;i++)
printf(" x[%d]=%.4f\n",i+1,a[i][n]);
printf("函數為：p(x) = %.4f + (%.4f)*x + (%.4f)*x*x\n",a[0][n],a[1][n],a[2][n]);
getchar();}
『肆』 怎麼用C語言實現最小二乘法
最小二乘法常用於根據實測數據求線性方程的最近似解。根據如圖（圖片引用於網路）的描述，利用C語言求，使用最小二乘法演算法求線性方程的解，程序如下：
#include<stdio.h>
#defineN4//共有4個記錄，根據需要增加記錄
typedefstructData{//定義實驗記錄結構
intw;//實驗次數
doublex;
doubley;
}DATA;
//根據d中的n個DATA記錄，計算出線性方程的a,b兩值
voidgetcs(DATA*d,intn,double&a,double&b){
doublefi11=0,fi12=0,fi21=0,fi22=0,f1=0,f2=0;
inti;
for(i=0;i<n;i++){
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);
}
intmain(){
DATAd[N]={//定義時賦初值，共4個記錄
{2,0.1,1.1},
{1,0.2,1.9},
{1,0.3,3.1},
{1,0.4,3.9}
};
doublea,b;
getcs(d,N,a,b);//計算線性方程參數a,b
printf("線性方程是：Y=%.4lf+%.4lfX
",a,b);
}
『伍』 急求用maple編寫最小二乘法的程序。
其實最小二乘法也可以寫得不像我原來的那麼麻煩：
fid=fopen('input.txt','r+');
U=fscanf(fid,'%f',a);
fclose(fid);
fid1=fopen('output1.txt','r+');
Y=fscanf(fid1,'%f',a);
fclose(fid1);
U1=U(1:a-1);
U2=U(1:a-2);
Y1=Y(1:a-1);
Y2=Y(1:a-2);
D1=-1*[0;Y1];
D2=-1*[0;0;Y2];
D3=[0;U1];
D4=[0;0;U2];
D=[D1 D2 D3 D4];
Q=inv(D'*D)*D'*Y
下面是遞推最小二乘的演算法：
m=input('m=')
fid=fopen('input.txt','r+');
U=fscanf(fid,'%f',a);
fclose(fid);
fid1=fopen('output1.txt','r+');
Y=fscanf(fid1,'%f',a);
fclose(fid1);
U1=U(1:m-1);
U2=U(1:m-2);
Y1=Y(1:m-1);
Y2=Y(1:m-2);
Y3=Y(1:m);
D1=-1*[0;Y1];
D2=-1*[0;0;Y2];
D3=[0;U1];
D4=[0;0;U2];
D=[D1 D2 D3 D4];
Q=inv(D'*D)*D'*Y3;
P=inv(D'*D);
%以上程序是用最小二乘法計算的初值，取前m個數據%
for i=m:a-1;
x=[-1*Y(i);-1*Y(i-1);U(i);U(i-1)];
y=Y(i+1);
p=1/(1+x'*P*x);
Z=Q+P*x*p*(y-x'*Q);
P=P-P*x*p*x'*P;
q=norm(Z-Q)/norm(Q);
Q=Z;
if q<10e-6;
k=i-m
q
Q
break
end
end
ppp='help me!'
『陸』 用c語言編寫最小二乘法
#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 = 0; i<n; i++)
{
tempx[i] = 1;
tempy[i] = y[i];
}
for (i = 0; i<2 * poly_n + 1; i++)
for (sumxx[i] = 0, j = 0; j<n; j++)
{
sumxx[i] += tempx[j];
tempx[j] *= x[j];
}
for (i = 0; i<poly_n + 1; i++)
for (sumxy[i] = 0, j = 0; j<n; j++)
{
sumxy[i] += tempy[j];
tempy[j] *= x[j];
}
for (i = 0; i<poly_n + 1; i++)
for (j = 0; j<poly_n + 1; j++)
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 = 0; k<n - 1; k++)
{
max = fabs(A[k*n + k]); /*find maxmum*/
r = k;
for (i = k + 1; i<n - 1; i++)
if (max<fabs(A[i*n + i]))
{
max = fabs(A[i*n + i]);
r = i;
}
if (r != k)
for (i = 0; i<n; i++) /*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 + 1; i<n; i++)
{
for (j = k + 1; j<n; j++)
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 - 1; i >= 0; x[i] /= A[i*n + i], i--)
for (j = i + 1, x[i] = b[i]; j<n; j++)
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 = 0; i < poly_n + 1; i++)/*這里是升序排列，Matlab是降序排列*/
cout << "a[" << i << "]=" << a[i] << endl;
}
運行結果，我是擬合的2次的，你可以擬合多次。
方程式：
0.267571*x*x-3.60531*x+13.4597=0
這個2次多項式的最低點還用我說嗎，在那個區間上，自己代入；