Keyword: 重回帰分析
概要
本サンプルは重回帰分析を行うFortranによるサンプルプログラムです。 本サンプルは以下に示される独立変数の観測値と従属変数の観測値を用いて重回帰分析を行い、独立変数の偏回帰係数のパラメータ推定と標準誤差、観測値の残差と対角要素のてこ値を出力します。
※本サンプルはnAG Fortranライブラリに含まれるルーチン g02daf() のExampleコードです。本サンプル及びルーチンの詳細情報は g02daf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はg02daf のマニュアルページを参照）| このデータをダウンロード |
G02DAF Example Program Data 12 4 'U' 'M' 1.0 0.0 0.0 0.0 33.63 0.0 0.0 0.0 1.0 39.62 0.0 1.0 0.0 0.0 38.18 0.0 0.0 1.0 0.0 41.46 0.0 0.0 0.0 1.0 38.02 0.0 1.0 0.0 0.0 35.83 0.0 0.0 0.0 1.0 35.99 1.0 0.0 0.0 0.0 36.58 0.0 0.0 1.0 0.0 42.92 1.0 0.0 0.0 0.0 37.80 0.0 0.0 1.0 0.0 40.43 0.0 1.0 0.0 0.0 37.89 1 1 1 1
- １行目はタイトル行で読み飛ばされます。
- ２行目に観測値の数（n=12）、独立変数の数（m=4）、重みづけをするかどうかを示すパラメータ（weight='U'）、切片を含むかどうかを示すパラメータ（mean='M'）を指定しています。"U"は最小二乗推定が使用されることを意味しており、"M"はモデルに切片が含まれることを意味しています。
- ３〜１４行目には左から縦１列から縦４列までに独立変数の観測値(x)を、最後の列に従属変数の観測値(y)を指定しています。
- １５行目にどの独立変数がモデルに含まれるかを示すフラグ(isx)を指定しています。"1"は含まれることを意味します。
出力結果
（本ルーチンの詳細はg02daf のマニュアルページを参照）| この出力例をダウンロード |
G02DAF Example Program Results
Model not of full rank, rank = 4
Residual sum of squares = 0.2223E+02
Degrees of freedom = 8
R-squared = 0.7004E+00
Adjusted R-squared = 0.5881E+00
AIC = 0.1540E+02
Variable Parameter estimate Standard error
1 0.3056E+02 0.3849E+00
2 0.5447E+01 0.8390E+00
3 0.6743E+01 0.8390E+00
4 0.1105E+02 0.8390E+00
5 0.7320E+01 0.8390E+00
Obs Residuals H
1 -0.2373E+01 0.3333E+00
2 0.1743E+01 0.3333E+00
3 0.8800E+00 0.3333E+00
4 -0.1433E+00 0.3333E+00
5 0.1433E+00 0.3333E+00
6 -0.1470E+01 0.3333E+00
7 -0.1887E+01 0.3333E+00
8 0.5767E+00 0.3333E+00
9 0.1317E+01 0.3333E+00
10 0.1797E+01 0.3333E+00
11 -0.1173E+01 0.3333E+00
12 0.5900E+00 0.3333E+00
- ３行目に回帰モデルのランクが出力されています。
- ５行目に残差平方和が出力されています。
- ６行目に自由度が出力されています。
- ７行目に決定係数が出力されています。
- ８行目に調整済み決定係数が出力されています。
- ９行目にAIC（Akaike's informaiton criteria:赤池情報量基準）が出力されています。
- １３〜１７行目に各変数のパラメータ推定（偏回帰係数）と標準誤差が出力されています。
- ２１〜３２行目に各観測値の残差と対角要素のてこ値が出力されています。
ソースコード
（本ルーチンの詳細はg02daf のマニュアルページを参照）
※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
PROGRAM g02dafe
! G02DAF Example Program Text
! Mark 23 Release. nAG Copyright 2011.
! .. Use Statements ..
USE nag_library, ONLY : g02buf, g02daf, nag_wp
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Parameters ..
INTEGER, PARAMETER :: nin = 5, nout = 6
! .. Local Scalars ..
REAL (KIND=nag_wp) :: aic, arsq, en, mult, rsq, rss, sw, tol
INTEGER :: i, idf, ifail, ip, irank, ldq, ldx, &
lwt, m, n
LOGICAL :: svd
CHARACTER (1) :: mean, weight
! .. Local Arrays ..
REAL (KIND=nag_wp), ALLOCATABLE :: b(:), cov(:), h(:), p(:), q(:,:), &
res(:), se(:), wk(:), wt(:), x(:,:), &
y(:)
REAL (KIND=nag_wp) :: c(1), wmean(1)
INTEGER, ALLOCATABLE :: isx(:)
! .. Intrinsic Functions ..
INTRINSIC count, log, real
! .. Executable Statements ..
WRITE (nout,*) 'G02DAF Example Program Results'
WRITE (nout,*)
! Skip heading in data file
READ (nin,*)
! Read in the problem size
READ (nin,*) n, m, weight, mean
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
ldq = n
ALLOCATE (b(ip),cov((ip*ip+ip)/2),h(n),p(ip*(ip+ &
2)),q(ldq,ip+1),res(n),se(ip),wk(ip*ip+5*(ip-1)))
! Use suggested value for tolerance
tol = 0.000001E0_nag_wp
! Fit general linear regression model
ifail = -1
CALL g02daf(mean,weight,n,x,ldx,m,isx,ip,y,wt,rss,idf,b,se,cov,res,h,q, &
ldq,svd,irank,p,tol,wk,ifail)
IF (ifail/=0) THEN
IF (ifail/=5) THEN
GO TO 20
END IF
END IF
! Calculate (weighted) total sums of squares, adjusted for mean if required
! If in G02DAF, an intercept is added to the regression by including a
! column of 1's in X, rather than by using the MEAN argument then
! MEAN = 'M' should be used in this call to G02BUF.
ifail = 0
CALL g02buf(mean,weight,n,1,y,n,wt,sw,wmean,c,ifail)
! Get effective number of observations (=N if there are no zero weights)
en = real(idf+irank,kind=nag_wp)
! Calculate R-squared, corrected R-squared and AIC
rsq = 1.0_nag_wp - rss/c(1)
IF (mean=='M' .OR. mean=='m') THEN
mult = (en-1.0E0_nag_wp)/(en-real(irank,kind=nag_wp))
ELSE
mult = en/(en-real(irank,kind=nag_wp))
END IF
arsq = 1.0_nag_wp - mult*(1.0_nag_wp-rsq)
aic = en*log(rss/en) + 2.0_nag_wp*real(irank,kind=nag_wp)
! Display results
IF (svd) THEN
WRITE (nout,99999) 'Model not of full rank, rank = ', irank
WRITE (nout,*)
END IF
WRITE (nout,99998) 'Residual sum of squares = ', rss
WRITE (nout,99999) 'Degrees of freedom = ', idf
WRITE (nout,99998) 'R-squared = ', rsq
WRITE (nout,99998) 'Adjusted R-squared = ', arsq
WRITE (nout,99998) 'AIC = ', aic
WRITE (nout,*)
WRITE (nout,*) 'Variable Parameter estimate ', 'Standard error'
WRITE (nout,*)
IF (ifail==0) THEN
WRITE (nout,99997) (i,b(i),se(i),i=1,ip)
ELSE
WRITE (nout,99996) (i,b(i),i=1,ip)
END IF
WRITE (nout,*)
WRITE (nout,*) ' Obs Residuals H'
WRITE (nout,*)
WRITE (nout,99997) (i,res(i),h(i),i=1,n)
20 CONTINUE
99999 FORMAT (1X,A,I4)
99998 FORMAT (1X,A,E12.4)
99997 FORMAT (1X,I6,2E20.4)
99996 FORMAT (1X,I6,E20.4)
END PROGRAM g02dafe