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

Time series handling

R Time series
  • Time series 時系列？
    • Pair of Time & Values 時刻と値のペア
      • Fixed interval vs. Varied interval 固定間隔 対 不定間隔
  • Make a time series data. 時系列データ
# time
n.t <- 100
my.time <- 1:n.t
# There are three variables x,y
phi <- 0.1
x <- cos(phi * my.time)
y <- sin(phi * my.time)
par(mfcol=c(1,2))
plot(x,y)
matplot(cbind(x,y),type="l")
par(mfcol=c(1,1))
  • A bit more complicated 少し複雑に
    • Stock market　株価
ndate<-100
ncol<-3
tui<-matrix(0,ndate,ncol)
# Fixed-interval time　固定間隔時刻
tui[,1]<-1:ndate 
# Monotonic increase but variable 単調増減成分、ただし値はばらばら
tui[,2]<-sort(rnorm(ndate,2,1))
plot(tui[,1],tui[,2],type="l")
# Randoms 乱雑項
tui[,3]<-tui[,2]+rnorm(ndate)*0.5
plot(tui[,1],tui[,3],type="l")
  • In time series analysis, you are interested in the increase/decrease. 時系列解析では、増減に興味がある
plot(diff(tui[,3]),type="l")
  • difference is following normal distribution?? 増減は正規分布？
hist(diff(tui[,3]),prob=T)
lines(density(diff(tui[,3])))
# mean ,sd
mu<-mean(diff(tui[,3]))
sigma<-sd(diff(tui[,3]))
# Normal distribution with mu and sigma
x<-seq(min(h$breaks),max(h$breaks),length=100)
lines(x,dnorm(x,mean=mu,sd=sigma),col="blue")
  • test it
# data is in the shape of normal dist?
ks.test(diff(tui[,3]),"pnorm",mu,sigma)
# Another way
# Shapiro-Wilk Normality Test specific for normal dist
shapiro.test(diff(tui[,3]))
  • Smoothing 
plot(tui[,3],type="l")
# filter represent how to smoothe the values
# 5 in 1/5 means 5 = 2+1+2-> (x[(i-2):(i+2)]/5
tui.1<-filter(tui[,3],filter=rep(1/5,5))
tui.2<-filter(tui[,3],filter=rep(1/25,25))
tui.3<-filter(tui[,3],filter=rep(1/50,50))
lines(tui.1,col="red")
lines(tui.2,col="blue")
lines(tui.3,col="green")
  • Decompose into periodic component and trend and residual
# 12 is twelve months
# Four panels are (1) raw (2) annual (every 12) (3) without annual (4) residual (1)「生プロット」(2)「年周期」(3)「全体の動き」(4)「残差」の
tt<-ts(tui[,3],freq=12)
plot(stl(tt,s.window="periodic"))
  • Fourier Transform
n<-500
x <- seq(from=0,to=10,length=500)
y <- sin(x) + 3*sin(runif(1)*x) + rnorm(n,0,0.5)
# Fourier and inverse-Fourier
# Identical
fft.out <- fft(y)
z <- fft(fft.out,inverse=TRUE)
plot(x,y)
points(x,Re(z)/length(z),pch=20,col=2,cex=0.6)
  • Smoothing with Fourier: Take out "BIG COMPONENT" only
# BIGNESS is Mod of complex values
hist(Mod(fft.out))
fft.out2 <- fft.out
fft.out2[which(Mod(fft.out) < 400)] <- 0
z2 <- fft(fft.out2,inverse=TRUE)
plot(x,y)
points(x,Re(z2)/length(z2),pch=20,col=2,cex=0.6)