R Time Series Analysis 時系列解析(3) (filterによる指数平滑化法)

R Program Mathematics Time Series Analysis 時系列解析 
####filter関数を使った指数平滑法####
#データは経済産業省総合原指数【月次】出荷（平成17年＝100.0）
#資本財を対象とする
dat <- read.csv("http://www.meti.go.jp/statistics/tyo/iip/result/h2afdldj/csv/ha2zom3j.csv",skip=2,header=TRUE) #skip=2は上位2行を省く意味
dat_ts <- ts(dat[,4],start=c(2003,1),frequency=12)
bunkai <- decompose(dat_ts,type="multiplicative")

dat_TC <- bunkai$trend
dat_TC2 <- dat_TC[is.na(dat_TC)==FALSE]
dat_TC3 <- ts(dat_TC2,start=c(2003,7),frequency=12)
n <- length(dat_TC3)
x <- rep(0,6);b <- 0;p <- 6
#xはダミー変数、bは予測値と実績値のバランス係数(0なら前時刻の実績重視、1なら前時刻の予測重視)
#pは過去何期を参照するかを示す
(w <- (1-b)*(cumprod(c(1,rep(b,p-1)))))
#移動平均ウェイト。指数平滑化法の場合、ウェイトが等比数列となる。

datINT <- rev(dat_TC3[(n-p+1):n]) #予測値の初期値を求めるため
Pre <- filter(x,filter=w,method="recursive",init=datINT)
Pre2 <- ts(Pre,start=c(2011,4),frequency=12)
title <- "指数平滑法による予測"
ts.plot(dat_TC3,Pre2,type="l",lty=c(1,2),main=title)

#漸化式による予測系列
#Y予測T = (1-β)YT-1+(1-β)βYT-2+…+(1-β)β^T-2*Y1において、
#Y予測T+1 - β×Y予測Tを計算すると
#Y予測T+1 = Y予測T +(1-β)×(YT-Y予測T)という関係式が導ける
#これは"予測値と実績値の差の1-β倍分を調整し、次の予測とする"という意味
#また、Y予測T+1 = (1-β)YT + βY予測T　…☆　とも書ける

b <- 0.8;w <- b　#bは上記のβにあたる定数
x <- (1-b)*dat_TC3[-1]　#xはトレンドの0.2倍値の系列
datHAT <- filter(x,filter=w,method="recursive",init=dat_TC3[1])
#☆式において、x：(1-β)YTに対応、w：βに対応、init：Y1に対応
#datHAT[1] == dat_TC3[1]*0.8 + x[1]
#datHAT[2] == datHAT[1]*0.8 + x[2] …のように計算が行われていることに注意

dathat <- ts(datHAT,start=c(2003,7),frequency=12)
title <- "漸化式による予測系列の作成"
ts.plot(dat_TC3,dathat,type="l",lty=c(1,2),main=title)

 f:id:nakhirot:20130614143023p:plain

【まとめ】

T期における実績値と予測値の内分点を、T+1期の予測とする方法が指数平滑化法であるとわかる。

Combining peak- and chromatogram-based retention time alignment algorithms for multiple chromatography-mass spectrometry data sets

Mathematics 論文 bioinformatics 遺伝統計学

論文はこれ

Background:

  • GC-MS,LC-MSデータのRetention Time(RT)のアラインメントは大まかに2つのカテゴリに分けられる
  • Peak-based algorithms：事前のpeak detectionに対してとてもsensitiveである。ピークモデルの形やSNR(signal-to-noise ratio)などのクライテリアとの一致でpeak detectする。
  • Raw data-based algorithms：uniform matrixと名付けられたbinned MS dataに基づいて実行。Uniform matrixはsignal mapsまたはbinned mass spectra間のペアワイズ類似度でDTW (Dynamic Time Warping) をするuniform matrix。多くのmass tracesを使うことはコンピュータリソースを必要とし、ノイズの傾向を増加させるかもしれない。

Methods:

  • この論文では、両者のhybrid methodを2つ紹介
  • BIPACEでよく対応するpeaks(Biderectional Best Hits)を抽出し、CEMAPP-DTWでglobal-alignmentを実行
  • この流れで実行

Best hIts Peak Assignment and Cluster Extesion(BIPACE)

  • Given a chromatogram C = {p1, p2, …, pl} as an ordered set of peaks
  • m(m/z), i(intensity)をvectorとし、RT = tのピークをp = (m,i,t)とする
  • つまり、RT毎に比較単位となるpを抽出している
  • 2つのchromatogram の各々のpeak p,qについてsimilarity function※を定義する
{f(p,q) = s(p,q) \exp(t_p-t_q)^{2} /2D^{2}, s(p,q) = \cos(arg(i_p,i_q))}
  • Dはretention time tolerance：peak matchingの対象となる領域を限定することで計算時間の節約
  • CとC'のピークの各々について、最も類似度の高いピークを抽出し、それをBBH (Biderectional Best Hits)とする
  • 次にBBHをmergeして集団(clique)にしていく
  • Minimal Clique Size(MCS)はパラメータ設定
  • peakを頂点、similality valueを辺としたグラフを作り、BBHをグラフ化する。complete subgraphとなるようにmergeする

CEnter-star Multiple Alignment by Pairwise Partitioned Dynamic Time Warping (CEMAPP-DTW)

  • Pairwise DTW is a global alignment of two series A=(a1,a2,…,am) and B=(b1,b2,…,bn), where ai, bi∈R^L
  • optimal alignmentを行うために(M+1)×(N+1)のmatrix Qを用意
  • Q[i, j]には(a1,a2,…,ai) and (b1,b2…,bj)のalignmentの最適なsimilarity valueが格納される
  • すると、このようなパスを作ることができる
  • 最適なパスは、pairwise similaritiesの和を最大化する。
  • つまり、DTW(A,B) := max (Σ pi∈P Q(pi)) (Pはパスの集合)
  • 計算時は、i=j=0から始めて、垂直方向に進むか、水平方向に進むか、対角線方向に進むのうち、最もsimilality valueが増加する方向を選択し、Q(M,N)まで辿り着けばそれがDTW(A,B)となる
  • パスに柔軟性を持たせるため、垂直方向、水平方向、対角線方向にはウェイトを設定する
  • symmetricなDTWの場合、これらのウェイトは過適合の問題を効果的に減らすためにも使われる