User defined pdfs (WALD and student bivariate) and failure to initialize

Modeling
, ,
1

I’m trying to fit a hierarichal model with a continuous latent variable (distest) coming from a WALD dispersal kernel and a higher order variable (dist2) coming from a student bivariate distribution. Both of these distributions are custom so I define the whole model as follows.

functions{
real wald_lpdf(real x, real mu, real lambda){ //defining wald pdf
	real prob;
	real lprob;
	prob = ((lambda/(2.00*pi()*(x^3.00)))^0.5)*exp(-(lambda*((x-mu)^2.00))/(2.00*(mu^2.00)*x));
	lprob  = log(prob);
	return lprob;
}

real student2dt_lpdf(real r, real p, real u){ //defining 2dt pdf
	real prob;
	real lprob;
	real area; 
	real kernel; 
	area = 2.00*pi()*r;
	kernel = (p/(pi()*u*(1.00+(r^2.00)/u)^(p+1.00)));
	prob = area * kernel;
	lprob = log(prob);
	return lprob;
}
}

data {
int <lower = 0> N;      //number of observations
real <lower = 0.0000000001> dist[N]; //trap distance
real <lower = 0> httrap[N]; //height on trap
real <lower = 0> mu_wald[N];
real <lower = 0> lambda[N];
real <lower = 0> sigma2[N];
}

parameters {
real<lower=0> p;	//additional predictor effects					
real<lower=0> u;          
real<lower=0.000000000001> distest[N]; //latent dispersal distance from WALD kernel
}

transformed parameters{
real <lower=0> dist2[N]; 

//calculate estimated distances
for(s in 1:N){
if(httrap[s] <= 0.2){
dist2[s] = dist[s]; //if height <= 0.2m, est distance is measured distance
}
else {
dist2[s] = dist[s] + distest[s]; //latent ground distance from wald kernel
}
}
}

model{
p ~ gamma(0.001,0.001);
u ~ gamma(0.001,0.001);

//likelihood holding out some data
for(n in 1:N){
target += wald_lpdf(distest[n]|mu_wald[n],lambda[n]);
target += student2dt_lpdf(dist2[n]|p,u);
}
}

I think I’m not correctly specifying the pdfs for both of these custom probability functions. The first is a wald dispersal kernel pdf (from which a distance, distest, is being drawn) and the second is a student 2dt pdf, from which a distance, dist2 is being drawn. Distest is latent, dist2 is the sum of distest and a measured value dist.

I’m getting the following errors.

“Chain 1: Rejecting initial value:
Chain 1: Log probability evaluates to log(0), i.e. negative infinity.
Chain 1: Stan can’t start sampling from this initial value.”

I’m including a subsection of my data (randomized).
MockData.RData (2.1 KB)

I’ve tried; specifying initial values (p = 1, u = 1), setting gamma priors on p and u to constrain them above 0 (since they cannot be 0 in the pdf), exp mu_wald and lamba and p and u (I’m a bit unclear on this since the lpdfs are specified via log probability). None of these has worked. I’ve also tried just separately running the wald and student bivariate distributions separately. The student bivariate will run, the wald still gives me the same error.

I edited this post to reflect recent updates I made.

2

Try adding “.0” behind each of the numbers in “1/2” so Stan understands you mean reals and not integers

1 Like
3

Fixed that problem. Thanks. Still have the other.

4

mu_wald and lambda both need to be > 0. See Inverse Gaussian distribution - Wikipedia. If I add a small value to those the program samples. Not sure you want this distribution if you have 0s for those values.

As for the program. It’s better if you just have everything on the log scale as well. That 0.000001 lb can just 0. The int x[N] syntax is going to be deprecated in favor of array[N] int x;. You can vectorize those distributions once you have nailed the program down if you need more speed too.

functions {
  real wald_lpdf(real x, real mu, real lambda) {
    //defining wald pdf
    return 0.5 * (log(lambda) - (log2() + log(pi()) + 3 * log(x)))
           - (lambda * ((x - mu) ^ 2.00)) / (2.00 * (mu ^ 2.00) * x);
  }
  
  real student2dt_lpdf(real r, real p, real u) {
    //defining 2dt pdf
    real logpi = log(pi());
    return log2() + logpi + log(r) + log(p)
           - (logpi + log(u) + (1.0 + p) * log1p((r ^ 2.00) / u));
  }
}
data {
  int<lower=0> N; //number of observations
  array[N] real<lower=0> dist; //trap distance
  array[N] real<lower=0> httrap; //height on trap
  array[N] real<lower=0> mu_wald;
  array[N] real<lower=0> lambda;
  array[N] real<lower=0> sigma2;
}
parameters {
  real<lower=0> p; //additional predictor effects					
  real<lower=0> u;
  array[N] real<lower=0> distest; //latent dispersal distance from WALD kernel
}
transformed parameters {
  array[N] real<lower=0> dist2;
  
  //calculate estimated distances
  for (s in 1 : N) {
    if (httrap[s] <= 0.2) {
      dist2[s] = dist[s]; //if height <= 0.2m, est distance is measured distance
    } else {
      dist2[s] = dist[s] + distest[s]; //latent ground distance from wald kernel
    }
  }
}
model {
  p ~ gamma(0.001, 0.001);
  u ~ gamma(0.001, 0.001);
  //likelihood holding out some data
  for (n in 1 : N) {
    target += wald_lpdf(distest[n] | mu_wald[n], lambda[n]);
    target += student2dt_lpdf(dist2[n] | p, u);
  }
}
1 Like