一変量時系列の平滑化標本スペクトル

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

Keyword: 一変量時系列, ラグウィンドウ, 平滑化標本スペクトル

概要

本サンプルはパルザンのラグウィンドウを用いた一変量時系列の平滑化標本スペクトルの計算を行うFortranによるサンプルプログラムです。 本サンプルは以下に示される時系列データを分析し、スペクトルを出力します。

時系列のデータ 

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

入力データ

(本ルーチンの詳細はg13caf のマニュアルページを参照)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

このデータをダウンロード
G13CAF Example Program Data
 256 100                                              :: NX,NC
1  0  0.1  4  100  200  0                             :: MTX,IC,PX,IW,MW,L,LG
360                                                   :: KC
  5.0  11.0  16.0  23.0  36.0  58.0  29.0  20.0  10.0   
  8.0   3.0   0.0   0.0   2.0  11.0  27.0  47.0  63.0
 60.0  39.0  28.0  26.0  22.0  11.0  21.0  40.0  78.0
122.0 103.0  73.0  47.0  35.0  11.0   5.0  16.0  34.0
 70.0  81.0 111.0 101.0  73.0  40.0  20.0  16.0   5.0
 11.0  22.0  40.0  60.0  80.9  83.4  47.7  47.8  30.7
 12.2   9.6  10.2  32.4  47.6  54.0  62.9  85.9  61.2
 45.1  36.4  20.9  11.4  37.8  69.8 106.1 100.8  81.6
 66.5  34.8  30.6   7.0  19.8  92.5 154.4 125.9  84.8
 68.1  38.5  22.8  10.2  24.1  82.9 132.0 130.9 118.1
 89.9  66.6  60.0  46.9  41.0  21.3  16.0   6.4   4.1
  6.8  14.5  34.0  45.0  43.1  47.5  42.2  28.1  10.1
  8.1   2.5   0.0   1.4   5.0  12.2  13.9  35.4  45.8
 41.1  30.1  23.9  15.6   6.6   4.0   1.8   8.5  16.6
 36.3  49.6  64.2  67.0  70.9  47.8  27.5   8.5  13.2
 56.9 121.5 138.3 103.2  85.7  64.6  36.7  24.2  10.7
 15.0  40.1  61.5  98.5 124.7  96.3  66.6  64.5  54.1
 39.0  20.6   6.7   4.3  22.7  54.8  93.8  95.8  77.2
 59.1  44.0  47.0  30.5  16.3   7.3  37.6  74.0 139.0
111.2 101.6  66.2  44.7  17.0  11.3  12.4   3.4   6.0
 32.3  54.3  59.7  63.7  63.5  52.2  25.4  13.1   6.8
  6.3   7.1  35.6  73.0  85.1  78.0  64.0  41.8  26.2
 26.7  12.1   9.5   2.7   5.0  24.4  42.0  63.5  53.8
 62.0  48.5  43.9  18.6   5.7   3.6   1.4   9.6  47.4
 57.1 103.9  80.6  63.6  37.6  26.1  14.2   5.8  16.7
 44.3  63.9  69.0  77.8  64.9  35.7  21.2  11.1   5.7
  8.7  36.1  79.7 114.4 109.6  88.8  67.8  47.5  30.6
 16.3   9.6  33.2  92.6 151.6 136.3 134.7  83.9  69.4
 31.5  13.9   4.4  38.0                               :: End of XG

  • 1行目はタイトル行で読み飛ばされます。
  • 2行目に時系列の長さ(nx)、計算される共分散の数(nc)を指定しています。
  • 3行目にデータが最初に切片の補正をされるのか傾きの補正をされるのか(mtx=1:切片)、共分散行列が計算されるのか与えられるのか(ic=0:計算される)、最初に減らされるデータの割合(px=0.1)、ラグウィンドウの選択(iw=4:Parzen)、ラグウィンドウのカットオフポイント(mw=100)、スペクトル推定の周波数分割(l=200)、ログありのスペクトル推定かログなしのスペクトル推定か(lg=0:ログなし)を指定しています。
  • 4行目に共分散の計算に使用されるFFT(高速フーリエ変換)の次数(kc)を指定しています。
  • 5~33行目に時系列データ(xg)を指定しています。

出力結果

(本ルーチンの詳細はg13caf のマニュアルページを参照)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

