Keyword: 一般化線形モデル, 予測
概要
本サンプルはフィッティングされた一般化線形モデルを用いた予測を行うFortranによるサンプルプログラムです。 本サンプルは以下に示されるデータについて予測を行います。
※本サンプルはnAG Fortranライブラリに含まれるルーチン g02gpf() のExampleコードです。本サンプル及びルーチンの詳細情報は g02gpf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はg02gpf のマニュアルページを参照）| このデータをダウンロード |
G02GPF Example Program Data Training Data 'R' 'M' 'N' 'U' 5 1 0.0 :: LINK,MEAN,OFFSET,WEIGHT,N,M,S 1.0 25.0 2.0 10.0 3.0 6.0 4.0 4.0 5.0 3.0 :: End of X,Y 1 :: ISX 0 1.0E-6 5.0E-5 0 :: IPRINT,EPS,TOL,MAXIT Prediction Data 2 .TRUE. 'N' 'U' :: N,VFOBS,OFFSET,WEIGHT 32.0 18.0 :: End of X
- １行目はタイトル行で読み飛ばされます。
- ２〜１０行目はトレーニングデータを指定しています。
- ３行目はどのリンク関数が使用されるか（link='R':相互リンク（reciprocal link））、一般化線形モデルの引数に切片が含まれるかどうか（mean='M':切片を含める）、オフセットが必要かどうか（offset='N':不要）、重みづけをするかどうか（weight='U':重みづけをしない）、観測値の数（n=5）、独立変数の数（m=1）、モデルのスケール引数（s=0.0:スケール引数が残差平均平方から推定される）を指定しています。
- ４〜８行目に独立変数の観測値(x)と従属変数の観測値(y)を指定しています。
- ９行目はモデルに独立変数が含まれるかどうか(isx)を指定しています。"1"は含まれることを意味します。
- １０行目は反復の情報の出力が必要かどうか（iprint=0:不要）、独立変数がフルランクかどうかそうでない場合ランクは何か（eps=1.0E-6）、モデルのフィットに必要な正確さ（tol=5.0E-5）、最大反復数（maxit=0:デフォルト値10を使用）を指定しています。
- １１〜１４行目は予測データを指定しています。
- １２行目は観測値の数（n=2）、予測変数の標準誤差に将来の観測値の分散が含まれるかどうか（vfobs=.TRUE.:含まれる）、観測値のオフセットが必要かどうか（offset='N':不要）、重みづけをするかどうか（weight='U':重みづけをしない）を指定しています。
- １３〜１４行目は観測値(x)を指定しています。
出力結果
（本ルーチンの詳細はg02gpf のマニュアルページを参照）| この出力例をダウンロード |
G02GPF Example Program Results
Residual sum of squares = 0.3872E+00
Degrees of freedom = 3
Estimate Standard error
-0.0239 0.0028
0.0638 0.0026
I ETA SE(ETA) Predicted SE(Predicted)
1) 2.01807 0.08168 0.49552 0.35981
2) 1.12472 0.04476 0.88911 0.36098
- ３行目に残差平方和と自由度が出力されています。
- ６〜９行目に一般化線形モデルの引数の推定値と標準誤差が出力されています。
- １１〜１４行目に線形予測子、線形予測子の標準誤差、予測値、予測値の標準誤差が出力されています。
ソースコード
（本ルーチンの詳細はg02gpf のマニュアルページを参照）
※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
PROGRAM g02gpfe
! G02GPF Example Program Text
! Mark 23 Release. nAG Copyright 2011.
! .. Use Statements ..
USE nag_library, ONLY : g02gaf, g02gpf, nag_wp
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Parameters ..
INTEGER, PARAMETER :: nin = 5, nout = 6
! .. Local Scalars ..
REAL (KIND=nag_wp) :: a, eps, rss, s, tol
INTEGER :: i, idf, ifail, ip, iprint, irank, &
ldv, ldx, loff, lt, lwk, lwt, m, &
maxit, n, tldx
LOGICAL :: vfobs
CHARACTER (1) :: error, link, mean, offset, weight
! .. Local Arrays ..
REAL (KIND=nag_wp), ALLOCATABLE :: b(:), cov(:), eta(:), off(:), &
pred(:), se(:), seeta(:), sepred(:), &
t(:), twt(:), tx(:,:), ty(:), &
v(:,:), wk(:), wt(:), x(:,:), y(:)
INTEGER, ALLOCATABLE :: isx(:)
! .. Intrinsic Functions ..
INTRINSIC count
! .. Executable Statements ..
WRITE (nout,*) 'G02GPF Example Program Results'
WRITE (nout,*)
! Skip headings in data file
READ (nin,*)
READ (nin,*)
! Read in training data for model that will be used for prediction
READ (nin,*) link, mean, offset, weight, n, m, s
IF (weight=='W' .OR. weight=='w') THEN
lwt = n
ELSE
lwt = 0
END IF
tldx = n
ALLOCATE (tx(tldx,m),ty(n),twt(lwt),isx(m))
! Read in data
IF (lwt>0) THEN
READ (nin,*) (tx(i,1:m),ty(i),twt(i),i=1,n)
ELSE
READ (nin,*) (tx(i,1:m),ty(i),i=1,n)
END IF
! Read in variable inclusion flags
READ (nin,*) isx(1:m)
! Calculate IP
ip = count(isx(1:m)>0)
IF (mean=='M' .OR. mean=='m') THEN
ip = ip + 1
END IF
! Read in power for exponential link
IF (link=='E' .OR. link=='e') THEN
READ (nin,*) a
END IF
ldv = n
lwk = (ip*ip+3*ip+22)/2
ALLOCATE (b(ip),se(ip),cov(ip*(ip+1)/2),v(ldv,ip+7),wk(lwk))
! Read in the offset
IF (offset=='Y' .OR. offset=='y') THEN
READ (nin,*) v(1:n,7)
END IF
! Read in control parameters
READ (nin,*) iprint, eps, tol, maxit
! Call routine to fit generalized linear model, with Normal errors,
! to training data
ifail = -1
CALL g02gaf(link,mean,offset,weight,n,tx,tldx,m,isx,ip,ty,twt,s,a,rss, &
idf,b,irank,se,cov,v,ldv,tol,maxit,iprint,eps,wk,ifail)
IF (ifail/=0) THEN
IF (ifail<6) THEN
GO TO 20
END IF
END IF
! Display parameter estimates for training data
WRITE (nout,99999) 'Residual sum of squares = ', rss
WRITE (nout,99998) 'Degrees of freedom = ', idf
WRITE (nout,*)
WRITE (nout,*) ' Estimate Standard error'
WRITE (nout,*)
WRITE (nout,99997) (b(i),se(i),i=1,ip)
! Skip second lot of headings in data file
READ (nin,*)
! Read in size of prediction data
READ (nin,*) n, vfobs, offset, weight
! Used G02GAF to fit training model, so error structure is normal is
! do not require T array
error = 'N'
lt = 0
ldx = n
IF (weight=='W' .OR. weight=='w') THEN
lwt = n
ELSE
lwt = 0
END IF
IF (offset=='Y' .OR. offset=='y') THEN
loff = n
ELSE
loff = 0
END IF
ALLOCATE (x(ldx,m),y(n),wt(lwt),off(loff),eta(n),seeta(n),pred(n), &
sepred(n),t(lt))
! Read in predication data
IF (lwt>0) THEN
READ (nin,*) (x(i,1:m),wt(i),i=1,n)
ELSE
READ (nin,*) (x(i,1:m),i=1,n)
END IF
! Read in offset for the prediction data
IF (offset=='Y' .OR. offset=='y') THEN
READ (nin,*) off(1:n)
END IF
! Call prediction routine
ifail = -1
CALL g02gpf(error,link,mean,offset,weight,n,x,ldx,m,isx,ip,t,off,wt,s, &
a,b,cov,vfobs,eta,seeta,pred,sepred,ifail)
IF (ifail/=0) THEN
WRITE (nout,99995) ifail
IF (ifail/=22) THEN
GO TO 20
END IF
END IF
! Display predicted values
WRITE (nout,*)
WRITE (nout,*) ' I ETA SE(ETA) Predicted ', &
' SE(Predicted)'
WRITE (nout,*)
WRITE (nout,99996) (i,eta(i),seeta(i),pred(i),sepred(i),i=1,n)
20 CONTINUE
99999 FORMAT (1X,A,E12.4)
99998 FORMAT (1X,A,I0)
99997 FORMAT (1X,2F14.4)
99996 FORMAT (1X,I3,')',F10.5,3F14.5)
99995 FORMAT (1X,A,I5)
END PROGRAM g02gpfe