Transcript 数値計算法I - 九州工業大学

数値計算法
第7回
最小２乗法近似
2006年度
九州工業大学工学部電気電子工学コース用講義資料
講師：趙孟佑
[email protected]
1
関数近似
• n組の与えられたデータ(xi,yi)(i=1,2,--,n)を最もよ
く近似する式（関数）を探し出す。
3
2
y
1
0
0
2
4
6
8
10
12
x
4次の多項式による近似例
2
3
4
y  1.68  0.54 x  0.18 x  0.024 x  0.0012 x
2
多項式近似
• n組のデータをm次の多項式で近似する。
– 普通はn≥m
f ( x )  am x
m
 L a 2 x 2  a1 x  a 0
• この近似式にxiを代入して得られる値とデータyiの差
– 偏差
yi  f ( x i )
3
多項式近似
• n個全てのデータで、偏差が0となるには、次式が成
り立つ必要
a m x1  L a 2 x1  a1 x1  a 0

y1
a m x 2  L a 2 x 2  a1 x 2  a 0

y2
M

M

yn 1

yn
m
2
m
式の数
n個
2
1
1
a m x n  1  L a 2 x n  1  a1 x n  1  a 0
m
am x
2
m
n
1
 L a 2 x  a1 x  a 0
2
n
1
n
未知数の数m個
n≠m
普通は、n個全ての点でぴったりと左辺と右辺が
一致することはない。
4
多項式近似
• 偏差の2乗和(residual sum of square)S
2
n
S 
 y
i
 f ( xi )
i 1
Sの値が最小になる近似式f(x)、即ち、最適な係
数a0,a1,a2,----,amの組み合わせを見つける
5
多項式近似
2
n
S 
 y
i
 f ( xi )
i 1
• Sの値が最小になる近似式f(x)を探す
– 最適な係数a0,a1,a2,----,amを見つける
a m x1  L a 2 x1  a1 x1  a 0

y1
 L a 2 x  a1 x  a 0

y2

M
a m x n  1  L a 2 x n  1  a1 x n  1  a 0

yn 1
a m x n  L a 2 x n  a1 x n  a 0

yn
m
式の数
n個
am x
m
2
2
2
2
1
1
2
M
m
2
m
1
2
1
未知数の数m個
n≠m
6
多項式近似
• 偏差の2乗和(residual sum of square)S
2
n
S 
 y
i
 f ( xi )
i 1

m
1
 L a 2 x  a 1 x  a 0  y1


m
2
 L a 2 x  a1 x  a 0  y 2


m
n 1

m
n
 am x
 am x
2
1
1
1
2
2
1
2
2
2
M
 am x
 am x
 L a2 x
2
n 1
 a1 x
1
n 1
 a0  yn 1
 L a 2 x  a1 x  a 0  y n
2
n
1
n

2

2
x1,x2,----,xn-1,xn,y1,y2,----yn-1,ynは与えられた値なで、
Sは係数a1,a2,---,am-1,amのm個の変数からなる関数と
見做せる
7
多項式近似
• Sが最小とは、
dS 
S
a1
da1 
S
a 2
da 2  L 
S
a m  1
da m  1 
S
a m
da m
でdSがゼロとなることである。
つまり、Sのa1,a2---の偏微分がゼロになればよい。
8
S
a 0


2 a 0  2 a m x 1  L a 2 x 1  a 1 x 1  y1
m

2
1

 2 a 0  2 a m x 2  L a 2 x 2  a1 x 2  y 2
m
2
1

M

 2 a
 2 a 0  2 a m x n  1  L a 2 x n  1  a1 x n  1  y n  1
2 a0
m
2
1
x  L a 2 x n  a1 x n  y n
m n
m
n
2
1
n


n
 2 n a 0  2 a1  x i  2 a 2  x i  L  2 a m  1  x i
m 1
2
i 1
i 1
i 1
n
n
 2 a m  x i  2  yi
m
i 1
i1
9
S
a 1


2 a 1 x 1  2 x 1 a m x 1  L a 2 x 1  a 0  y1
2
1
m

