差分

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

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
r:markovchain [2018/12/12 08: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以上を必要とする。各自でダウンロードしてインストールする必要があるかもしれない。
  
 === マルコフ解析 === === マルコフ解析 ===
行 41: 行 43:
 </code> </code>
  
-<code> 
-mcList <- new("markovchainList", markovchains = list(mcWeather, defaultMc), 
- name = "A list of Markov chains") 
-</code> 
  
 === markovchainのメソッド === === markovchainのメソッド ===
行 136: 行 134:
 </code> </code>
  
 +関数statesはマルコフ連鎖のオブジェクトmcWeatherの中の状態の一覧を表示させる。
 <code> <code>
 states(mcWeather) states(mcWeather)
 </code> </code>
  
 +関数namesもstatesと同じ動作をする。
 <code> <code>
 names(mcWeather) names(mcWeather)
 </code> </code>
  
 +関数dimは状態数を返す。
 <code> <code>
 dim(mcWeather) dim(mcWeather)
 </code> </code>
  
 +関数nameはマルコフ連鎖の名前を返す。これはnewで作成するときにつけている。
 <code> <code>
- name(mcWeather)+name(mcWeather)
 </code> </code>
  
行 265: 行 264:
 |predict|Method to calculate predictions from markovchain or markovchainList objects.| |predict|Method to calculate predictions from markovchain or markovchainList objects.|
 |rmarkovchain|Function to sample 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=rain$rain[-1], after=rain$rain[-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=rain$rain[-1],r2=rain$rain[-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>