用 R 中的大马尔可夫转移矩阵计算从 s1:s400 到 sn 的概率需要永远

Posted 2023-03-29

【中文标题】用 R 中的大马尔可夫转移矩阵计算从 s1:s400 到 sn 的概率需要永远

【英文标题】：Calculating probability to get from s1:s400 to sn with large Markov transition matrix in R takes forever

【发布时间】：2016-08-19 01:06:43

【问题描述】：

我成功实现了以下任务，但我只让它与一个小型测试数据集一起工作。对于我的真实数据集，它只是一直在计算。

我在 R 中有一个 400x400 的转换概率矩阵。如果用户在马尔可夫步行中转换，则点击“转换”。对于所有用户，吸收状态为“Null”。 “开始”是我的开始状态。

我需要计算两件事：

在从“Start”开始的随机游走中命中状态 s_j 在从其他 397 个州开始的随机游走中点击“转化”

第一个在 R 中很容易：

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)

我成功实现了 2.，但我只让它与一个小矩阵一起工作。一个大的要花很长时间：

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)

我可能不应该使用循环。不幸的是，我没有成功矢量化上述操作。

【问题讨论】：

帮助我的是同时将两个矩阵保存在内存中 - 原始的一个和转置的一个，并相应地使用一个或另一个

虽然我不完全理解您要计算的内容，但我认为您无需重复乘以矩阵即可获得所需的表达式。有关快速概述，请参阅en.wikipedia.org/wiki/Absorbing_Markov_chain，有关完整分析，请参阅math.pku.edu.cn/teachers/yaoy/Fall2011/Kemeny-Snell1976.pdf。您需要计算基本矩阵，然后才能获得各种性能指标。

【参考方案1】：

感谢@Forzaa，借助您的链接，我能够在 MatLab (https://de.mathworks.com/matlabcentral/newsreader/view_thread/48984) 中找到解决方案，并将其转换为 R：

 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: [e-mail] 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中找到