I am trying to initialize a model. I tried to create a simple example below that illustrates the issue. Of course the example is ridiculous and the problem would be solved, if I just specified the proper range for the parameter. My real issue really arises from underflow when I try initial values, for whcih the log-posterior is is such that I end up taking the log of a posterior that is very, very close to 0.
In the example below, the third example for initialization tries a reasonable guess near a good guess and perturbs it. In particular, I am not 100% sure that the guess would lead to a valid log-likelihood, but I would guess that I perturb in that area, I should hopefully (eventually) be fine. But I may get it rejected, but if I tried often enough it should work. Is there a way to try random inits repeatedly?
Another really weird thing is that for my real model (way too complicated to show), I actually get “Initialization between (-2, 2) failed after 100 attempts.”, even though I provided a function for random initialization. I did not manage to reproduce this in the simple example below and am totally puzzled why rstan tries the initialization betwen -2 and 2. I did provide a function that gave an initial value for all parameters defined in the parameters block, should it then not do that, at all??
example_code ← "
functions {
real shiftedbeta_lpdf(real y, real alpha, real beta){
return( beta_lpdf(y-3.0 | alpha, beta) );
}
}
data {
int y;
int n;
}
parameters {
real piplusthree;
}
model {
piplusthree ~ shiftedbeta(0.5, 0.5);
y ~ binomial(n, piplusthree-3);
}
"
# Try 1: This does not work
stanfit4 ← stan(model_code=example_code, data=list(y=5, n=10))
# Try 2: This of course works fine
stanfit5 ← stan(model_code=example_code, data=list(y=5, n=10), init=list(list(piplusthree=3.5),list(piplusthree=3.5),list(piplusthree=3.5),list(piplusthree=3.5)))
# Try 3: As far as I can tell this just generates a single random proposal that gets rejected
initf1 ← function(){
return( list(piplusthree=rnorm(n=1, mean=3.5,sd=1)) )
}
stanfit6 ← stan(model_code=example_code, data=list(y=5, n=10), init=initf1, seed=123)