How can I match an element of a vector after some matrix calculations in R?-CodePudding

I have a stochastic matrix (transition first order probabilities) in R:

P = matrix(c(1,0,0,0,0,
             0.4,0,0.6,0,0,
             0,0.4,0,0.6,0,
             0,0,0.4,0,0.6,
             0,0,0,0,1),5,5,byrow = TRUE) ;P

I have created a function in order to calculate the probability of first hitting a state before another ,given that the chain has started in a specific state.

For example the matrix is 5x5 so I ask what is the probability of first hitting state 5 before hit state 1, given that the chain starts from state 3. (here 5 and 1 are absorbing states).

The problem that I face is that the Omega vector inside the function contains 3 elements of states 2,3,4. But I want the element corresponding to state 3 not the third element as I have done. The trick with start -1 as I have done is only for this specific problem.

In general I want this function to take any arbitrary matrix, do the matrix calculations as they are written in my function and at the end to match somehow the started state with its corresponding element in the omega vector .


first_hit_before = function(mat,hit,before,start){
  Q = P[-c(hit,before), -c(hit,before)];Q
  R = P[-c(hit,before),c(hit,before)];R
  I = diag(dim(Q)[1]) ;I
  Omega = (solve(I-Q)%*%R)[,1];Omega
  Time = Omega[start-1];Time
  return(Time)
}
first_hit_before(P,5,1,3)

other way is by doing this :

first_hit_before <- function(mat,hit,before,start)
{
  Q = P[-c(hit,before), -c(hit,before)]
  R = P[-c(hit,before),c(hit,before)]
  I = diag(dim(Q)[1]) ;I
  Omega = (solve(I-Q)%*%R)[,1]
  names(Omega) = c(2, 3, 4)
  Time = Omega[as.character(start)]
  return(unname(Time))
}
first_hit_before(P,5,1,3)

But it doesn't help when the matrix is huge and the arguments are different resulted to bigger vector omega.

How this can be done in R ?

Any help ?

CodePudding user response：

A few changes:

Identify the columns at the beginning of the function; then remove Ids associated with removed columns:

id <- 1:ncol( mat )

id <- id[ !id %in% c( hit, before )]

Use mat instead of P in these two lines so the function performs calculations on the given matrix:

change Q = P[-c(hit,before), -c(hit,before)]

to Q = mat[-c(hit,before), -c(hit,before)

and

change R = P[-c(hit,before),c(hit,before)]

to R = mat[-c( hit,before),c( hit,before)]

Use something other than I as the name of the diagonal matrix, since I is a base R function to use a parameter AsIs, inhibiting some conversions.

first_hit_before <- function( mat, hit, before, start ){
  id <- 1:ncol( mat )                                # Identify columns
  id <- id[ !id %in% c( hit, before )]               # Remove names of removed columns
  Q  <- mat[ -c( hit, before ),  -c( hit, before ) ] # mat instead of P
  R  <- mat[ -c( hit, before ),   c( hit, before ) ] # mat instead of P
  D  <- diag( dim( Q )[ 1 ] )                        # D instead of I
  Omega <- ( solve( D - Q ) %*% R )[ , 1 ]           # D instead of I
  Omega[ id == start ]                               # Select start state
}