2016-04-25 4 views
0

私は以下のタスクをうまく実装しましたが、小さなテストデータセットでしか動作しませんでした。私の実際のデータセットは、永遠に計算し続けます。s1:s400からsnの大きな確率でRの大きなマルコフ遷移行列を得る確率を計算する

私はRで400x400の遷移確率行列を持っています。ユーザーがマルコフウォークを変換すると、「変換」をヒットします。吸収状態は、すべてのユーザーにとって「Null」です。 「スタート」は私の最初の状態です。

2つのことを私は計算する必要があります。

  1. ヒット状態s_jをランダムウォークの「スタート」から始まる
  2. 397の他の状態のそれぞれから始まるランダムウォークの
  3. ヒット「変換」

第1は、Rで簡単です:私は首尾よく2を実装

v <- numeric(length = ncol(transitionMatrix1)) 
v[1] <- 1 
i <- 2 
R0 <- v%*%(transitionMatrix1 %^% 1) 
R <- R0 
repeat { 
    R1 <- v%*%(transitionMatrix1 %^% i) 
    R <- rbind(R, R1) 
    if (rowSums(R[nrow(R)-1,] - R1) == 0) { 
    #if (rowSums(R[nrow(R)-1,] - R1) < epsilon) { 
    break 
    } 
    else { 
    i <- i+1 
    } 
} 
visit1 <- colSums(R) 

が、私小さな行列でしか動作しませんでした。それは大きなもので永遠にかかる:

w <- i 
C1 <- matrix(nrow = w, ncol = ncol(transitionMatrix1)) 
for (i in 1:ncol(transitionMatrix1)) { 
    x <- numeric(length = ncol(transitionMatrix1)) 
    x[i] <- 1 
    for (j in 1:w) { 
    C1[j,i] <- x%*%(transitionMatrix1 %^% j)[,ncol(transitionMatrix1)-1] 
    } 
} 
convert1 <- colSums(C1) 

私はおそらくループを使用しないでください。残念ながら、私は上記の操作をベクトル化することに成功しませんでした。 @Forzaaは、あなたのリンクの助けを借りて、私はRに変換されたと、MATLAB(https://de.mathworks.com/matlabcentral/newsreader/view_thread/48984)で解決策を見つけることができた

+0

元のものと入れ替えたものを同時に2つの行列に保存して、それに応じて1つまたは別のものを使用するのに役立つもの –

+0

計算しようとしていることを完全には理解していませんが、繰り返し行列を乗算することなく、あなたが望む表現を得ることができると思います。クイック概要については、https://en.wikipedia.org/wiki/Absorbing_Markov_chainを、完全な分析についてはhttp://www.math.pku.edu.cn/teachers/yaoy/Fall2011/Kemeny-Snell1976.pdfを参照してください。基本的な行列を計算する必要がありますし、すべての種類のパフォーマンス測定値を取得することができます。 – Forzaa

答えて

0

ありがとう:

f.absorb <- function(P) { 
    # input: P an absorbing Markov transition matrix in 'from rows to column' form. Rows are 
    # probability vectors and must sum to 1.0. At least one main diagonal 
    # element=1 output: 
    # tau: Time before absorption from a transient state. 
    # B: P(ending up in absorbing state from a transient state) 
    # tau2: variance matrix for tau; 
    # N: fundamental matrix 
    # N2: covariance matrix for N; 
    # CP: transition matrix in canonical form. 
    # Q: Submatrix of transient states, 
    # R: from transient to absorbing state submatrix, 
    # H: hij the probability process will ever go to transient state j starting 
    # in transient state i (not counting the initial state (K & S, p. 61)) 
    # MR: expected number of times that a Markov process remains in a transient 
    # state once entered (including entering step) 
    # VR: Variance for MR (K & S, Theorem 3.5.6, p. 
    # 61) written by E. Gallagher, Environmental Sciences Program UMASS/Boston 
    # Internet: [email protected] written 9/26/93, revised: 10/26/94 
    # refs: Kemeny, J. G. and J. L. Snell. 1976. Finite Markov chains. 
    # Springer-Verlag, New York, New York, U.S.A. Roberts, F. 1976. Discrete 
    # Mathematical Models with applications to social, biological, and 
    # environmental problems. Prentice-Hall 
    # R adaption from MatLab Code found at 
    # https://de.mathworks.com/matlabcentral/newsreader/view_thread/48984 
    pdim = nrow(P) 
    if (rowSums ((colSums(t(P))-matrix(1,1,pdim)) * (colSums(t(P))-matrix(1,1,pdim))) > 0.1) { 
    stop('Rows of P must sum to 1.0') 
    } 
    dP=diag(P); 
    ri=which(dP==1.0); 
    if (is.null(ri)) { 
    stop('No absorbing states (reenter P or use ergodic.m)') 
    } 
    rdim=length(ri); 
    qi=which(dP!=1.0); 
    qdim=length(qi); 
    I=diag(qdim) 
    Q=P[qi,qi]; 
    N=solve(I-Q); # the fundamental matrix 
    tau=data.matrix(colSums(t(N))); 
    CP=P[c(t(ri),t(qi)),c(t(ri),t(qi))]; # the canonical form of P 
    R=CP[(rdim+1):pdim,1:rdim]; 
    B=N%*%R; 
    Ndg=diag(diag(N)); 
    N2=N%*%(2*Ndg-I)-N*N 
    tau2=(2*N-I)%*%tau-tau*tau; 
    H=solve(Ndg,N-I); 
    dQ=diag(Q); 
    oneq=matrix(1,qdim,1); 
    MR=oneq/(oneq-dQ); 
    VR=(dQ/(oneq-dQ)) * (dQ/(oneq-dQ)); 
    list(tau=tau, N=N, B=B, CP=CP, Q=Q, R=R, N2=N2, H=H, MR=MR, VR=VR) 
} 

私が探していた値を見つけることができます行列の中でH

関連する問題