用C語(yǔ)言實(shí)現(xiàn)的曲線擬合的最小二乘法

1、精選優(yōu)質(zhì)文檔-傾情為你奉上實(shí)驗(yàn)名稱(chēng)：曲線擬合的最小二乘法實(shí)驗(yàn)?zāi)康牧私馇€擬合的最小二乘法實(shí)驗(yàn)類(lèi)型設(shè)計(jì)型實(shí)驗(yàn)環(huán)境Windows XP TC實(shí)驗(yàn)內(nèi)容相關(guān)知識(shí)：已知Ca,b中函數(shù)f(x)的一組實(shí)驗(yàn)數(shù)據(jù)（xi,yi）（i=0,1,m），其中yi=f(xi)。設(shè)是Ca,b上線性無(wú)關(guān)函數(shù)族。在中找函數(shù)f(x) 曲線擬合的最小二乘解，其法方程（組）為： 其中， k=0,1,n特別是，求函數(shù)f(x) 曲線擬合的線性最小二乘解的計(jì)算公式為：數(shù)據(jù)結(jié)構(gòu)：兩個(gè)一維數(shù)組或一個(gè)二維數(shù)組算法設(shè)計(jì)：（略）實(shí)驗(yàn)用例： 已知函數(shù)y=f(x)的一張表：x0102030405060708090y6867.166.465.664.66

2、1.861.060.860.460試驗(yàn)要求：利用曲線擬合的線性最小二乘法求被逼近函數(shù)f(x)在點(diǎn)x=55處的近似值，并畫(huà)出實(shí)驗(yàn)數(shù)據(jù)和直線。編寫(xiě)代碼：#include<stdio.h>#include<stdlib.h> #include<graphics.h>double qiuhe1(double a102,int p)int i;double y;y=0;for(i=0;i<10;i+)y=y+aip;return y;double qiuhe2(double a102,int p)int i;double y=0;for(i=0;i<10;

3、i+)y=y+ai0*aip;return y;double nihe(double a102,double x)double a1,b,y;a1=(10*qiuhe2(a,1)-qiuhe1(a,0)*qiuhe1(a,1)/(10*qiuhe2(a,0)-qiuhe1(a,0)*qiuhe1(a,0);b=(qiuhe2(a,0)*qiuhe1(a,1)-qiuhe1(a,0)*qiuhe2(a,1)/(10*qiuhe2(a,0)-qiuhe1(a,0)*qiuhe1(a,0);y=a1*x+b;return y;int main()double a102=0,68,10,67.1,20

4、,66.4,30,65.6,40,64.6,50,61.8, 60,61.0,70,60.8,80,60.4,90,60; double x,x1,q=1; char c12; int i; long n; int arw6=515,235,520,240,515,245; int arw16=315,45,320,40,325,45; int gdriver=IBM8514; int gmode=IBM8514HI; initgraph(&gdriver, &gmode, "c:TC20BGI");cleardevice();printf("in

5、put x:n");scanf("%lf",&x);printf("%fn",nihe(a,x);n=nihe(a,x)*+1;c0='y'c1='=' c4='.' for(i=10;i>1;i-) if(i!=4) ci=n%10+48; n=n/10; c11='0'x1=x;setbkcolor(7); setcolor(14);setlinestyle(0,0,3);drawpoly(3,arw);drawpoly(3,arw1);line(120,240,

6、520,240);line(320,40,320,440);x=0;setcolor(2);setlinestyle(0,0,1); line(0+320),(int)(240-nihe(a,0)*q),(90+320),(int)(240-nihe(a,90)*q); setcolor(3); outtextxy(320,30,"Y"); outtextxy(310,245,"O"); outtextxy(525,240,"X"); outtextxy(x1+330),(240-nihe(a,x1)-10),c); settexts

7、tyle(4,0,4); outtextxy(450,400,"Nihe.");for(i=0;i<=9;i+) putpixel(ai0+320),(240-ai1*q),11);setcolor(4);setlinestyle(1,0,1);line(x1+320),(240-nihe(a,x1)*q)-80),(x1+320),(240-nihe(a,x1)*q)+120);getch();closegraph();實(shí)驗(yàn)結(jié)果（測(cè)試用例、實(shí)驗(yàn)結(jié)果）實(shí)驗(yàn)總結(jié)與心得通過(guò)本次實(shí)驗(yàn)，對(duì)曲線擬合的最小二乘法有了更深刻的了解！概念最小二乘法多項(xiàng)式曲線擬合，根據(jù)給定的m個(gè)點(diǎn),并

