MCMC sampling does not work when execute

Modeling
61

Thank you so much. By

but you have to do that for each observation.

you mean for each observed y[i] value? But what I want to get is the posterior distribution of target variable though. So basically, I have to put a for-loop around this line?

stanfit <- sampling(stanmodel1, data = list(ncol = ncol,nrow = nrow,
                                        yH = yH, x=x,
                                        A = A, sigma_x = sigma_x, sigma_y = sigma_y, epsilon = epsilon)
                ,iter=iterations, warmup = 1000, chains = 4, cores = 2, refresh = 250);
62

Yes, you have to loop over both observations and possible values for the discrete integer.

1 Like
63

But I just realize sampling already indicates the number of iterations, which essentially is a for loop? So do we need to loop it again?

64

Nested for loops: The outer one over all observations and the inner one over all possible values of the discrete unknown.

1 Like
65

So do you mean this:

for (i in 0: ncol) {
      x <- rpois(1)
      y <- rnorm(1, mean = x, sd = 1)
      vector[big_integer] log_weights;
       vector[big_integer] log_dens;
       for (i in 0:big_integer) {
           log_weights[i] = poisson_log_lpmf(i | log_lambda);
            log_dens[i] = normal_log_lpdf(y | i, 1);
       }
      target += log_sum_exp(log_weights + log_dens);
}

This part should be entirely in the model part, correct? And the parameters{} block should contain only vect[ncol] target? Btw, sampling() does not work with target variable, because of certain default in name.

66

Sorry, but this is what I got so far:

Stan Model

 data {
  int<lower=1> ncol;
  int<lower=1> nrow;
  vector[ncol] yH;
  int x[nrow];
  matrix[nrow, ncol] A;
  vector[ncol] log_sigmay;
  vector[nrow] epsilon;
  vector[ncol] log_lambda;
  vector[nrow] k;
}

parameters {
   vector[ncol] targ;
}

model { 
   targ += log_sum_exp(log_weights + log_dens);
}

Fitting MCMC in R

library(rstan)
library(MASS)
library(truncnorm)
library(bayesplot)
library(bridgesampling)
library(ggplot2)
library(MCMCpack)
library(bayesplot)
library(clusterGeneration)
require(MASS)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

nrow = 16;
ncol = 24;
iterations = 1500;
A ā† matrix(c(1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_A^1 - Done
0,0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_A^2 - Done
0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_B^1 - Done
1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_B^2 - Done
0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_C^1 - Done
0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0, #x_C^2 - Done
0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0, #x_D^1 - Done
0,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, #x_D^2 - Done
0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0,0,0,0, #x_E^1 - Done
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,1,0, #x_E^2 - Done
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0, #x_F^1 - Done
0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,1,0,0, #x_F^2 - Done
0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,0,0,0, #x_G^2 - Done
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1, #x_G^1 - Done
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,0,0,0, #x_H^1 - Done
0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,1), #x_H^2 - Done
nrow, ncol, byrow= TRUE);

for (i in 1:1000) # first for-loop
{

yH ā† as.integer(rpois(ncol, 100));
x ā† as.integer(A%%yH + rtruncnorm(nrow, a = 0, b = Inf, mean = 10, sd = 5))
sigmay ā† array(12, ncol);
epsilon ā† array(rtruncnorm(nrow, a = 3, b = 8, mean = 6, sd = 3), nrow);
lambda ā† runif(ncol, 20, 100)
log_lambda ā† log(lambda)
k ā† A%%yH
log_sigmay ā† log(sigmay)

  stan_data = list(nrow = nrow, ncol = ncol, yH = yH, x = x, k = k, A = A, log_sigmay = log_sigmay, log_lambda = log_lambda, epsilon = epsilon);

  stanmodel1 <- stan_model(file = "Poisson_Normal.stan", model_name = "stanmodel1");

}

stanfit ā† sampling(stanmodel1, data = stan_data, iter=iterations, warmup = 1000, chains = 4, cores = 2, refresh = 250);

I keep getting this error (the one you wrote also informed this error). Could you please help me fix this model so that it could run? Although I am not even sure if it would do what i wantā€¦

targ = targ + log_sum_exp(log_weights + log_dens);
PARSER EXPECTED: <expression assignable to left-hand side>
Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model 'stanmodel1' due to the above error.

Another error tht I got, after creating a transformed parameter block like following:

transformed parameters{
vector[ncol] log_weights;
vector[ncol] log_dens;
for (i in 0:ncol){
log_weights[i] = poisson_log_lpmf(i|log_lambda[i]);
log_dens[i] = normal_log_lpdf(targ[i]|(log(k)+epsilon)[i], log_sigmay[i]);
}
}

is on the normal_log_lpdf() line:

No matches for:

normal_log_lpdf(real, real, real)

Function normal_log_lpdf not found.

I am so sorry, but I really have no idea how to fix this issue.

