Keyword: 最小二乗, 双３次, スプライン, フィット
概要
本サンプルは最小二乗双３次スプライン曲面フィットを行うFortranによるサンプルプログラムです。 本サンプルは以下に示されるデータについて最小二乗双３次スプライン曲面フィットを行い、スプラインの値を求めて出力します。
※本サンプルはnAG Fortranライブラリに含まれるルーチン e02daf() のExampleコードです。本サンプル及びルーチンの詳細情報は e02daf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はe02daf のマニュアルページを参照）| このデータをダウンロード |
E02DAF Example Program Data
0.000001
30
8
10
-0.52 0.60 0.93 10.
-0.61 -0.95 -1.79 10.
0.93 0.87 0.36 10.
0.09 0.84 0.52 10.
0.88 0.17 0.49 10.
-0.70 -0.87 -1.76 10.
1.00 1.00 0.33 1.
1.00 0.10 0.48 1.
0.30 0.24 0.65 1.
-0.77 -0.77 -1.82 1.
-0.23 0.32 0.92 1.
-1.00 1.00 1.00 1.
-0.26 -0.63 8.88 1.
-0.83 -0.66 -2.01 1.
0.22 0.93 0.47 1.
0.89 0.15 0.49 1.
-0.80 0.99 0.84 1.
-0.88 -0.54 -2.42 1.
0.68 0.44 0.47 1.
-0.14 -0.72 7.15 1.
0.67 0.63 0.44 1.
-0.90 -0.40 -3.34 1.
-0.84 0.20 2.78 1.
0.84 0.43 0.44 1.
0.15 0.28 0.70 1.
-0.91 -0.24 -6.52 1.
-0.35 0.86 0.66 1.
-0.16 -0.41 2.32 1.
-0.35 -0.05 1.66 1.
-1.00 -1.00 -1.00 1.
-0.5
0.0
- １行目はタイトル行で読み飛ばされます。
- ２行目に一次方程式の有効ランクを決定するための閾値(eps)を指定しています。
- ３行目にデータ点の数(m)を指定しています。
- ４行目に変数xのノット数の合計(px)を指定しています。
- ５行目に変数yのノット数の合計(py)を指定しています。
- ６〜３５行目にデータ点yの座標、xの座標、fの座標、重み(w)を指定しています。
- ３６〜３７行目に変数yに関する内部ノット(mu)を指定しています。
出力結果
（本ルーチンの詳細はe02daf のマニュアルページを参照）| この出力例をダウンロード |
E02DAF Example Program Results
Interior Y-knots
-0.5000
0.0000
Interior X-knots
None
Sum of squares of residual RHS 1.47E+01
Rank 22
X and Y have been interchanged
X Y Data Fit Residual
-0.9500 -0.6100 -1.7900 -1.7931 -0.31E-02
-0.8700 -0.7000 -1.7600 -1.7521 0.79E-02
-0.7700 -0.7700 -1.8200 -2.4301 -0.61E+00
-0.6300 -0.2600 8.8800 7.6346 -0.12E+01
-0.6600 -0.8300 -2.0100 -1.5815 0.43E+00
-0.5400 -0.8800 -2.4200 -2.6795 -0.26E+00
-0.7200 -0.1400 7.1500 7.5708 0.42E+00
-1.0000 -1.0000 -1.0000 -1.0228 -0.23E-01
-0.4000 -0.9000 -3.3400 -4.6955 -0.14E+01
-0.2400 -0.9100 -6.5200 -4.7072 0.18E+01
-0.4100 -0.1600 2.3200 2.7039 0.38E+00
-0.0500 -0.3500 1.6600 2.2865 0.63E+00
0.6000 -0.5200 0.9300 0.9441 0.14E-01
0.8700 0.9300 0.3600 0.3529 -0.71E-02
0.8400 0.0900 0.5200 0.5024 -0.18E-01
0.1700 0.8800 0.4900 0.4705 -0.20E-01
1.0000 1.0000 0.3300 0.6315 0.30E+00
0.1000 1.0000 0.4800 1.4910 0.10E+01
0.2400 0.3000 0.6500 0.9241 0.27E+00
0.3200 -0.2300 0.9200 -0.3692 -0.13E+01
1.0000 -1.0000 1.0000 1.0835 0.84E-01
0.9300 0.2200 0.4700 1.4912 0.10E+01
0.1500 0.8900 0.4900 0.4414 -0.49E-01
0.9900 -0.8000 0.8400 0.5495 -0.29E+00
0.4400 0.6800 0.4700 1.5862 0.11E+01
0.6300 0.6700 0.4400 0.6288 0.19E+00
0.2000 -0.8400 2.7800 1.7123 -0.11E+01
0.4300 0.8400 0.4400 0.6888 0.25E+00
0.2800 0.1500 0.7000 0.7713 0.71E-01
0.8600 -0.3500 0.6600 0.9347 0.27E+00
Sum of squared residuals 1.47E+01
Spline coefficients
-1.0228 115.4668 -433.5558 -68.1973
24.8426 -140.1485 258.5042 15.6756
-29.4878 132.2933 -173.5103 20.0983
9.9575 -51.6200 67.6666 -5.8765
10.0577 4.7543 -15.3533 -0.3260
1.0835 -2.7932 7.7708 0.6315
- ３〜５行目にyの内部ノットが出力されています。
- ７〜８行目にxの内部ノットがないことが示されています。
- １０行目に残差平方和が出力されています。
- １２行目にランクが出力されています。
- １４行目にxとyが入れ替えられたことが示されています。
- １６〜４６行目にx、y、f、フィットされた値、残差が出力されています。
- ４８行目に残差平方和が出力されています。
- ５０〜５６行目にスプライン係数が出力されています。
ソースコード
（本ルーチンの詳細はe02daf のマニュアルページを参照）
※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
PROGRAM e02dafe
! E02DAF Example Program Text
! Mark 23 Release. nAG Copyright 2011.
! .. Use Statements ..
USE nag_library, ONLY : e02daf, e02def, e02zaf, nag_wp
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Parameters ..
INTEGER, PARAMETER :: nin = 5, nout = 6
CHARACTER (1), PARAMETER :: label(2) = (/ 'X', 'Y'/)
! .. Local Scalars ..
REAL (KIND=nag_wp) :: eps, sigma, sum, temp
INTEGER :: i, iadres, ifail, itemp, j, m, &
nadres, nc, npoint, nws, px, py, rank
! .. Local Arrays ..
REAL (KIND=nag_wp), ALLOCATABLE :: c(:), dl(:), f(:), ff(:), lamda(:), &
mu(:), w(:), wrk(:), ws(:), x(:), y(:)
INTEGER, ALLOCATABLE :: adres(:), iwrk(:), point(:)
! .. Executable Statements ..
WRITE (nout,*) 'E02DAF Example Program Results'
! Skip heading in data file
READ (nin,*)
! Read data, interchanging X and Y axes if PX < PY
READ (nin,*) eps
READ (nin,*) m
READ (nin,*) px, py
IF (px<py) THEN
itemp = px
px = py
py = itemp
itemp = 1
ELSE
itemp = 0
END IF
nadres = (px-7)*(py-7)
npoint = m + (px-7)*(py-7)
nc = (px-4)*(py-4)
nws = (2*nc+1)*(3*py-6) - 2
ALLOCATE (lamda(px),mu(py),x(m),y(m),f(m),ff(m),w(m),dl(nc),c(nc), &
ws(nws),point(npoint),adres(nadres),wrk(py-4),iwrk(py-4))
IF (itemp==1) THEN
READ (nin,*) (y(i),x(i),f(i),w(i),i=1,m)
IF (py>8) THEN
READ (nin,*) mu(5:(py-4))
END IF
IF (px>8) THEN
READ (nin,*) lamda(5:(px-4))
END IF
ELSE
READ (nin,*) (x(i),y(i),f(i),w(i),i=1,m)
IF (px>8) THEN
READ (nin,*) lamda(5:(px-4))
END IF
IF (py>8) THEN
READ (nin,*) mu(5:(py-4))
END IF
END IF
! Sort points into panel order
ifail = 0
CALL e02zaf(px,py,lamda,mu,m,x,y,point,npoint,adres,nadres,ifail)
WRITE (nout,*)
WRITE (nout,99995) 'Interior ', label(itemp+1), '-knots'
IF (px==8) THEN
WRITE (nout,*) 'None'
ELSE
DO j = 5, px - 4
WRITE (nout,99996) lamda(j)
END DO
END IF
WRITE (nout,*)
WRITE (nout,99995) 'Interior ', label(2-itemp), '-knots'
IF (py==8) THEN
WRITE (nout,*) 'None'
ELSE
DO j = 5, py - 4
WRITE (nout,99996) mu(j)
END DO
END IF
! Fit bicubic spline to data points
ifail = 0
CALL e02daf(m,px,py,x,y,f,w,lamda,mu,point,npoint,dl,c,nc,ws,nws,eps, &
sigma,rank,ifail)
WRITE (nout,*)
WRITE (nout,99999) 'Sum of squares of residual RHS', sigma
WRITE (nout,*)
WRITE (nout,99998) 'Rank', rank
! Evaluate spline at the data points
ifail = 0
CALL e02def(m,px,py,x,y,lamda,mu,c,ff,wrk,iwrk,ifail)
sum = 0.0E0_nag_wp
IF (itemp==1) THEN
WRITE (nout,*)
WRITE (nout,*) 'X and Y have been interchanged'
END IF
! Output data points, fitted values and residuals
WRITE (nout,*)
WRITE (nout,*) &
' X Y Data Fit Residual'
DO i = 1, nadres
iadres = i + m
LOOP: DO
iadres = point(iadres)
IF (iadres<=0) THEN
EXIT LOOP
END IF
temp = ff(iadres) - f(iadres)
WRITE (nout,99997) x(iadres), y(iadres), f(iadres), ff(iadres), &
temp
sum = sum + (temp*w(iadres))**2
END DO LOOP
END DO
WRITE (nout,*)
WRITE (nout,99999) 'Sum of squared residuals', sum
WRITE (nout,*)
WRITE (nout,*) 'Spline coefficients'
DO i = 1, px - 4
WRITE (nout,99996) (c((i-1)*(py-4)+j),j=1,py-4)
END DO
99999 FORMAT (1X,A,1P,E16.2)
99998 FORMAT (1X,A,I5)
99997 FORMAT (1X,4F11.4,E11.2)
99996 FORMAT (1X,6F11.4)
99995 FORMAT (1X,A,A1,A)
END PROGRAM e02dafe