Keyword: ポアソン誤差, 一般化線形モデル, フィット
概要
本サンプルはポアソン誤差をもつ一般化線形モデルのフィットを行うFortranによるサンプルプログラムです。 本サンプルは以下に示されるデータについて一般化線形モデルのフィットを行います。
※本サンプルはnAG Fortranライブラリに含まれるルーチン g02gcf() のExampleコードです。本サンプル及びルーチンの詳細情報は g02gcf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はg02gcf のマニュアルページを参照）| このデータをダウンロード |
G02GCF Example Program Data 'L' 'M' 'N' 'U' 15 8 :: LINK,MEAN,OFFSET,WEIGHT,N,M 1.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 141.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 67.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 114.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 79.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 39.0 0.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 131.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 0.0 66.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 0.0 143.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 72.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 35.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 0.0 36.0 0.0 0.0 1.0 0.0 1.0 0.0 0.0 0.0 14.0 0.0 0.0 1.0 0.0 0.0 1.0 0.0 0.0 38.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0 0.0 28.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 16.0 :: End of X, Y 1 1 1 1 1 1 1 1 :: ISX 0 1.0E-6 5.0E-5 0 :: IPRINT,EPS,TOL,MAXIT
- １行目はタイトル行で読み飛ばされます。
- ２行目はどのリンク関数が使用されるか（link='L':ログリンク）、一般化線形モデルの引数に切片が含まれるかどうか（mean='M':切片を含める）、オフセットが必要かどうか（offset='N':不要）、重みづけをするかどうか（weight='U':重みづけをしない）、観測値の数（n=15）、独立変数の数（m=8）を指定しています。
- ３〜１７行目に独立変数の観測値(x)と従属変数の観測値(y)を指定しています。
- １８行目はモデルに独立変数が含まれるかどうか(isx)を指定しています。"1"は含まれることを意味します。
- １９行目は反復の情報の出力が必要かどうか（iprint=0:不要）、独立変数がフルランクかどうかそうでない場合ランクは何か（eps=1.0E-6）、モデルのフィットに必要な正確さ（tol=5.0E-5）、最大反復数（maxit=0:デフォルト値10を使用）を指定しています。
出力結果
（本ルーチンの詳細はg02gcf のマニュアルページを参照）| この出力例をダウンロード |
G02GCF Example Program Results
Deviance = 0.9038E+01
Degrees of freedom = 8
Estimate Standard error
2.5977 0.0258
1.2619 0.0438
1.2777 0.0436
0.0580 0.0668
1.0307 0.0551
0.2910 0.0732
0.9876 0.0559
0.4880 0.0675
-0.1996 0.0904
Y FV Residual H
141.0 132.99 0.6875 0.604
67.0 63.47 0.4386 0.514
114.0 127.38 -1.2072 0.596
79.0 77.29 0.1936 0.532
39.0 38.86 0.0222 0.482
131.0 135.11 -0.3553 0.608
66.0 64.48 0.1881 0.520
143.0 129.41 1.1749 0.601
72.0 78.52 -0.7465 0.537
35.0 39.48 -0.7271 0.488
36.0 39.90 -0.6276 0.393
14.0 19.04 -1.2131 0.255
38.0 38.21 -0.0346 0.382
28.0 23.19 0.9675 0.282
16.0 11.66 1.2028 0.206
- ３行目にデビアンス（deviance）が出力されています。
- ４行目に自由度が出力されています。
- ６〜１６行目に一般化線形モデルの引数の推定値と標準誤差が出力されています。
- １８〜３４行目に従属変数の観測値、フィットされた値、残差、てこ値が出力されています。
ソースコード
（本ルーチンの詳細はg02gcf のマニュアルページを参照）
※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
PROGRAM g02gcfe
! G02GCF Example Program Text
! Mark 23 Release. nAG Copyright 2011.
! .. Use Statements ..
USE nag_library, ONLY : g02gcf, nag_wp
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Parameters ..
INTEGER, PARAMETER :: nin = 5, nout = 6
! .. Local Scalars ..
REAL (KIND=nag_wp) :: a, dev, eps, tol
INTEGER :: i, idf, ifail, ip, iprint, irank, &
ldv, ldx, lwk, lwt, m, maxit, n
CHARACTER (1) :: link, mean, offset, weight
! .. Local Arrays ..
REAL (KIND=nag_wp), ALLOCATABLE :: b(:), cov(:), se(:), v(:,:), wk(:), &
wt(:), x(:,:), y(:)
INTEGER, ALLOCATABLE :: isx(:)
! .. Intrinsic Functions ..
INTRINSIC count
! .. Executable Statements ..
WRITE (nout,*) 'G02GCF Example Program Results'
WRITE (nout,*)
! Skip heading in data file
READ (nin,*)
! Read in the problem size
READ (nin,*) link, mean, offset, weight, n, m
IF (weight=='W' .OR. weight=='w') THEN
lwt = n
ELSE
lwt = 0
END IF
ldx = n
ALLOCATE (x(ldx,m),y(n),wt(lwt),isx(m))
! Read in data
IF (lwt>0) THEN
READ (nin,*) (x(i,1:m),y(i),wt(i),i=1,n)
ELSE
READ (nin,*) (x(i,1:m),y(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
! Fit generalized linear model with Poisson errors
ifail = -1
CALL g02gcf(link,mean,offset,weight,n,x,ldx,m,isx,ip,y,wt,a,dev,idf,b, &
irank,se,cov,v,ldv,tol,maxit,iprint,eps,wk,ifail)
IF (ifail/=0) THEN
IF (ifail<7) THEN
GO TO 20
END IF
END IF
! Display results
WRITE (nout,99999) 'Deviance = ', dev
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)
WRITE (nout,*)
WRITE (nout,*) ' Y FV Residual H'
WRITE (nout,*)
WRITE (nout,99996) (y(i),v(i,2),v(i,5),v(i,6),i=1,n)
20 CONTINUE
99999 FORMAT (1X,A,E12.4)
99998 FORMAT (1X,A,I0)
99997 FORMAT (1X,2F14.4)
99996 FORMAT (1X,F7.1,F10.2,F12.4,F10.3)
END PROGRAM g02gcfe