EP Model sometime crashes immediately, other times continues

Modeling
#1

I am working with a 6 parameter model. Sometimes when I run it, it will crash immediately with

Exception: AdAttribEP_sg_model_namespace::write_array: eta1 is -1.66624, but must be greater than or equal to 0 (in ‘/home/ubuntu/adtech/ep-stan/src/models/AdAttribEP_sg.stan’ at line 213)

other times, it will start and dump out lots messages about rejecting the draw
Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Exception: validate transformed params: eta1 is -0.00213054, but must be greater than or equal to 0 (in ‘/home/ubuntu/adtech/ep-stan/src/models/AdAttribEP_sg.stan’ at line 213)

why the different behaviors? in particular, why stop running immediately?

As I add even more parameters, it consistently crashes right away. I need to figure this out to get it running. The only complexity really is the use of transformed parameters. I declared a number of parameters phi[dim], then declared transformed parameters like alpha1,2,3 and eta1, 2, 3, and set them equal to the phi’s. At the end, I also introduced arrays, like alpha[3], and set the array alpha parameters to alpha1, alpha2, etc. Is that detrimental to the model?

#2

I don’t know what’s happening with the different behaviors, but you can probably fix it by explicitly declaring it as a positive parameter:

real<lower=0> eta1;

This will work for arrays as well, so

real<lower=0> phi[dim];
#3

I do already have that when I declare eta1

real<lower=0> eta1;

Isn’t that why it complains about the constraint in the first place?

#4

It would probably be helpful to see the model if possible.

I’m not sure why you’re just copying the elements of phi into other variables; this isn’t really what the transformed parameters block is for. According to this:

https://groups.google.com/forum/#!topic/stan-users/RMlROwqHiug

your constraints in the transformed parameters block are checked every iteration. It sounds to me like you’re not declaring phi to be positive, so it samples negative values for phi and then tries to copy them into the positive eta parameters.

#5

I am already declaring the upper/lower bound on the transformed parameters, all the eta’s for sure. Since I declare eta1 = phi[4] and so on, I don’t suppose I need to redeclare phi[i] with these bounds?

Or maybe you are saying that’s what I should do? But I don’t see this done in sample code.

#6

Stan directly samples the parameters. The transformed parameters are computed from these. So you do need to declare the bounds on phi.

#7

Ok, then I will. But that’s wierd, because I am following this type of sample code
data {
int<lower=1> N;
int<lower=1> D;
matrix[N,D] X;
int<lower=0,upper=1> y[N];
vector[D+1] mu_phi;
matrix[D+1,D+1] Omega_phi;
}
parameters {
vector[D+1] phi;
real eta;
vector[D] etb;
}
transformed parameters {
real alpha;
real<lower=0> sigma_a;
vector[D] beta;
vector<lower=0>[D] sigma_b;
sigma_a <- exp(phi[1]);
alpha <- eta * sigma_a;
sigma_b <- exp(tail(phi, D));
beta <- etb .* sigma_b;
}
model {
phi ~ multi_normal_prec(mu_phi, Omega_phi);
eta ~ normal(0, 1);
etb ~ normal(0, 1);
y ~ bernoulli_logit(alpha + X * beta);
}

You can see the phi’s do not have bounds. Their derivative parameters, like sigma_a does.

#8

wait, how would I declare bounds on the phi array? each element has a different bound requirements. I tried

real<upper=0> phi[7];

but it says I am redeclaring parameter. Which is true. but how would I set bounds for the 7th element of phi[dim]?

#9

Every element of an array has to have the same bounds. It sounds like you should be declaring each eta, alpha, etc. as a different array in the parameters block, and getting rid of phi altogether. The example you are following here works because the elements of phi are unconstrained, then they are exponentiated so sigma_s are positive.

#10

oh boy. The use of phi array is a requirement to run EP. I don’t have the option of declaring separate arrays. I think this is a question for tuomassivula then.

tuomassivula, is there a way to get around this?

#11

ok, I got around it using the same trick. set to exp() where I need a positive parameter. I suppose I could add a negative sign if I ever need a negative parameter.

#12

That’s how our internal transforms work. Just make sure to include the log Jacobian if you’re trying to put a distribution on the transformed parameter.

#13

Bob brought up a point I missed earlier. I was under the impression you always declare a distribution on the basest variable, which would by my phi array, and not on the transformed variables. Is it in fact ok to declare a distribution on the declared variables?

#14

Suppose you define a parameter on the log scale and then take its exponent and give it a distribution:

parameters {
  real log_x;
}
transformed parameters {
  real x = exp(log_x);
}
model {
  x ~ normal(0, 2);
}

This won’t give you the same distribution for x as

parameters {
  real<lower = 0> x;
}
model {
  x ~ normal(0, 2);
}

Stan tries to tell you this by giving you a Jacobian warning. You need to add the Jacobian to give x a half-normal distribution:

parameters {
  real log_x;
}
transformed parameters {
  real x = exp(log_x);
}
model {
  x ~ normal(0, 2);
  target += x;  // log(abs(d/dx exp(x)))
}