ここから本文です

数値データ(x,y)があり、それを関数y=a+b*x+c*x*x+d*x*x*x+e*x*x*x*xの多項式関数...

yuc********さん

2012/6/618:43:59

数値データ(x,y)があり、それを関数y=a+b*x+c*x*x+d*x*x*x+e*x*x*x*xの多項式関数でフィッティングし、パラメータa,b,c,d,eを抽出したいと思っています。プログラム言語は現在Fortranを使用中です。

gnuplotではできるのですが、プログラムでパラメータ抽出を行いたいと思っています。参考になるようなソースやサイト等があれば教えていただきたいと思います。
アドバイスよろしくお願いします。

閲覧数:
232
回答数:
2
お礼:
25枚

違反報告

ベストアンサーに選ばれた回答

tra********さん

2012/6/810:36:58

私の使っているのは最小自乗法です。
マトリクス演算と連立一次方程式を解くことで係数が求まります。

昔作ったのでFortran77ですが参考までにソースを添付します。
c -----------------------------
subroutine aproxl(n,m,x,y,cc)
c -----------------------------
parameter(mxnw=61)
dimension x(1),y(1),cc(1)
dimension amat(3000,mxnw),a(mxnw,mxnw),c(mxnw)
c
if( n.gt.3000 ) stop 'Maximum Number of Samples is 3000'
if( m.gt. 60 ) stop 'Maximum Jisuu is 19'
c
m1=m+1
do i=1,n
do k=0,m
amat(i,k+1)=x(i)**k
end do
end do
call multtn(amat,amat,a,n,m1,m1,3000,3000,mxnw)
call multtn(amat,y ,c,n,m1, 1,3000,3000,mxnw)
call psimeq(m1,a,c,mxnw)
do i=1,m1
cc(i)=c(i)
end do
c
return
end

c ---------------------------------------
subroutine multtn(a,b,c,l,m,n,l1,m1,n1)
c ---------------------------------------
dimension a(l1,1), b(m1,1), c(n1,1)
c
do i=1,m
do j=1,n
c(i,j)=0.0
do k=1,l
c(i,j)=c(i,j)+a(k,i)*b(k,j)
end do
end do
end do
c
return
end

c ---------------------------
subroutine psimeq(n,a,c,n1)
c ---------------------------
parameter(mxnw=150)
dimension a(n1,n1), c(n1)
dimension b(mxnw), km(mxnw), kl(mxnw), ker(mxnw)
c
m=0
mer=0
nn=n+1
c
do 100 k=1,n
do 101 i=k,n
if( a(i,k).eq.0.0 ) goto 101
aa=a(i,k)
do j=k,n
a(i,j)= a(i,j)/aa
end do
c(i)= c(i)/aa
101 continue
c
if( a(k,k).ne.0.0 ) go to 290
kk=k
c ll=0
260 kk=kk+1
if( kk.eq.nn ) go to 285
if(a(kk,k).eq.0.0) go to 260
ll=kk
do j=1,n
b(j)= a(ll,j)
end do
b(nn)= c(ll)
c
do j=1,n
a(ll,j)=a(k,j)
a(k ,j)=b(j)
end do
c(ll)=c(k)
c(k) =b(nn)
m=m+1
km(m)=k
kl(m)=ll
go to 290
285 mer=mer+1
ker(mer)=k
go to 100
290 if( k.eq.n ) go to 100
kk=k+1
do 102 i=kk,n
if( a(i,k).eq.0.0 ) goto 102
do j=k,n
a(i,j)=a(i,j)-a(k,j)
end do
c(i)= c(i)-c(k)
102 continue
100 continue
c
if( mer.eq.0 ) then
do 103 k=1,n
kk=n+1-k
if( kk.eq.n ) goto 103
l=kk+1
do j=l,n
c(kk)= c(kk)-a(kk,j)*c(j)
end do
103 continue
return
else
write(*,600) mer
do i=1,mer
write(*,3000) ker(i)
end do
do i=1,m
write(*,2000) km(i),kl(i)
end do
stop
end if
c
600 format(/6x,'* Simeq Error *'
& /11x,'Number of Zero Pivot =',i3)
2000 format(11x,'Change Row',i3,' to Row',i3)
3000 format(11x,'Zero Pivot Row =',i3)
c
end

この回答は投票によってベストアンサーに選ばれました!

ベストアンサー以外の回答

1〜1件/1件中

spo********さん

2012/6/702:00:29

多項式の計算には「掃き出し法」と言う方法が古くから使われています。
私もスプライン曲線による補間などで使ったことがありますが、行列式を使った方法です。

ググればすぐ出てきますが
http://www.asahi-net.or.jp/~uc3k-ymd/Lesson/Section03/SweepOut.html

のサイトなどはわかり易かったです。

みんなで作る知恵袋 悩みや疑問、なんでも気軽にきいちゃおう!

Q&Aをキーワードで検索:

Yahoo! JAPANは、回答に記載された内容の信ぴょう性、正確性を保証しておりません。
お客様自身の責任と判断で、ご利用ください。
本文はここまでです このページの先頭へ

「追加する」ボタンを押してください。

閉じる

※知恵コレクションに追加された質問は選択されたID/ニックネームのMy知恵袋で確認できます。

不適切な投稿でないことを報告しました。

閉じる