calculation of expected time on random walk

model.stan

functions {
  real expected_time(int I, int J, int[,] can_goal, int[] di, int[] dj, vector p) {
    matrix[I*J,I*J] A;
    vector[I*J] b;
    vector[I*J] res;

    for (k1 in 1:(I*J)) {
      b[k1] = 0;
      for (k2 in 1:(I*J)) A[k1,k2] = 0;
    }

    for (i in 1:I) {
      for (j in 1:J) {
        if (i == I && j == J || can_goal[i,j] == 0) {
          A[(i-1)*J+j, (i-1)*J+j] = 1;
        } else {
          real move;
          move = 0;
          for (k in 1:4) {
            int i_next;
            int j_next;
            i_next = i + di[k];
            j_next = j + dj[k];
            if (1 <= i_next && i_next <= I && 1 <= j_next && j_next <= J && can_goal[i_next, j_next] == 1) {
              A[(i-1)*J+j, (i_next-1)*J+j_next] = -p[k];
              move = move + p[k];
            }
          }
          A[(i-1)*J+j, (i-1)*J+j] = move;
          b[(i-1)*J+j] = move;
        }
      }
    }
    res = A \ b;
    return res[1];
  }
}

data {
  int I;
  int J;
  int Can_Goal[I,J];
  int N;
  vector[N] Y;
}

transformed data {
  int Di[4];
  int Dj[4];
  Di[1] = -1; Di[2] = 1; Di[3] = 0; Di[4] = 0;
  Dj[1] = 0; Dj[2] = 0; Dj[3] = -1; Dj[4] = 1;
}

parameters {
  simplex[4] p;
  real<lower=0> sigma;
}

transformed parameters {
  real mu;
  mu = expected_time(I, J, Can_Goal, Di, Dj, p);
}

model {
  Y ~ normal(mu, sigma);
}

run-dfs.R

I <- 3
J <- 10
Grid <- matrix(c(c(T, F, T, T, T, F, T, T, T, F), 
                 c(T, F, T, F, T, F, T, F, T, F), 
                 c(T, T, T, F, T, T, T, F, T, T)), I, J, byrow=TRUE)
can_goal <- matrix(FALSE, I, J)
di <- c(-1, 1, 0, 0)
dj <- c(0, 0, -1, 1)

dfs <- function(i, j) {
  can_goal[i, j] <<- TRUE
  for (k in 1:4) {
    i_next <- i + di[k]
    j_next <- j + dj[k]
    if (1 <= i_next && i_next <= I && 1 <= j_next && j_next <= J && !can_goal[i_next, j_next] && Grid[i_next, j_next]) {
      dfs(i_next, j_next)
    }
  }
}

dfs(I, J)
print(can_goal)

run-model.R

library(rstan)

stanmodel <- stan_model(file='model.stan')

I <- 3
J <- 10
Can_Goal <- matrix(c(c(1, 0, 1, 1, 1, 0, 1, 1, 1, 0), 
                     c(1, 0, 1, 0, 1, 0, 1, 0, 1, 0), 
                     c(1, 1, 1, 0, 1, 1, 1, 0, 1, 1)), I, J, byrow=TRUE)
set.seed(1234)
N <- 10
Y <- round(rnorm(N, mean=80, sd=10), 1)

data <- list(I=I, J=J, Can_Goal=Can_Goal, N=N, Y=Y)
fit <- sampling(stanmodel, data=data, seed=123)