67

@bgoodri: could you please help me resolve the issue above? I tried my best to work out the marginalization under your instruction, but I really could not get the code above running despite spending hours on fixing it;p

68

These are just basic errors with the Stan language, which you should be more familiar with before trying a mixture model. I accidentally copy and pasted an error with normal_log_lpdf, which does not exist. It is just normal_lpdf. The target keyword is what you have to increment to get your log-likelihood; it is not an unknown that you declare in the parameters block. The way you have written things, everything is known except the discrete integer so I donā€™t even know what you are trying to infer.

This is about right, but you have to define big_integer in the data block and pass it from R.

for (i in 0: nrow) {
      vector[big_integer] log_weights;
       vector[big_integer] log_dens;
       for (i in 0:big_integer) {
           log_weights[i] = poisson_log_lpmf(i | log_lambda);
            log_dens[i] = normal_lpdf(y | i, 1);
       }
      target += log_sum_exp(log_weights + log_dens);
}
1 Like
69

Thank you for your patience. My apology for being so dumb that I have to ask these questions. I did change it to normal_lpdf, but it also did not work. Actually, what I want to make inference on, is the posterior distribution of target. Because my problem is given as:
Given:
yT|yH ~ Poisson(lambda_i) (prior distribution, with the Poisson rate ~ Uniform(0, 100))
x|yT ~ N(A*yT, sigma) (sigma ~ Uniform(0,50), and A = matrix defined above)

Use Bayesian (and MCMC to deal with the intractable normalizing constant) to make inference on the posterior distribution of $yT|x$.

I was pondering with the way you define normal_lpdf(y | i, 1). First of all, y is a vector in my problem, so I am not sure if normal_lpdf(y | i, 1) is appropriate. What I was thinking is normal_lpdf(y| A%*%yT, sigma). But then Stan gave an error message on the expression A%*%yT, which is weird because their dimensions match (A*yT only gives element-wise multiplication, so it is not what we need).

70

@bgoodri: could you please be patient and help me with this really weird error (because this is what is shown in the manual)

targ +=  log_sum_exp(log_weights + log_dens);
      ^
PARSER EXPECTED: <expression>

Itā€™s really weird that when I tried fitting the stan_model() function, it announces the error with the common expressions like += or %*%. Very unsure what is going on with stan.

71

The keyword is target, which can be incremented with the += operator. On the other hand, you were trying to declare targ in the parameters block, in which case in cannot be incremented.

72

The reason I declared targ in the parameters block is because I want to make inference on that variable. I thought thatā€™s what marginalization implies? Because we take the logarithm on both sides of p(y|yT,x) ~ p(yT|yH)p(x|yT, A). So, the targ variable should stand for log(p(y|yT,x)), isnā€™t it?

73

We donā€™t have unlimited time to help out with basic errors and teach people Stan one error message at a time. As @bgoodri suggested, you might want to get some practice writing simpler Stan programs before trying something complex.

+= isnā€™t defined in R and %*% isnā€™t defined in Stan. We have a manual with an index defining all the functions and operators defined in Stan and only those may be used in a Stan program. You canā€™t use arbitrary R.

1 Like
74

True. I am trying very hard to make it work, but this modeling problem sucks so much. Could you give a hand, as there is some serious problem with the above marginalization? We all know we have to use marginalization, but the specific details have to be correct, and so far, they arenā€™t:(

75

You might consider asking somebody locally with more math background how to go about figuring this out.

1 Like
77

I donā€™t understand why the math is hard. I formulated the paperā€™s version above, but then Mr. @bgoodriā€™s code seems not to work exactly the same as the paperā€™s formulation. The only difficult thing is to code it up with the marginalization trick. The prior and likelihood are already given (Poisson with mean(lambda) where lambda ~ U(0,100), while I got the likelihood part already, which follows normal with mean = A*yT and constant variance (yT = the variable we try to find the posterior conditional distribution over). Why is it so hard in Stan for the prior part??

78

This is math, it has nothing to do with Stan.

79

So this is my math: p(y|yT, x) ~ p(x|yT)p(yT|y, A) ( ~ defines proportional up to a constant). Now, by marginalization for the poisson distribution p(d|dH, A), we have: p(y|yT, A) = sum_{x in all possible values of x} p(x,y|yT, A) = sum_{x in all possible values of x} p(x|yT, A)p(y|x, yT, A). But y|x,yT,A is the same as y|x,yT, A, so marginalization does not work in this case;p

80

When you say ā€œmarginalization does not work in this caseā€, what do you mean? Marginalization is just an operation on distributions allowing you to define p(u) in terms of p(u, v). Can you write down the joint density and how you comptue the marginal? I couldnā€™t follow your derivation of p(y | yT, x) ā€” I donā€™t see how you derive

p(x, yT | y, A) = p(x|yT) p(yT|y, A)

to be proportional to p(y | yT, x).