この出力例をダウンロード
 G13CAF Example Program Results

 Covariances
   1152.9733   937.3289   494.9243    14.8648  -342.8548  -514.6479
   -469.2733  -236.6896   109.0608   441.3498   637.4571   641.9954
    454.0505   154.5960  -136.8016  -343.3911  -421.8441  -374.4095
   -241.1943   -55.6140   129.4067   267.4248   311.8293   230.2807
     56.4402  -146.4689  -320.9948  -406.4077  -375.6384  -273.5936
   -132.6214    11.0791   126.4843   171.3391   122.6284   -11.5482
   -169.2623  -285.2358  -331.4567  -302.2945  -215.4832  -107.8732
     -3.4126    73.2521    98.0831    71.8949    17.0985   -27.5632
    -76.7900  -110.5354  -126.1383  -121.1043  -103.9362   -67.4619
    -10.8678    58.5009   116.4587   140.0961   129.5928    66.3211
    -35.5487  -135.3894  -203.7149  -216.2161  -152.7723   -30.4361
     99.3397   188.9594   204.9047   148.4056    34.4975  -103.7840
   -208.5982  -252.4128  -223.7600  -120.8640    23.3565   156.0956
    227.7642   228.5123   172.3820    87.4911   -21.2170  -117.5282
   -176.3634  -165.1218   -75.1308    67.1634   195.7290   279.3039
    290.8258   225.3811   104.0784   -44.4731  -162.7355  -207.7480
   -165.2444   -48.5473   118.8872   265.0045

 Degrees of freedom = 9.0      Bandwidth = 0.1165

 95 percent confidence limits -     Lower = 0.4731  Upper = 3.3329

       Spectrum       Spectrum      Spectrum       Spectrum
       estimate       estimate      estimate       estimate
    1  210.4696    2  428.2020    3  810.1419    4  922.5900
    5  706.1605    6  393.4052    7  207.6481    8  179.0657
    9  170.1320   10  133.0442   11  103.6752   12  103.0644
   13  141.5173   14  194.3041   15  266.5730   16  437.0181
   17  985.3130   18 2023.1574   19 2681.8980   20 2363.7439
   21 1669.9001   22 1012.1320   23  561.4822   24  467.2741
   25  441.9977   26  300.1985   27  172.0184   28  114.7823
   29   79.1533   30   49.4882   31   27.0902   32   16.8081
   33   27.5111   34   59.4429   35   97.0145   36  119.3664
   37  116.6737   38   87.3142   39   54.9570   40   42.9781
   41   46.6097   42   53.6206   43   50.6050   44   36.7780
   45   25.6285   46   24.8555   47   30.2626   48   31.5642
   49   27.3351   50   22.4443   51   18.5418   52   15.2425
   53   12.0207   54   12.6846   55   18.3975   56   19.3058
   57   12.6103   58    7.9511   59    7.1333   60    5.4996
   61    3.4182   62    3.2359   63    5.3836   64    8.5225
   65   10.0610   66    7.9483   67    4.2261   68    3.2631
   69    5.5751   70    7.8491   71    9.3694   72   11.0791
   73   10.1386   74    6.3158   75    3.6375   76    2.6561
   77    1.8026   78    1.0103   79    1.0693   80    2.3950
   81    4.0822   82    4.6221   83    4.0672   84    3.8460
   85    4.8489   86    6.3964   87    6.4762   88    4.9457
   89    4.4444   90    5.2131   91    5.0389   92    4.6141
   93    5.8722   94    7.9268   95    7.9486   96    5.7854
   97    4.5495   98    5.2696   99    6.3893  100    6.5216
  101    6.2129

  • 3~20行目に共分散が出力されています。
  • 22行目に自由度と平滑化ウィンドウの帯域幅が出力されています。
  • 24行目に95%信頼限界であることと、信頼限界の下限と上限が出力されています。
  • 26~53行目にはスペクトル推定値が出力されています。

ソースコード

(本ルーチンの詳細はg13caf のマニュアルページを参照)

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80

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

!      G13CAF Example Program Text

!      Mark 23 Release. nAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : g13caf, nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER              :: nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)              :: px
       INTEGER                         :: i, ic, ifail, iw, kc, l, lg, lxg,    &
                                          mtx, mw, nc, ng, nx, nxg
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE :: c(:), xg(:)
       REAL (KIND=nag_wp)              :: stats(4)
!      .. Intrinsic Functions ..
       INTRINSIC                          max
!      .. Executable Statements ..
       WRITE (nout,*) 'G13CAF Example Program Results'
       WRITE (nout,*)

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

!      Read in the problem size
       READ (nin,*) nx, nc

!      Read in smoothing parameters
       READ (nin,*) mtx, ic, px, iw, mw, l, lg
       IF (ic==0) THEN
          READ (nin,*) kc
       END IF

       IF (ic==0) THEN
          nxg = max(kc,l)
       ELSE
          nxg = l
       END IF
       lxg = max(nxg,nx)
       ALLOCATE (xg(lxg),c(nc))

!      Read in the data
       READ (nin,*) xg(1:nx)

!      Calculate smoothed spectrum
       ifail = -1
       CALL g13caf(nx,mtx,px,iw,mw,ic,nc,c,kc,l,lg,nxg,xg,ng,stats,ifail)
       IF (ifail/=0) THEN
          IF (ifail<4) THEN
             GO TO 20
          END IF
       END IF

!      Display results
       WRITE (nout,*) 'Covariances'
       WRITE (nout,99999) c(1:nc)
       WRITE (nout,*)
       WRITE (nout,99998) 'Degrees of freedom =', stats(1), &
          '      Bandwidth =', stats(4)
       WRITE (nout,*)
       WRITE (nout,99997) '95 percent confidence limits -     Lower =', &
          stats(2), '  Upper =', stats(3)
       WRITE (nout,*)
       WRITE (nout,*) &
          '      Spectrum       Spectrum      Spectrum       Spectrum'
       WRITE (nout,*) &
          '      estimate       estimate      estimate       estimate'
       WRITE (nout,99996) (i,xg(i),i=1,ng)

20     CONTINUE

99999  FORMAT (1X,6F11.4)
99998  FORMAT (1X,A,F4.1,A,F7.4)
99997  FORMAT (1X,A,F7.4,A,F7.4)
99996  FORMAT (1X,I4,F10.4,I5,F10.4,I5,F10.4,I5,F10.4)
    END PROGRAM g13cafe


関連情報
© 日本ニューメリカルアルゴリズムズグループ株式会社 2024
Privacy Policy  /  Trademarks