ロバスト回帰のM推定量の計算

ホーム > 製品 > nAG数値計算ライブラリ > nAG Fortranライブラリ > サンプルソースコード集 > ロバスト回帰のM推定量の計算

Keyword: ロバスト回帰, M推定量

概要

本サンプルはロバスト回帰のM推定量の計算を行うFortranによるサンプルプログラムです。本サンプルは以下に示されるデータについてロバスト回帰のM推定量の計算を行います。

ロバスト回帰のデータ　

※本サンプルはnAG Fortranライブラリに含まれるルーチン g02haf() のExampleコードです。本サンプル及びルーチンの詳細情報は g02haf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで

入力データ

このデータをダウンロード

G02HAF Example Program Data
 8  3                       : N,M
 1. -1. -1. 2.1  
 1. -1.  1. 3.6
 1.  1. -1. 4.5
 1.  1.  1. 6.1
 1. -2.  0. 1.3
 1.  0. -2. 1.9
 1.  2.  0. 6.7
 1.  0.  2. 5.5             : End of X1 X2 X3 and Y values
 1  2  1 0 50 1.0E-5        : INDW,IPSI,ISIGMA,NITMON,MAXIT,TOL
 3.0  0                     : CUCV,INDC
 1.5 3.0 4.5                : H1,H2,H3
 1.5                        : DCHI
 1.0                        : Initial value for SIGMA
 0.0 0.0 0.0                : Initial values for THETA

１行目はタイトル行で読み飛ばされます。
２行目は観測値の数(n)と独立変数の数(m)を指定しています。
３～１０行目に独立変数の観測値(x)と従属変数の観測値(y)を指定しています。
１１行目は実行される回帰の種類（indw=1:Schweppeタイプの回帰）、どのΨ関数が使用されるか（ipsi=2:Hampelの区分線形関数）、σの推定方法（isigma=1:χ関数を使用して推定）, 出力される反復の情報量（nitmon=0:出力なし）、最大反復数（maxit=50）、Aの計算の相対精度（tol=1.0E-5）を指定しています。
１２行目はKrasker-Welschの重みのu関数の定数（cucv=3.0）、共分散行列の推定で使用される近似法（indc=0:期待値を観測値で置き換え）を指定しています。
１３行目はHampelの区分線形関数のパラメータ(h1,h2,h3)を指定しています。
１４行目はχ関数の定数(dchi)を指定しています。
１５行目はσの初期値(sigma)を指定しています。
１６行目はθの初期値(theta)を指定しています。

出力結果

（本ルーチンの詳細はg02haf のマニュアルページを参照）

この出力例をダウンロード

 G02HAF Example Program Results

 Sigma =     0.2026

        THETA      Standard
                    errors
       4.0423       0.0384
       1.3083       0.0272
       0.7519       0.0311

      Weights     Residuals
       0.5783       0.1179
       0.5783       0.1141
       0.5783      -0.0987
       0.5783      -0.0026
       0.4603      -0.1256
       0.4603      -0.6385
       0.4603       0.0410
       0.4603      -0.0462

３行目に σ の最終推定値が出力されています。
７～９行目にθと標準誤差が出力されています。
１２～１９行目に各観測値の重みと残差が出力されています。

ソースコード

（本ルーチンの詳細はg02haf のマニュアルページを参照）

※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法

このソースコードをダウンロード

    PROGRAM g02hafe

!      G02HAF Example Program Text

!      Mark 23 Release. nAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : g02haf, nag_wp, x04abf
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: iset = 1, nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: cpsi, cucv, dchi, h1, h2, h3, sigma, &
                                          tol
       INTEGER                         :: i, ifail, indc, indw, ipsi, isigma,  &
                                          ldc, ldx, lwork, m, maxit, n, nadv,  &
                                          nitmon
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: c(:,:), rs(:), theta(:), wgt(:),     &
                                          work(:), x(:,:), y(:)
!      .. Executable Statements ..
       WRITE (nout,*) 'G02HAF Example Program Results'
       WRITE (nout,*)

!      Skip heading in data file
       READ (nin,*)

!      Read in the problem size
       READ (nin,*) n, m

       ldx = n
       ldc = m
       lwork = 4*n + m*(n+m)
       ALLOCATE (x(ldx,m),y(n),theta(m),c(ldc,m),work(lwork),rs(n),wgt(n))

!      Read in data
       READ (nin,*) (x(i,1:m),y(i),i=1,n)

!      Read in control parameters
       READ (nin,*) indw, ipsi, isigma, nitmon, maxit, tol

!      Read in appropriate weight function parameters
       IF (indw/=0) THEN
          READ (nin,*) cucv, indc
       END IF
       IF (ipsi>0) THEN
          IF (ipsi==1) THEN
             READ (nin,*) cpsi
          ELSE IF (ipsi==2) THEN
             READ (nin,*) h1, h2, h3
          END IF
          IF (isigma>0) THEN
             READ (nin,*) dchi
          END IF
       END IF

!      Set the advisory channel to NOUT for monitoring information
       IF (nitmon/=0) THEN
          nadv = nout
          CALL x04abf(iset,nadv)
       END IF

!      Read in initial values
       READ (nin,*) sigma
       READ (nin,*) theta(1:m)

!      Perform M-estimate regression
       ifail = -1
       CALL g02haf(indw,ipsi,isigma,indc,n,m,x,ldx,y,cpsi,h1,h2,h3,cucv,dchi, &
          theta,sigma,c,ldc,rs,wgt,tol,maxit,nitmon,work,ifail)
       IF (ifail/=0) THEN
          IF (ifail<7) THEN
             GO TO 20
          ELSE
             WRITE (nout,*) &
                '       Some of the following reslts may be unreliable'
          END IF
       END IF

!      Display results
       WRITE (nout,99999) 'Sigma = ', sigma
       WRITE (nout,*)
       WRITE (nout,*) '       THETA      Standard'
       WRITE (nout,*) '                   errors'
       WRITE (nout,99998) (theta(i),c(i,i),i=1,m)
       WRITE (nout,*)
       WRITE (nout,*) '     Weights     Residuals'
       WRITE (nout,99998) (wgt(i),rs(i),i=1,n)

20     CONTINUE

99999  FORMAT (1X,A,F10.4)
99998  FORMAT (1X,F12.4,F13.4)
    END PROGRAM g02hafe

Fortranによるサンプルソースコード : 使用ルーチン名：g02haf

概要

入力データ

出力結果

ソースコード