相関行列のロバスト推定の計算

ホーム > 製品 > nAG数値計算ライブラリ > nAG Fortranライブラリ > サンプルソースコード集 > 相関行列のロバスト推定の計算

Keyword: 相関行列, ロバスト推定

概要

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

ロバスト推定のデータ　

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

入力データ

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

G02HMF Example Program Data
    10    3                   : N,M
  3.4  6.9  12.2              
  6.4  2.5  15.1
  4.9  5.5  14.2
  7.3  1.9  18.2
  8.8  3.6  11.7
  8.4  1.3  17.9
  5.3  3.1  15.0
  2.7  8.1   7.7
  6.1  3.0  21.9
  5.3  2.2  13.9              : End of X
  1.0 0.0 1.0 0.0 0.0 1.0     : A
  0.0 0.0 0.0                 : THETA
  4.0 2.0                     : RUSER
1 0 0.9 0.9 50 5.0E-5         : INDM,NITMON,BL,BD,MAXIT,TOL

１行目はタイトル行で読み飛ばされます。
２行目は観測値の数（n=10）と独立変数の数（m=3）を指定しています。
３～１２行目に独立変数の観測値(x)を指定しています。
１３行目は下三角行列Aの初期推定値を指定しています。
１４行目は位置パラメータθの初期推定値を指定しています。
１５行目はサブルーチンUCVのパラメータ(ruser)を指定しています。
１６行目は使用される関数vの形式（indm=1:v=1）、出力される反復の情報量（nitmon=0:出力なし）、非対角要素の境界の大きさ（bl=0.9）、対角要素の境界の大きさ（bd=0.9）、最大反復数（maxit=50）、共分散行列の最終推定値の相対精度（tol=5.0E-5）を指定しています。

出力結果

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

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

 G02HMF Example Program Results

 G02HMF required 34 iterations to converge

 Robust covariance matrix
             1          2          3
 1      3.2779    -3.6918     4.7391
 2                 5.2841    -6.4087
 3                           11.8373

 Robust estimates of THETA
      5.700
      3.864
     14.704

３行目には収束するのに34回の反復が必要だったことが示されています。
５～９行目にロバスト共分散行列が出力されています。
１１～１４行目に位置パラメータθのロバスト推定値が出力されています。

ソースコード

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

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

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

!   G02HMF Example Program Text
!   Mark 23 Release. nAG Copyright 2011.

    MODULE g02hmfe_mod

!      G02HMF Example Program Module:
!             Parameters and User-defined Routines

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: iset = 1, nin = 5, nout = 6
    CONTAINS
       SUBROUTINE ucv(t,ruser,u,w)

!         u function

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: t
          REAL (KIND=nag_wp), INTENT (OUT)    :: u, w
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (INOUT)  :: ruser(*)
!         .. Local Scalars ..
          REAL (KIND=nag_wp)                  :: cu, cw, t2
!         .. Executable Statements ..
          cu = ruser(1)
          u = 1.0_nag_wp
          IF (t/=0.0_nag_wp) THEN
             t2 = t*t
             IF (t2>cu) THEN
                u = cu/t2
             END IF
          END IF
!         w function
          cw = ruser(2)
          IF (t>cw) THEN
             w = cw/t
          ELSE
             w = 1.0_nag_wp
          END IF
       END SUBROUTINE ucv
    END MODULE g02hmfe_mod
    PROGRAM g02hmfe

!      G02HMF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : g02hmf, nag_wp, x04abf, x04ccf
       USE g02hmfe_mod, ONLY : iset, nin, nout, ucv
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: bd, bl, tol
       INTEGER                             :: i, ifail, indm, la, lcov, ldx,   &
                                              lruser, m, maxit, n, nadv, nit,  &
                                              nitmon
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: a(:), cov(:), ruser(:),          &
                                              theta(:), wk(:), wt(:), x(:,:)
!      .. Executable Statements ..
       WRITE (nout,*) 'G02HMF Example Program Results'
       WRITE (nout,*)

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

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

       ldx = n
       lruser = 2
       la = ((m+1)*m)/2
       lcov = la
       ALLOCATE (x(ldx,m),ruser(lruser),cov(lcov),a(la),wt(n),theta(m), &
          wk(2*m))

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

!      Read in the initial value of A
       READ (nin,*) a(1:la)

!      Read in the initial value of THETA
       READ (nin,*) theta(1:m)

!      Read in the values of the parameters of the ucv functions
       READ (nin,*) ruser(1:lruser)

!      Read in the control parameters
       READ (nin,*) indm, nitmon, bl, bd, maxit, tol

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

!      Compute robust estimate of variance / covariance matrix
       ifail = 0
       CALL g02hmf(ucv,ruser,indm,n,m,x,ldx,cov,a,wt,theta,bl,bd,maxit,nitmon, &
          tol,nit,wk,ifail)

!      Display results
       WRITE (nout,99999) 'G02HMF required ', nit, ' iterations to converge'
       WRITE (nout,*)
       FLUSH (nout)
       ifail = 0
       CALL x04ccf('Upper','Non-Unit',m,cov,'Robust covariance matrix',ifail)
       WRITE (nout,*)
       WRITE (nout,*) 'Robust estimates of THETA'
       WRITE (nout,99998) theta(1:m)

99999  FORMAT (1X,A,I0,A)
99998  FORMAT (1X,F10.3)
    END PROGRAM g02hmfe

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

概要

入力データ

出力結果

ソースコード