RTisean-tseriesChaos.R

# 岡田昌史(著者代表)『Rパッケージガイドブック』東京図書
# http://www.tokyo-tosho.co.jp/books/ISBN978-4-489-02097-1.html
#
# Google http://www.google.com/search?q=ISBN978-4-489-02097-1
# Amazon http://www.amazon.co.jp/dp/448902097X
#
# オリジナルのままであればソースファイルの再配布は自由です。

# RTisean-tseriesChaos.R
# 20110425版

# まず，TISEAN をダウンロードしましょう。TISEAN の公式ページ
#（http://www.mpipks-dresden.mpg.de/~tisean/Tisean_3.0.1/index.html）の，
# 左側下方にある“Downloads, updates, announcements”をクリックして
# 開くページから，お使いのOS に応じたTISEAN をダウンロードします


library(RTisean)
library(tseriesChaos)
laser <- scan("http://www-psych.stanford.edu/~andreas/Time-Series/SantaFe/A.dat")


## Theiler 窓の推定
op <- par(no.readonly=TRUE)
par(mfrow=c(2,2))
# 時間遅れd=2,...,5, 埋め込み次元m=1,...,6,
# 時刻t=1,...,50 に対して時空間分離プロットを実行
laser.stp <- sapply(2:5, function(d) stplot(laser,
                    m=6, d=d, mdt=50))
par(op)


## アトラクタの再構成
laser.mu <- mutual(laser, plot=FALSE)

laser.mu.d <- diff(laser.mu)
d <- min(which(laser.mu.d > 0)) - 1
d

# 誤り最近傍法
# ( 埋め込み次元m は10 まで, 時間遅れd=2,Theiler 窓t=5, 許容誤差eps=5)
laser.fnn <- false.nearest(laser, m=10, d=d, t=5, eps=5)
m <- min(which(laser.fnn[1, ] < 0.01))
m


## Lyapunov 指数の推定
# 近傍半径
eps.all <- c(1, 2, 4, 8, 16)
# 埋め込み次元m(=4), 時間遅れd(=2), Theiler 窓5,
# 計算に使用する点の数500, 時刻1 ～ 100, 近傍半径eps=1, 2, 4, 8, 16
laser.lk <- sapply(eps.all,
                   function(eps) {
                                  lyap_k(laser, m=m, d=d, t=5, ref=500, s=100,
                                  eps=eps)
                   }
)
colnames(laser.lk) <- paste("eps=", eps.all, sep="")
matplot(laser.lk, type="l", xlab="t", ylab="S(t)")
legend(60, 1, paste("eps=", eps.all), lty=1:5, col=1:5,
       cex=0.8)

laser.lyap <- apply(laser.lk, 2,
                    function(x) lyap(x, start=30, end=60))
laser.lyap