2009-3-23:時系列解析①　Rubyを用いて線形補間をする

1　Time Series Analysis（時系列解析）とは

　　時間観測されたデータは，連続的ではなく離散的なものです．そのため，このようなデータから時系列を推測するにはアルゴリズムが必要になります．これによって推測された波の集合である時系列を分解解析し，ある現象のメカニズムを解明していこうというものがTime Series Analysis（時系列解析）になります．時系列解析は，医学や経済の分野でも応用されているようです．
　　このワークショップでは，対象データとして，過去およそ42万年までさかのぼる氷床コアの酸素安定同位体δ18Oデータを使用します．このデータより，Rで時系列解析をおこなっていきます．

2　Interpolation

　　解析をする前に，離散的なデータを補間する必要があります．ここでは，２種類の補間法（interpolation）を紹介します.　　

1)　linear Interpolation

　　いろんな補間法があるなか，今回は一次関数のアルゴリズムを用いた補間（linear interpolation）をおこないます.以下の手順により，RとRubyで作業していきます.

a. Rで元のデータをプロットする

まず，エディタを立ち上げ，

(O:)ドライブ　→　研究室ワークショップ　→　DATA　→　ftp ftp.ncdc.noaa.gov pub data paleo･･･

より，データを落としてきてください．
ここで保存する際，ファイル名を「ice.txt」としておきます．

次に，Rを立ち上げ，以下のように打ってください．

ice <- read.table("ice（fileのパス）",sep = "\t",header = T)　…タブ区切りを分割，一行目のヘッダを認識して読み込んだものをiceと定義
ls()　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 …リストをとる（iceができているか確認）
ice　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　 …iceを表示
plot(ice)　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　…iceをプロット

そうすると，下のようなグラフができます．　　　　　　

つづけて，

lines(ice,type="l")

と打つと，下のグラフができます．

ここでtype="l"のlの部分をb，c，o，h，s，S，nと変えてみてください．色んな種類のグラフができます．
ここから，グラフに修正を加えていきます．

lines(ice,col="blue") 　　　　　…線に色をつける
lines(ice,col="red",type="p")　　　　…点に色をつける
plot(ice,type="o",col="blue")　　　　　　…type="o"で上書き
plot(ice,type="o",col="blue",xlim = c(420000,0))　　 …x軸を逆転する
plot(ice,type="o",col="blue",xlim = c(420000,0),ylim = c(1.5,-0.5))　　…y軸を逆転する

以上よりできあがったものがこちらのグラフです．

b. データをサンプリングする

　今回扱うデータは，時間軸が不等間隔なので改ざん式を使って等間隔になおしていきます.ここでは，下に示す簡単なアルゴリズムで，1000年ごとの物理量（ここではδ18Oの値）を補間していきます．
$(1)\hspace{3.5em}y_{t} = y_{i} + \frac{y_{i+1}-y_{i}}{x_{i+1}-x{i}} \times (x_{t}-x_{i})$

$(2)\hspace{3.5em}x_{t} \geq x_{i}$ 　　　 $x_{t} < x_{i+1}$

　(1)式は，元データの時間ｘ（i）と物理量ｙ（i）から1000年ごとの時間ｘ（ｔ）における物理量ｙ（ｔ）を求めるものになります．(2)式は，元データの不等間隔データに，1000年間隔の時間をあてはめるときに使用する式になります．これらの作業をRubyでおこないます．

エディタを立ち上げ，以下のようなスクリプトを作ってください．

################## DEFINITION #########################
file = "ice.txt"                           #file open
f = open(file)
f.gets                                     #一行読み飛ばす
x = []                                     #時間の配列の初期化（xの箱を作る）
y = []                                     #物理量の配列の初期化（yの箱を作る）
　　　　　　　　　　　　　　　　　　　　　　      #x = y = []でもよい 
while str = f.gets　　　　　　　　　　　　　　 #データがあるところまで読みとる
 v = str.chop.split("\t")                  #タブと改行コードを切り取る、vは二次元の配列
 x << v[0].to_i           　　　　　        　#xを文字列から整数に読み取る
 y << v[1].to_f                            #yを文字列から少数に読み取る
end
f.close                                    #開いたfileを閉じる

#p x                                       #xをとめる
#p y                                       #yをとめる
#__END__                                   #ここまでをとめる


 ################ LINEAR INTERPOLATION #################
xt = []; yt = []; n = x.size
(1..414).each do |k|                        #“i”を414回繰り返す
  xt[k] = 1000 * k                          #1000の倍数をさがす
  (n - 1).times do |i|                      #条件を満たすところをさがすまで(n-1)回繰り返す
    if xt[k] >= x[i] && xt[k] < x[i+1]      #もし式(2)なら
      yt[k] = y[i] + (y[i+1] - y[i]) / (x[i+1] - x[i]) * (xt[k] - x[i]) #式(1)の計算をする
      break
    end
  end
end

#p xt
#p yt


#################### OUTPUT FILE #######################
xt.shift; yt.shift                          #1列目の数字をとばす
p [xt.size,yt.size]                         #xt,ytの配列の中の数を表示
#__END__
out = open("ice1000","w")                   #ファイルを出力
out.puts "Age   D18O"                       #ヘッダ行
xt.each do |t|
  k = xt.index(t)                           #kの添え字は何番目になるかを書く
  out.printf "%d %.3f\n",t,yt[k]            #%d：数センチ,3f：少数点3桁,\n：改行,t：k番目の値がt
end
out.close

保存する際，出来上がったファイル名を「neq2eq.rb」としておきます．

コマンドプロンプトを立ち上げ，nec2eq.rbがあるディレクトリまで移動してください．
そこで　ruby nec2eq.rb でRubyを走らせてください．
上のRubyのスクリプトで作った「ice1000」というファイルができているか確認します．

c.　補間した値をプロットする

ここまでできたら，再びRでの作図になります．
上の“a)Rで元のデータをプロットする”で作ったグラフと同じものを，今度はice1000の補正データでつくります．
そうして，できたのが下のグラフになります．

今回は，ここまでになります．お疲れ様でした．

タグ：

+ タグ編集

「2009-3-23:時系列解析①　Rubyを用いて線形補間をする」をウィキ内検索

最終更新：2009年04月12日 19:14

添付ファイル

BioEoS::環境解析学研究室::

メニュー

更新履歴