2

 2 a1 x 2  2 x 2 a m x 2  L a 2 x 2  a 0  y 2
2
1
m
2

M

 2 a1 x n  1  2 x n  1 a m x n  1  L a 2 x n  1  a 0  y n  1
2
1

m
2
 2 a1 x n  2 x n a m x n  L a 2 x n  a 0  y n
2
1
n
m
2
n
n


n
n
 2 a 0  x i  2 a1  x i  2 a 2  x i  L  2 a m  1  x i  2 a m  x i
2
i 1
i 1
3
i 1
m 1
m
i 1
i1
n
 2  x i yi
i 1
10
S
a 2


2 a 2 x 1  2 x 1 a m x 1  L a 1 x 1  a 0  y1
4
2
m

1

 2 a 2 x 2  2 x 2 a m x 2  L a1 x 2  a 0  y 2
4
2
m
1

M

 2 a 2 x n  1  2 x n  1 a m x n  1  L a1 x n  1  a 0  y n  1
4
2
m

1
 2 a 2 x n  2 x n a m x n  L a1 x n  a 0  y n
4
2
n
m
1
n
n


n
 2 a 0  x i  2 a1  x i  2 a 2  x i  L  2 a m  1  x i
2
i 1
3
i 1
m 1
4
i 1
i 1
n
 2 am  xi
m2
i 1
n
 2  x i yi
2
i 1
11
S
a m

2m
2 a m x1
a
 2 x a
 2 x1
2 am x2
2m
m
m 1
x
m 1 1
m
2
 L a 1 x 1  a 0  y1
1
m 1
x
m 1 2

 L a1 x 2  a 0  y 2
1

M

m 1
 2 a m x n  1  2 x n  1 a m  1 x n  1  L a1 x n  1  a 0  y n  1
2m
2 am xn
2m
m

m 1
 2 xn a m 1xn
m
n
n
 2 a 0  x i  2 a1  x i
m 1
m
i 1
i 1
1
 L a1 x n  a 0  y n
1
n
 2 a2  xi
m2
i 1


n
 L  2 am 1  xi
2 m 1
i 1
n
 2 am  xi
2m
i 1
n
 2  x i yi
m
i1
12
多項式近似
• 結局、dS=0となる係数aの組み合わせは、以下の
式を解いて得られる。
XA  Y

n


 n
  xi
 i 1
M
X  
 n
m 1

xi

 i 1
 n

m
xi


 i 1
n
x
n
i
L
i 1
x
x
n
2
i
L
x
i 1
O
n
x
M
n
m
i
L
i 1
x
2m2
i
i 1
n
i 1
m
i
i 1
M
x
m 1
i
i 1
n
n
m 1
i
L
x
i 1

x 
i 1

n

m 1
 xi 
i1

M 

n
2 m 1
 x i 
i 1

n
2m 
x
 i 

i 1
n
2 m 1
i
m
i
 a0 


a1


A   M 


 am 1 
 a m 


y
  i 

 i1

 n
  x i yi 

 i1

M
Y  

 n
m 1

x i yi 


 i1

 n


m
x i yi




 i1
n
13
多項式近似
XA  Y
行列Xは(m+1)x(m+1)行列である。この連立一次方程式を
解いて、Aを求める。
行列の解法には、Gauss-Seidel法や反復法等を使う
14
多項式近似（例：直線近似）
3
2
y
1
0
0
2
4
6
8
10
12
x
f ( x )  a 0  a1 x で近似
15
多項式近似（例：直線近似）

 n

 n
  xi
 i 1



x
y
 i  a    i 
i 1
 0   i 1

 
n
n
  a1  

2
 xi 
  x i yi 

 i 1

i 1
n
n
I
x
y
1
0.843
1.301
2
1.896
1.353
3
2.018
0.962
4
3.068
1.293
5
4.276
0.919
6
5.536
1.433
7
6.290
1.117
8
7.712
1.251
9
8.499
1.764
10
9.559
1.362
11
10.766
2.458
n=11,m=1
 11

 60.46
60.46   a 0   15.21 
   

445.2   a1   92.57 
16
直線近似
 11

 60.46
