issues using custom pdf; chains not moving

[Peter Smits]

Hi list,

I'm trying to implement the Reflected Beta distribution as described here http://paleobiol.geoscienceworld.org/content/42/2/240.abstract and I'm running into very strange issues when trying to fit the distribution to simulated data.

The pdf of the reflected beta is defined

f(x | lambda, theta) = (1 - lambda) / theta * (1 - x/theta)^(-lambda)  if lambda <= 0, 0 <= x <= 1 and 

f(x | lambda, theta) = (1 + lambda) / theta * (x / theta)^(lambda) if lambda > 0, 0 <= x <= 1

where x are individual occurrence times, lambda describes the "tilt" of the distribution and theta is maximum age of an occurrence. The function I've written for the reflected beta pdf is as follows:

/**

 * return probability of observation y given lambda and theta

 *

 * @param y value between 0 and theta of occurrence time

 * @param lambda shape of the sampling profile

 * @param theta duration of taxon

 * @return probability of occurrence time

 */

real reflected_beta(real y, real lambda, real theta) {

  real p;

   if(lambda > 0) {

    p <- ((1 + lambda) / theta) * pow((y / theta), lambda);

  } else if(lambda <= 0) {

    p <- ((1 - lambda) / theta) * pow((1 - y / theta), (lambda * -1));

  }

  return p;

}

I've followed https://groups.google.com/forum/#!msg/stan-users/4gv3fNCqSNk/4dwmJj1YwbsJ except I've written a function to perform the pdf calculation. I then use increment_log_probability on a vector of the log transformed probabilities; the increment and the log transform are both done in the model block.

In my simulation, all samples are generated with the same value of lambda but each "species" as a unique theta value generated from a shared Weibull distribution with constant shape and scale. My stan file reflects my attempts to estimate the key values (shape and scale) from the simulated data. When I try fitting the model the chains remain at their initial condition and don't move.  I've tried changing step size and adapt_delta to no avail.

I've attached the stan file and associated simulation data.

Are there any recommendations for how to proceed this?

Cheers,

Peter

[andre.p...@googlemail.com]

Peter,

I hope I'm not the only one who would like to see your input data 

to your provided output.

I checked your reflected_beta function and if I'm not hopeless wrong it

returns:

pow<-function(x,k) x**k

> theta

[1] 0.3589448

> lambda

[1] 0.3906438

> y

[1] 0.4934321

> ((1 + lambda) / theta) * pow((y / theta), lambda)

[1] 4.387076

This is not a valid probability. I struggled to download your paper to

see the reflected_beta definition of the paper. 

Andre

[Peter Smits]

Hi Andre,

Thank you for the reply. Sorry about the difficulties surrounding the paper. I've attached a copt of the relevant paper to this email.

I thought I provided the input data? The .data.R file is a dump from R formatted for use with the Stan model. I'll make sure to attach it this time.

I need to correct a mistake in my original definition: x values are constrained to be less than or equal to theta, not 1. This, however, does not explain your problem. This appears to happen when values of theta are approximately less than 1. I'll contact the authors and try and figure out what is going on.

Thanks for helping me notice that. Hopefully resolving that goes a long way to helping with this problem.

Cheers,

Peter

[Herra Huu]

I think the problem is here. The support of pdf is [0,theta], so theta has to be larger than the largest observed value of x. So, at least, those constraints should be added to the model (and the priors should be changed accordingly). Also, it's probably a good idea to write down the log-density directly and avoid those pow() calculations, that way it should be more numerically stable and maybe a bit more efficient as well.

[Peter Smits]

Herra and others,

Yes, you're absolutely right. I totally forgot to include the individual constraints and while it wasn't straight forward in Stan, I appear to have found a solution. I've amended my code to constrain the minimums for each individual theta following https://groups.google.com/forum/#!topic/stan-users/ibEPXA7apAU. It seems to have done the trick. New code is like the following:

parameters {

  ...

 vector<lower=1>[S] theta_raw;  // duration

  ...

}

transformed parameters {

  vector<lower=0>[S] theta;

  theta <- theta_raw .* d;

}

model {

  ...

  increment_log_prob(theta, alpha, sigma);

  ...

  

  for(n in 1:N) {

    y[n] ~ reflected_beta(lambda, theta[sp[n]]);

  }