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)