8、不要求這條曲線精確地經(jīng)過(guò)這些點(diǎn)，而是曲線y=f(x)的近似曲線y= (x)。原理原理部分由個(gè)人根據(jù)互聯(lián)網(wǎng)上的資料進(jìn)行總結(jié)，希望對(duì)大家能有用     給定數(shù)據(jù)點(diǎn)pi(xi,yi)，其中i=1,2,m。求近似曲線y= (x)。并且使得近似曲線與y=f(x)的偏差最小。近似曲線在點(diǎn)pi處的偏差i= (xi)-y，i=1,2,.,m。 常見(jiàn)的曲線擬合方法:     1.使偏差絕對(duì)值之和最小          2.使偏差絕對(duì)值最大的最小     

9、     3.使偏差平方和最小          按偏差平方和最小的原則選取擬合曲線，并且采取二項(xiàng)式方程為擬合曲線的方法,稱(chēng)為最小二乘法。推導(dǎo)過(guò)程：     1. 設(shè)擬合多項(xiàng)式為：               2. 各點(diǎn)到這條曲線的距離之和，即偏差平方和如下：              &

10、#160;3. 為了求得符合條件的a值，對(duì)等式右邊求ai偏導(dǎo)數(shù)，因而我們得到了：                                              .               4. 

11、將等式左邊進(jìn)行一下化簡(jiǎn)，然后應(yīng)該可以得到下面的等式：                                         .               5. 把這些等式表示成矩陣的形式，就可以得到下面的矩陣：   

12、           6. 將這個(gè)范德蒙得矩陣化簡(jiǎn)后可得到:               7. 也就是說(shuō)X*A=Y，那么A = (X'*X)-1*X'*Y，便得到了系數(shù)矩陣A，同時(shí)，我們也就得到了擬合曲線。實(shí)現(xiàn)運(yùn)行前提:1. Python運(yùn)行環(huán)境與編輯環(huán)境；2. Matplotlib.pyplot圖形庫(kù)，可用于快速繪制2D圖表，與matlab中的plot命令類(lèi)似，而且用法也基本相同。代碼:python

13、  1. # coding=utf-8  2.   3. ''''' 4. 作者:Jairus Chan 5. 程序:多項(xiàng)式曲線擬合算法 6. '''  7. import matplotlib.pyplot as plt  8. import math  9. import numpy  10

14、. import random  11.   12. fig = plt.figure()  13. ax = fig.add_subplot(111)  14.   15. #階數(shù)為9階  16. order=9  17.   18. #生成曲線上的各個(gè)點(diǎn)  19. x = numpy.arange(-1,1,0.02)  2

15、0. y = (a*a-1)*(a*a-1)*(a*a-1)+0.5)*numpy.sin(a*2) for a in x  21. #ax.plot(x,y,color='r',linestyle='-',marker='')  22. #,label="(a*a-1)*(a*a-1)*(a*a-1)+0.5"  23.   24. #生成的曲線上的各個(gè)點(diǎn)偏移一下，并放入到xa,ya中去

16、  25. i=0  26. xa=  27. ya=  28. for xx in x:  29.     yy=yi  30.     d=float(random.randint(60,140)/100  31.     #ax.plot(xx*d,yy*d,color='m',line

17、style='',marker='.')  32.     i+=1  33.     xa.append(xx*d)  34.     ya.append(yy*d)  35.   36. '''''for i in range(0,5): 37.  &#

18、160;  xx=float(random.randint(-100,100)/100 38.     yy=float(random.randint(-60,60)/100 39.     xa.append(xx) 40.     ya.append(yy)'''  41.   42. ax.plot(xa,ya,color='m',lines

19、tyle='',marker='.')  43.   44.   45. #進(jìn)行曲線擬合  46. matA=  47. for i in range(0,order+1):  48.     matA1=  49.     for j in range(0,order+1): 

20、; 50.         tx=0.0  51.         for k in range(0,len(xa):  52.             dx=1.0  53.    

21、60;        for l in range(0,j+i):  54.                 dx=dx*xak  55.             tx+

22、=dx  56.         matA1.append(tx)  57.     matA.append(matA1)  58.   59. #print(len(xa)  60. #print(matA00)  61. matA=numpy.array(matA)  62.   63. matB= &#

23、160;64. for i in range(0,order+1):  65.     ty=0.0  66.     for k in range(0,len(xa):  67.         dy=1.0  68.       

24、60; for l in range(0,i):  69.             dy=dy*xak  70.         ty+=yak*dy  71.     matB.append(ty)  72.    73. matB=numpy.array(matB)  74.   75. matAA=numpy.linalg.solve(matA,matB)  76.   77. #畫(huà)出擬合后的曲線  78. #print(matAA)  79. xxa= numpy.arange(-1,1.06,0.01)  80. yya=