遺伝学・遺伝統計学関連の姉妹ブログ『ryamadaの遺伝学・遺伝統計学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

パターン検出

自己相関
  • 上で書いたこと
  • 時間の流れの中で、繰り返しパターンがあるかどうかは、ごくごく単純に考えれば、「繰り返しのタイムラグ」を適当に選びながら、「重ね合わせ」がどの程度よいか、で判断できる
  • Wiki記事
  • あとは、その軸の取り替え・軸の多数化・バリエーションで、ある程度のことはできる(はず)
  • Rの時系列解析
  • Rの回帰一般
  • Moran's spacial autocorrelation coefficient
  • 上記「Rの時系列解析　資料１−２」は、株取引データを扱っているようなので、それを、作ってみる(６カラムが日付、取引開始値、日最大値、日最小値、取引終了値、取引量。データとしての不適切などはこの際目をつぶるとする)
ndate<-100
ncol<-6
tui<-matrix(0,ndate,ncol)
tui[,1]<-rnorm(ndate)
tui[,1]<-1:ndate
tui[,2]<-sort(rnorm(ndate))
tui[,3]<-tui[,2]+rnorm(ndate)
tui[,3]<-tui[,2]+rnorm(ndate)+10
tui[,4]<-tui[,2]+rnorm(ndate)-10
tui[,5]<-tui[,2]+rnorm(ndate)*0.5
tui[,6]<-runif(ndate)

#まず値、の時系列をプロットする
plot(tui[,5],type="l")
# 増減を見る
plot(diff(tui[,5]),type="l")
  • 和の関数"sum",積の関数"prod"はそれぞれ二項演算子だけれど、差の関数"diff"はベクトルに対して作用することを前提に作られた関数なのは、ここのように使うことをデフォルトとして考えているというのが理由かもしれない。
  • 「差分」の関数→微分的に評価。それを確率分布に使って確かめれば
x<-seq(from=0,to=10,by=0.1)
q<-pexp(x,1)
d<-dexp(x,1)
dif<-diff(q)
par(mfcol=c(1,3))
plot(q,type="l")
plot(d,type="l")
plot(dif,type="l")
  • 差分の分布を評価する
    • ヒストグラムを描く
      • ヒストグラムは総標本数に対する割合でプロットすると、density()関数や、確率分布関数との重ね合わせに便利なので、そのために"prob=TRUE"のオプションを用いている
h<-hist(diff(tui[,4]),prob=T)

hist(diff(tui[,4]),prob=T)
lines(density(diff(tui[,4])))

mu<-mean(diff(tui[,4]))
sigma<-sd(diff(tui[,4]))

x<-seq(min(h$breaks),max(h$breaks),length=100)
lines(x,dnorm(x,mean=mu,sd=sigma),col="blue")
  • 上記作業は、分布を描いて、それを「正規分布だったら」という分布と重ね合わせて比較する作業
  • 別のやり方でこれをやる
    • 正規分布の累積確率がy=xの直線になるようなスケールで描いて、y=xの直線を付け加える処理
qqnorm(diff(tui[,4]))
qqline(y, col = 2)
  • 分布が正規分布様だと思えたら、それを検定したくなる
x<-diff(tui[,4])
# 分布xを累積分布関数(ここでは正規分布)と比較している
# "pnorm"には累積分布関数名を入れる。それに続く引数は、累積分布関数を規定する引数
ks.test(x,"pnorm",mean(x),sd(x))
# Shapiro-Wilk Normality Testは正規分布に近いかどうかのテスト
shapiro.test(x)

> ks.test(x,"pnorm",mean(x),sd(x))

        One-sample Kolmogorov-Smirnov test

data:  x 
D = 0.0467, p-value = 0.975
alternative hypothesis: two-sided 

> shapiro.test(x)

        Shapiro-Wilk normality test

data:  x 
W = 0.995, p-value = 0.976
  • 時間の変化を追う(株式ニュースなどでよく見るグラフ)
    • 均す(ならす）
      • 指定の幅の平均をとる
      • 均す幅が広くなると、グラフの両端のプロットのない部分が長くなることに注意
plot(tui[,5],type="l")
# 1/5の5は2+1+2なので、自分を中心に前後2ずつ、合わせて5つの値の平均を取って、その値を、自身の位置の値に置きなおす、という処理
# 以前はtsというパッケージを読み込む必要があったようだが、今は標準装備されているらしい
tui.1<-filter(tui[,5],filter=rep(1/5,5))
tui.2<-filter(tui[,5],filter=rep(1/25,25))
tui.3<-filter(tui[,5],filter=rep(1/50,50))
lines(tui.1,col="red")
lines(tui.2,col="blue")
lines(tui.3,col="green")
  • 周期性を見よう
    • ts()関数とstl()関数
    • ts()関数は、時系列データを周期的な時系列(毎月ｘ何年)に変える
    • stl()関数は、１年間の変動繰り返し要素と、全体としての変化と、それらで説明できない要素に分ける
    • それをプロットすると、「生プロット」「年周期」「全体の動き」「残差」の４つを縦に並べてプロットする
tt<-ts(tui[,5],freq=12)
plot(stl(tt,s.window="periodic"))