60.46   a 0   15.21 
   

445.2   a1   92.57 
 a 0   0.9468 
  

a
0.0793
 1 

よって、
f ( x )  0.9468  0.0793 x
で近似する。
2
11
S 
 y
i 1
i
 a 0  a1 x i   1.0906
17
直線近似
3
回帰直線
2
y
1
0
0
2
4
6
8
10
12
x
f ( x )  0.9468  0.0793 x
回帰係数
18
多項式近似
I
x
y
1
0.843
1.301
2
1.896
1.353
3
2.018
0.962
4
3.068
1.293
5
4.276
0.919
6
5.536
1.433
7
6.290
1.117
8
7.712
1.251
9
8.499
1.764
10
9.559
1.362
11
10.766
2.458
n=11,m=4
f ( x )  a 0  a1 x  a 2 x  a 3 x  a 4 x
2
3
4
で近似
19
多項式近似


y
  i 
 i1

n


  x i yi 
 i 1

 n 2 
Y    x i yi 
 i1

 n

3
  x i yi 
 i 1

 n

4
  x i yi 
 i  1

n

 n

 n
  xi
 i 1
 n 2
X    xi
 i 1
 n
3
  xi
 i 1
 n
  x i4
 i  1
n

n
xi

n
2
xi

n
3
xi

i 1
i 1
i 1
i 1
n
n
n
n
x
2
i
x
3
i
x
4
i

i 1
i 1
i 1
i 1
n
n
n
n
x
3
i
i 1
4
i
i 1
n
x
x
x
5
i
i 1
n
4
i
x
i 1
n
5
i
x

n
6
i

i 1
i 1
i 1
i 1
n
n
n
n
x
i 1
5
i
x
i 1
6
i
x
i 1
7
i

i 1
4 
xi 


5
xi 


6
xi 


7
xi 


8
xi 

 a0 
 
a1
 
A   a2 
 
a
 3
 a 4 
20
多項式近似
 a 0   1 6 .7 7 2 
  

a1
 0 .5 4 2 2 8
  

A   a 2    0 .1 8 3 3 7 
  

 a 3    0 .0 2 4 6 8 
 a 4   0 .0 0 1 1 9 
2
11
S 

y i  a 0  a1 x i  a1 x i  a1 x i  a1 x i
2
3
4

 0.4711
i 1
21
多項式近似
3
2
y
1
0
0
2
4
6
8
10
12
x
22
指数近似
n組の与えられたデータ(xi,yi)(i=1,2,--,n)を最もよく近似する
指数関数
f ( x )  ae
bx
を探しだす。
偏差の2乗和Sは
n
S 
 y
i
 ae
bx i

2
i 1
23
指数近似
6
4
y
2
0
0
2
4
6
8
10
12
x
24
指数近似
データの自然対数をとって考える。
( x i , yi )
( x i , log( y i ))
10
結局は
直線近似
y
1
0
2
4
6
x
8
10
12
25
指数近似
Y i  log( y i ) として、
( x i , yi )
f ( x )  ae
( x i ,Y i )
bx
の対数をとって、
F ( x )  log a  bx  A  bx
結局、 ( x i ,Y i ) というデータの組を A  b x
という直線で近似するのと同じ。
26
指数近似
n
S 
 log y
 ( A  bxi )
2
i
i 1
を最小にするAとbを見つける。
S
A
 0,
S
b
 0
を満たせばよい。
27
指数近似
n
S 
 log y
 ( A  bxi )
2
i
i 1

 n

 n
  xi
 i 1



 x i   A    lo g y i 
i 1
i 1




 
n
n
b



 
2
 xi 
  x i lo g y i 

 i 1

i 1
n
n
28
指数近似

2
x
 i
 i 1
 n
xi    xi
 i 1
n
A
1
n
n
n
 
b
2
  n
x
 i   xi 
i 1
i 1
i 1

2
x
  i  log y i   x i  x i log y i
i 1
i 1
i 1
 i 1
n
n
n

2
n  xi   xi  xi

i 1
i 1
i 1
 
