Specifying a custom prior distribution

Modeling
1

I have already make the following stan model which is compiling and I obtain good results (I use cmdstanr).
In this model, the prior distributions for the logit of the components of \lambda are Normal(0,10).

// The input data: 
// 'n0' an array of integers giving the numbers of seronegative persons of the different ages in the sample, 
// 'n' an array if integers giving the numbers of persons of the different ages in the sample, 
// 'n0' and 'n'  of length 'N' the number of ages observed in the sample.
data {
  int<lower=0> N;
  array[N] int<lower=0> n0; //int<lower=0> n0[N];
  array[N] int<lower=0> n;//int<lower=0> n[N];
}

parameters {
    // specify logitlambda as a parameter
    vector[N] logitlambda; 
}

transformed parameters {
    // treat lambda as a transformed parameter
    vector[N] lambda;
    lambda[1:N] = inv_logit(logitlambda[1:N]);
    vector[N] theta;    
    vector[N] delta; 
    theta = cumulative_sum(lambda);
    delta = exp(-theta);
}

model {
  for (i in 1:N){
    logitlambda[i] ~ normal(0, 10);
  }
  n0 ~ binomial(n,delta);
}

But now I want to change the prior distributions for the components of my vector \lambda. Indeed, I want a priori a kind of Bernoulli-Logit distribution:
with two parameters s and max, lambda[i] should be equal to zero with probability s, and should follow a uniform(0,max) with probability 1-s.

Then I modified my model in the following way, but it is not compiling anymore:

data {
  int<lower=0> N;
  array[N] int<lower=0> n0; //int<lower=0> n0[N];
  array[N] int<lower=0> n;//int<lower=0> n[N];
  real<lower=0> max1;
  real<lower=0> s1;
}

parameters {
    // 'u' of length 'N', to compute the prior for lambda.
    vector[N] lambda;
}

transformed parameters {
    vector[N] theta;    
    vector[N] delta; 
    theta = cumulative_sum(lambda);
    delta = exp(-theta);
    vector[N] u; 
}

model {
  for (i in 1:N){
    u[i] ~ uniform(0,1);
    if (u[i] < s1)
      lambda[i] = 0;
    else 
      lambda[i] ~ uniform(0,max1);
    }
  n0 ~ binomial(n,delta);
}

The error is in the model block, and is the following:

Cannot assign to global variable ‘lambda’ declared in previous blocks.

Any help would be welcome. I understand that I can not assign lambda[i]=0 in the model block, but I do not know where to assign it. I tried different things without success. Maybe I should create a new function to define a custom distribution but I do not know how to do. I read this part of the manual concerning a zero-inflated Poisson (5.6 Zero-Inflated and Hurdle Models | Stan User’s Guide), but it does not concern a prior distribution, in this example it is the data which follows the new zero inflated poisson distribution, not a parameter (hence it concerns the likelihood, not the prior).

You can use the following data if you want to make some tests:

n <- rep(100,10)
n0 <- c(97,91,85,77,68,58,49,40,30,24)
data_list <- list(N = 10, n0 = n0 ,n=n, max1=0.5,s1=0.05)

# Ajusting model
fit <- mod$sample(
  data = data_list, 
  seed = 123, 
  chains = 4, 
  parallel_chains = 4,
  refresh = 0,
  iter_warmup=200,
  iter_sampling=2000,
  save_warmup = TRUE
)
2

A little update.
I have try the following: to specify directly the prior through target +=

data {
  int<lower=0> N;
  array[N] int<lower=0> n0; //int<lower=0> n0[N];
  array[N] int<lower=0> n;//int<lower=0> n[N];
  real<lower=0> max1;
  real<lower=0, upper=1> s1;
}

parameters {
    // 'u' of length 'N', to compute the prior for lambda.
    vector[N] lambda;
}

transformed parameters {
    vector[N] theta;    
    vector[N] delta; 
    theta = cumulative_sum(lambda);
    delta = exp(-theta);
    vector[N] u; 
}

model {
  for (i in 1:N){
    if (lambda[i] == 0)
      target += bernoulli_lpmf(1 | s1);
    else 
      target += bernoulli_lpmf(0 | s1) + uniform_lpdf(lambda[i] | 0, max1);
    }
  n0 ~ binomial(n,delta);
}

