How to estimate a truncated, user-defined distribution?

Modeling
, , ,
1

Hello,
I want to estimate the parameter of a truncated distribution that is not prebuilt in Stan. In order to understand the process, I follow Annis, Miller, Palmeri (2016) in their approach to user-define the Exponential distribution. The prebuilt code works fine for me:

exp_prebuilt <- stan_model(model_code= '
  data {
  int N;
  vector[N] y;
  real L;
}
parameters {
  real<lower=0> lambda;
}
model {
  real alpha;
  real beta;
  alpha =1;
  beta =1;
  lambda ~ gamma(alpha, beta);
  for (i in 1:N) {
    y[i] ~ exponential(lambda) T[L,];
  }
}
')

I then tried to build the exponential distribution by myself, later I want to replace the Exponential with a user-defined function:

exp_user <-stan_model(model_code= '
functions{
  real newexp_log(vector x, real lam) {
      vector[num_elements(x)] prob;
      real lprob;
      for(i in 1:num_elements(x)){
        prob[i] =  lam*exp(-lam*x[i]);
      }
      lprob = sum(log(prob));
      return lprob;
  }
}
data{
  int<lower=0> N; 
  real y[N];
  real L;
}
parameters {
  real<lower=0> lambda;
}
model {
  real alpha;
  real beta;
  alpha =1;
  beta =1;
  lambda ~ gamma(alpha, beta);
  for (s in 1:N){
  y[s] ~ newexp(lambda) T[L,];
}
}
')

I tried to work with this post but am not sure if I implemented it correctly.

I receive two error messages for my code:

  1. No matches for: real ~ newexp(real) Available argument signatures for newexp: vector ~ newexp(real)

  2. Real return type required for probability function.

Any help in solving this problem, explanations, hints are appreciated!

2

Stan is quite strict about types. You’ve declared that the first argument to newexp_log is a vector but in the model you visit each element individually, meaning the type passed to newexp is a scalar. The error is trying to say that scalars and vectors are not compatible.
The solution is: change the function argument type to real and remove the loop from the fuction

real newexp_log(real x, real lam) {
    real prob;
    prob =  lam*exp(-lam*x);
    return log(prob);
}

The next problem is that truncation T[L,] requires a cumulative probability distribution function in addition to the probability density function.

real newexp_lcdf(real x, real lam) {
    real prob;
    prob = exp(-lam*x);
    return log(prob);
}

Stan also wants the complementary cumulative distribution function, CCDF.

real newexp_lccdf(real x, real lam) {
    real prob;
    prob = exp(-lam*x);
    return log(1-prob);
}

With these the model compiles.

3

Niko,
thank you so much for taking the time to help me!
I really appreciate your support and have a great day!