n
n
n

 n  x i log y i   x i  log y i
i 1
i 1
i 1

n
n
n
2

n  xi   xi  xi


i 1
i 1
i 1
n
n
n
n


  x i    log y i 
i 1
  i 1

 n

n    x i log y i 
  i 1

n












n
ae
A
よりaを求める
29
指数近似(例）
I
1
x
0.35849
y
1.7937
2
1.358
2.0022
3
2.9812
2.4203
4
3.8039
2.2281
5
4.3609
2.6195
6
5.0844
2.5143
7
6.155
2.9978
8
7.3467
3.4408
9
8.9731
4.0702
10
11
9.702
10.102
4.014
4.4716
30
指数近似
6
4
y
2
0
0
2
4
6
8
10
12
x
f ( x )  1.72 e
0.091 x
31
直線近似のプログラム例
implicit real*8 (a-h,o-z)
parameter(nmax=100,mmax=10)
dimension xdata(nmax),ydata(nmax)
dimension acof(mmax)
dimension xx(mmax,mmax),yy(mmax)
c
c* input of x and y data
c
n=11
open(unit=10,file='input.dat',status='old')
do i=1,n
read(10,*)idum,xdata(i),ydata(i)
end do
close(10)
c
c* polynominal dimension
c
m=1
c
c* definition of dS/da
c
do i=1,m+1
yy(i)=0
do ii=1,n
yy(i)=yy(i)+ydata(ii)*xdata(ii)**(i-1)
end do
c
c
do j=1,m+1
xx(i,j)=0
do ii=1,n
xx(i,j)=xx(i,j)+xdata(ii)**(i-2+j)
end do
end do
end do
c
call gauss(xx,yy,acof,m+1)
c
resid=0
do i=1,n
fvalue=0
do j=1,m+1
fvalue=fvalue+acof(j)*xdata(i)**(j-1)
end do
resid=resid+(ydata(i)-fvalue)*(ydata(i)-fvalue)
end do
write(6,*)(acof(j),j=1,m+1)
write(6,*)' residual=',resid
stop
end
32
直線近似のプログラム例
subroutine gauss(a,b,x,n)
implicit real*8 (a-h,o-z)
parameter(mmax=10)
dimension a(mmax,mmax),x(mmax),b(mmax)
dimension atmp(mmax)
c
c
c* scaling
c
do i=1,n
colmax=abs(a(i,1))
do j=2,n
if(abs(a(i,j)).gt.colmax)colmax=abs(a(i,j))
end do
do j=1,n
a(i,j)=a(i,j)/colmax
end do
b(i)=b(i)/colmax
end do
c
c* forward deleting
c
do j=1,n-1
c
c* pivotting
aamax=abs(a(j,j))
imax=j
do i=j+1,n
if(dabs(a(i,j)).gt.aamax)then
aamax=dabs(a(i,j))
imax=i
end if
end do
c
c* exchange
c
do jj=1,n
atmp(jj)=a(j,jj)
end do
btmp=b(j)
c
do jj=1,n
a(j,jj)=a(imax,jj)
end do
b(j)=b(imax)
c
do jj=1,n
a(imax,jj)=atmp(jj)
end do
b(imax)=btmp
c
do i=j+1,n
do jj=j+1,n
a(i,jj)=a(i,jj)-a(j,jj)*a(i,j)/a(j,j)
end do
b(i)=b(i)-b(j)*a(i,j)/a(j,j)
end do
end do
c
c* back substitution
c
x(n)=b(n)/a(n,n)
c
do i=n-1,1,-1
sum=0
do ii=i+1,n
sum=sum+a(i,ii)*x(ii)
end do
x(i)=(b(i)-sum)/a(i,i)
end do
c
return
end
33
直線近似のプログラム例
input.datの中身
1
2
3
4
5
6
7
8
9
10
11
0.84323 1.3011
1.8955 1.3532
2.0183 0.96193
3.0679 1.2927
4.2760 0.91855
5.5359 1.4333
6.2897 1.1165
7.7115 1.2508
8.4991 1.7638
9.5590 1.3617
10.766 2.4583
34