It is compiling, but not working, I have errors about the parameters computed:

Chain 4 Exception: binomial_lpmf: Probability parameter[1] is 4.3263, but must be in the interval [0, 1] (in ‘/tmp/RtmpwfVGWu/model-33a638e69c40.stan’, line 35, column 2 to column 25)
Chain 4 Initialization between (-2, 2) failed after 100 attempts.
Chain 4 Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.

3

Your implementation was mostly correct, the error is being caused by delta not being in the unit interval (I think you applied the wrong transformation).

I think it should already work if you fix that, but you should also provide parameter bounds for lambda or the sampler will be much less efficient (and I think the resulting prior will also not be the one you actually want to have).

Here’s the model with a couple more tweaks and it seems to work alright:

data {
  int<lower=0> N;
  array[N] int<lower=0> n0; //int<lower=0> n0[N];
  array[N] int<lower=0> n;//int<lower=0> n[N];
  real<lower=0> max1;
  real<lower=0, upper=1> s1;
}

parameters {
    array[N] real<lower=0, upper=max1> lambda;
}

transformed parameters {
    array[N] real theta;    
    theta = cumulative_sum(lambda);
}

model {
  for (i in 1:N){
    if (lambda[i] == 0)
      target += bernoulli_lpmf(1 | s1);
    else 
      target += bernoulli_lpmf(0 | s1) + uniform_lpdf(lambda[i] | 0, max1);
    }
  n0 ~ binomial_logit(n, theta);
}
4

In the latest implementation, lambda will never be equal to a literal zero, and so the code executes fine but the model is always on the second branch of the control flow.

What you are trying to achieve is in general hard. It’s not possible to include discrete parameters in Stan. If the length of lambda is quite short, it would be possible to marginalize over all of the discrete parameters in all of the 2^n possible combinations of zeros and nonzeros. Otherwise, this model cannot be fit in Stan. Software based on other MCMC samplers is nominally capable of fitting this model with discrete parameters, but it’s very difficult to verify that you’re recovering reliable inference, especially as n grows (if n is 20, then there are over a million combinatorial possibilities for zeros and non-zeros and it’s impossible for a sampler to visit them all even once).

There are other flavors of sparsity-inducing priors that you could consider using that are less extreme than the literal spike-and-slab that you’re discussing above, and that can be fit reliably using Stan. The best treatment of these that I’m aware of is Sparsity Blues

1 Like
5

Thanks a lot to both of you!
For Luna Fazio: your suggestion is working, and I did even keep the delta parameter without any problem, because I really need it in n0 ~ binomial(n,delta).
For Jacob Socolar: I keep thinking about your explanations. Indeed, the length of lambda can be large, more than 50. When I looked to the example of the zero-inflated poisson (5.6 Zero-Inflated and Hurdle Models | Stan User’s Guide) , I see a discretization, so I thought it was possible in stan to model a mixture of distributions with one of the distribution which is a dirac (like a spike-and-slab). But I realize that in the example of the zero-inflated Poisson, the if condition is on the observed data, and not on a parameter like in my case. So in the ezro-inflated Poisson example, the model can go on the first branch if y[n] == 0.
In my case, It is more difficult for me to understand why lambda will never be equal to zero. In fact I just wanted to specify a prior distribution, so for me I just wanted to modify the target accordingly.
Anyway, I will take your suggestions into account, and I will try to model a mixture for my prior which is not a spike-and-slab. Maybe a mixture of two Gaussians, with one centered on zero with a small variance.

6

Ignoring the HMC part of Stan and just thinking about a random walk metropolis-hastings MCMC routine: the sampler initializes lambda somewhere, and then makes a proposal to move lambda somewhere else. That proposal will never propose a literal zero, because the proposal is made with like 16 digits of precision. If the computation were truly continuous with infinite precision, note that the probability of a continuous parameter ever taking any particular exact value is zero.

In random walk MH, you could have another discrete parameter that is only allowed to be one or zero. That would get sampled as zero sometimes. But you can’t do a discrete parameter of this sort in Stan (and your code above doesn’t attempt to introduce such a parameter).