差分

このページの2つのバージョン間の差分を表示します。

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
r:markovchain [2018/12/12 07:48] – watalur:markovchain [2018/12/17 11:47] (現在) – watalu
行 7: 行 7:
   * [[https://cran.r-project.org/web/packages/markovchain/vignettes/markovchainCrashIntro.pdf|Crash introduction ]]   * [[https://cran.r-project.org/web/packages/markovchain/vignettes/markovchainCrashIntro.pdf|Crash introduction ]]
   * [[https://cran.r-project.org/web/packages/markovchain/vignettes/an_introduction_to_markovchain_package.pdf|An introduction]]   * [[https://cran.r-project.org/web/packages/markovchain/vignettes/an_introduction_to_markovchain_package.pdf|An introduction]]
 +
 +このパッケージはR 3.5.0以上を必要とする。各自でダウンロードしてインストールする必要があるかもしれない。
 +
 +=== マルコフ解析 ===
 +
 +マルコフ過程に基づいて現象の性質を明らかにしたり、未来の時点のおける状態の確率分布を推定することなどを、マルコフ解析という。
 +
 +=== マルコフ解析に必要なもの ===
 +
 +マルコフ解析に必要なものは、現在の状態もしくは初期の状態と、対象の状態遷移行列である。
 +
 +  * 状態
 +  * 状態遷移行列
 +
 +天気の移り変わりを例に説明を進める。
 +
 +<code>
 +weatherStates <- c("sunny", "cloudy", "rain")
 +</code>
 +
 +<code>
 +byRow <- TRUE
 +</code>
 +
 +<code>
 + weatherMatrix <- matrix(data = c(0.70, 0.2, 0.1,
 + 0.3, 0.4, 0.3,
 + 0.2, 0.45, 0.35), byrow = byRow, nrow = 3,
 + dimnames = list(weatherStates, weatherStates))
 +</code>
 +
 +<code>
 +mcWeather <- new("markovchain", states = weatherStates, byrow = byRow,
 + transitionMatrix = weatherMatrix, name = "Weather")
 +</code>
 +
 +
 +=== markovchainのメソッド ===
 +
 +|Method|Purpose|
 +|*|Direct multiplication for transition matrices.|
 +|[|Direct access to the elements of the transition matrix.|
 +|==|Equality operator between two transition matrices.|
 +|as|Operator to convert markovchain objects into data.frame and table object.|
 +|dim|Dimension of the transition matrix.|
 +|names|Equal to states.|
 +|names<-|Change the states name.|
 +|name|Get the name of markovchain object.|
 +|name<-|Change the name of markovchain object.|
 +|plot|plot method for markovchain objects.|
 +|print|print method for markovchain objects.|
 +|show|show method for markovchain objects.|
 +|sort|sort method for markovchain objects.|
 +|states|Name of the transition states.|
 +|t|Transposition operator (which switches byrow slot value and modifies the transition matrix coherently).|
 +
 +=== markovchainの使い方の基本 ===
 +
 +上で定義したマルコフ連鎖を確認しておく。
 +
 +<code>
 +weatherStates
 +weatherMatrix
 +mcWeather
 +</code>
 +
 +さて、初期状態は曇りcloudy。
 +
 +<code>
 +initialState = c(0, 1, 0)
 +</code>
 +
 +二日後の天気の確率分布を求めてみる。
 +
 +<code>
 +after2Days <- initialState * (mcWeather * mcWeather)
 +</code>
 +
 +求めた確率分布を確認する。
 +
 +<code>
 +after2Days
 +</code>
 +
 +７日後の天気の条件付き確率分布を推定する。
 +
 +<code>
 +after7Days <- initialState * (mcWeather ^ 7)
 +</code>
 +
 +求めた確率分布を確認する。
 +
 +<code>
 +after7Days
 +</code>
 +
 +こんなに長い桁は必要ないので、四捨五入して小数点以下第３位まで表示してみる。
 +
 +<code>
 +round(after7Days, 3)
 +</code>
 +
 +状態を列ベクトルで表すなら、すべてを転地すればいい。
 +
 +<code>
 +initialState <- c(0, 1, 0)
 +after2Days <- (t(mcWeather) * t(mcWeather)) * initialState
 +after7Days <- (t(mcWeather) ^ 7) * initialState
 +after2Days
 +</code>
 +
 +<code>
 +fvals<-function(mchain,initialstate,n) {
 +out<-data.frame()
 +names(initialstate)<-names(mchain)
 +for (i in 0:n)
 +{
 +iteration<-initialstate*mchain^(i)
 +out<-rbind(out,iteration)
 +}
 +out<-cbind(out, i=seq(0,n))
 +out<-out[,c(4,1:3)]
 +return(out)
 +}
 +fvals(mchain=mcWeather,initialstate=c(90,5,5),n=4)
 +</code>
 +
 +関数statesはマルコフ連鎖のオブジェクトmcWeatherの中の状態の一覧を表示させる。
 +<code>
 +states(mcWeather)
 +</code>
 +
 +関数namesもstatesと同じ動作をする。
 +<code>
 +names(mcWeather)
 +</code>
 +
 +関数dimは状態数を返す。
 +<code>
 +dim(mcWeather)
 +</code>
 +
 +関数nameはマルコフ連鎖の名前を返す。これはnewで作成するときにつけている。
 +<code>
 +name(mcWeather)
 +</code>
 +
 +
 +<code>
 +markovchain:::sort(mcWeather)
 +</code>
 +
 +
 +<code>
 +transitionProbability(mcWeather, "cloudy", "rain")
 +</code>
 +
 +
 +<code>
 +mcWeather[2,3]
 +</code>
 +
 +
 +<code>
 +print(mcWeather)
 +</code>
 +
 +
 +<code>
 +show(mcWeather)
 +</code>
 +
 +<code>
 +plot(mcWeather, package="diagram",box.size = 0.04)
 +</code>
 +
 +<code>
 +mcDf <- as(mcWeather, "data.frame")
 +mcNew <- as(mcDf, "markovchain")
 +mcDf
 +</code>
 +
 +
 +<code>
 +mcIgraph <- as(mcWeather, "igraph")
 +</code>
 +
 +<code>
 +require(msm)
 +Q <- rbind ( c(0, 0.25, 0, 0.25),
 +c(0.166, 0, 0.166, 0.166),
 +c(0, 0.25, 0, 0.25),
 +c(0, 0, 0, 0) )
 +cavmsm <- msm(state ~ years, subject = PTNUM, data = cav, qmatrix = Q, death = 4)
 +msmMc <- as(cavmsm, "markovchain")
 +msmMc
 +</code>
 +
 +<code>
 +myMatr<-matrix(c(.1,.8,.1,.2,.6,.2,.3,.4,.3), byrow=TRUE, ncol=3)
 +myMc<-as(myMatr, "markovchain")
 +myMc
 +</code>
 +
 +=== マルコフ連鎖の確率推論 ===
 +
 +|Method|Returns|
 +|absorbingStates|the absorbing states of the transition matrix, if any.|
 +|steadyStates|the vector(s) of steady state(s) in matrix form.|
 +|communicatingClasses|list of communicating classes. sj , given actual state si.|
 +|canonicForm|the transition matrix into canonic form.|
 +|is.accessible|verification if a state j is reachable from state i.|
 +|is.irreducible|verification whether a DTMC is irreducible.|
 +|period|the period of an irreducible DTMC.|
 +|recurrentClasses|list of recurrent classes.|
 +|steadyStates|the vector(s) of steady state(s) in matrix form.|
 +|summary|DTMC summary.|
 +|transientStates|the transient states of the transition matrix, if any.|
 +
 +<code>
 +conditionalDistribution(mcWeather, "sunny")
 +</code>
 +
 +<code>
 +steadyStates(mcWeather)
 +</code>
 +
 +<code>
 +absorbingStates(mcWeather)
 +</code>
 +
 +<code>
 +summary(mcWeather)
 +</code>
 +
 +<code>
 +transientStates(mcWeather)
 +</code>
 +
 +<code>
 +is.accessible(object = mcWeather, from = "sunny", to = "cloudy")
 +</code>
 +
 +<code>
 +period(mcWeather)
 +</code>
 +
 +<code>
 +firstPassage(object = mcWeather, state = "sunny",n = 10)
 +</code>
 +
 +=== マルコフ連鎖のための統計的推測 ===
 +
 +|Function|Purpose|
 +|markovchainFit|Function to return fitted Markov chain for a given sequence.|
 +|predict|Method to calculate predictions from markovchain or markovchainList objects.|
 +|rmarkovchain|Function to sample from markovchain or markovchainList objects.|
 +
 +状態遷移のデータがあるとする。
 +
 +Alofi島の１日あたりの降雨のデータがある。
 +
 +<code>
 +data(rain)
 +</code>
 +
 +このデータは、２つのフィールドV1とrainを持つ。使うのは rain の方。
 +一日ごとの推移を集計してみる。
 +
 +<code>
 +table(data.frame(before=rainrain[−1],after=rainrain[-length(rain$rain)]))
 +</code>
 +
 +出力は次の通り。
 +
 +<code>
 +      after
 +before   0 1-5  6+
 +   0   362 136  50
 +   1-5 126  90  79
 +   6+   60  68 124
 +</code>
 +
 +この推移をグラフに表してみる。Gmiscパッケージの中の関数transitionPlotを使うために、準備をする。
 +
 +<code>
 +install.packages(c("rms","sandwich","stringr","Hmisc"))
 +reps = c("http://ftp.sunet.se/pub/lang/CRAN","http://cran.gforge.se")
 +install.packages("Gmisc", repos=reps, dependencies=TRUE, type="source")
 +</code>
 +
 +そして、使ってみる。
 +
 +<code>
 +require(Gmisc)
 +transitionPlot(table(data.frame(r1=rainrain[−1],r2=rainrain[-length(rain$rain)])),
 +  overlap_add_width=1.2,type_of_arrow="gradient",
 +  min_lwd = unit(2, "mm"), max_lwd = unit(15, "mm"),cex=1.5)
 +</code>
 +
 +{{:r:rain-transitionplot.png|}}
 +
 +<code>
 +library(dplyr)
 +library(stringr)
 +library(DiagrammeR)
 +library(networkD3)
 +</code>
 +
 +この準備で頑張ると、[[https://kazutan.github.io/kazutanR/koneta_state_trans_viz.html|こんな]]グラフも描けるらしい。
 +
 +この推移から、遷移行列を推定する。
 +
 +<code>
 +markovchainFit(rain$rain, method="mle", byrow=TRUE)
 +</code>
 +
 +次の推定値などが表示される。
 +
 +<code>
 +$estimate
 +            0       1-5        6+
 +0   0.6605839 0.2299270 0.1094891
 +1-5 0.4625850 0.3061224 0.2312925
 +6+  0.1976285 0.3122530 0.4901186
 +
 +
 +$standardError
 +             0        1-5         6+
 +0   0.03471952 0.02048353 0.01413498
 +1-5 0.03966634 0.03226814 0.02804834
 +6+  0.02794888 0.03513120 0.04401395
 +
 +$confidenceLevel
 +[1] 0.95
 +
 +$lowerEndpointMatrix
 +            0       1-5         6+
 +0   0.6034754 0.1962346 0.08623909
 +1-5 0.3973397 0.2530461 0.18515711
 +6+  0.1516566 0.2544673 0.41772208
 +
 +$upperEndpointMatrix
 +            0       1-5        6+
 +0   0.7176925 0.2636194 0.1327390
 +1-5 0.5278304 0.3591988 0.2774279
 +6+  0.2436003 0.3700387 0.5625151
 +
 +$logLikelihood
 +[1] -1040.419
 +</code>
 +
 +rain.MCに代入しておく。
 +
 +<code>
 +rain.mcfit <- markovchainFit(rain$rain, method="mle", byrow=TRUE)
 +</code>
 +
 +返して貰った内容を見てみる。
 +
 +<code>
 +str(rain.mcfit)
 +</code>
 +
 +rain.MCはS3オブジェクトのリストだが、rain.MC$estimateはS4オブジェクトとややこしい。
 +遷移行列を行列として取り出すには、次のようにリストの名前とその中のスロットを指定する。
 +
 +<code>
 +rain.mcfit$estimate@transitionMatrix
 +</code>
 +
 +これによって新たなマルコフ連鎖を定義するには、次のように書く。
 +
 +<code>
 +mcrain = new("markovchain", states = rain.mcfit$estimate@states,
 +      byrow = TRUE,
 +      transitionMatrix = rain.mcfit$estimate@transitionMatrix, name = "rain")
 +</code>
 +
 +関数markovchainFitは、markovchainオブジェクトをリストの要素として出力しているので、そちらを使うのも良い。
 +
 +<code>
 +mcrain = rain.mcfit$estimate
 +</code>
 +
 +
 +<code>
 +rain.mcfit <- markovchainFit(rain$rain, method="laplace", laplacian=0.5, byrow=TRUE)
 +</code>