Weighted log-probabilities in generalised gamma distribution

Question

This question is related to the problems I mentioned in this question. I am not sure if there is a good solution, but am hoping someone more experienced with this type of thing can help out.

I am characterising my observations as forming a (lower) truncated generalised gamma distribution. The (usual, non-truncated) generalised gamma distribution can be parameterised like this: $$ f(m|H_s,\alpha,\beta) = \frac{\beta \left(\frac{m}{H_s}\right)^\alpha \exp\left(-\left(\frac{m}{H_s}\right)^\beta\right)}{H_s \Gamma\left(\frac{\alpha+1}{\beta}\right)}, $$ where the denominator is just the normalisation to form a PDF. This is usually defined for $m>0$, $\beta>0$, $H_s>0$ and $\alpha>-1$ (so that the argument to the gamma function is positive). However, my data has $\alpha \sim -2$ over the range of its observation. Usually this would be impossible, since $\alpha<-1$ means the integral cannot converge towards low $m$ (I use the variable $m$ because my data are actually masses). However, I only have masses larger than some $m_{\rm min}$, so the distribution is truncated. In that case, the distribution is $$ f(m|H_s,\alpha,\beta) = \frac{\beta \left(\frac{m}{H_s}\right)^\alpha \exp\left(-\left(\frac{m}{H_s}\right)^\beta\right)}{H_s \Gamma\left(\frac{\alpha+1}{\beta},\left(\frac{m_{\rm min}}{H_s}\right)^\beta\right)}, $$ and $\alpha$ is now free to be $<-1$ (in my case it would be unphysical for it to be $<-2$).

My problem comes in because most libraries only implement the incomplete gamma for positive shape parameter, ie. $\alpha>-1$. So, though the distribution 'makes sense', many libraries cannot actually calculate it. There are some that can, but I want to be more flexible, and specifically, I want to use Stan (which doesn't support this at the moment).

One way to circumvent the problem (with some caveats), is to "weight" each data point by its mass, and then implicitly estimate $\alpha+1$. Since my data must have $\alpha>-2$, then $\alpha+1$ is always greater than -1, and this will work. The downside of course is that data with more noise is upweighted, but there are additional benefits which I won't go in to.

To make things explicit, here is the Stan code (hopefully this question is not actually stan-specific, but I think the code should be easy to understand anyway). The actual distribution is implemented with the functions block:

functions {
  /**
  * g() is the un-normalised generalised gamma distribution (GGD)
  */
  vector g(vector m, real Hs, real alpha, real beta){
    vector[num_elements(m)] y;
    vector[num_elements(m)] x;
    y <- m/Hs;
    for (i in 1:num_elements(m)) x[i] <- pow(y[i],beta);
    return log(beta) + alpha*log(y) - x;
  }

  /**
  * q() is the normalisation of the truncated generalised gamma distribution
  */
  real q(real mmin, real Hs, real logHs, real alpha, real beta){
    real z;

    z <- (alpha+1)/beta;
    // If z<0, gamma_q won't work !!!!
    return log(10)*logHs + log(gamma_q(z,pow(mmin/Hs,beta))) + lgamma(z);
  }

  /**
  * truncated_GGD_log is the log PDF of the lower-truncated GGD.
  */
  real truncated_GGD_log(vector m, real Hs, real logHs, real alpha, real beta, real mmin, vector weight){
    return sum((g(m,Hs,alpha,beta)-q(mmin,Hs, logHs, alpha,beta)) .*weight);
  }
}

Assume this functions block is included in every other stan program in this question. Note that a weight vector is passed to the truncated_GGD_log, which is multiplied by the underlying log-likelihoods before summing.

The rest of the model looks like this:

data {
  int<lower=0> N; // number of halos
  vector<lower=0>[N] log_m_meas; // mass of halos

  //PASS IN RELEVANT BOUNDS ON PARAMETERS
  real<lower=0> hs_min; //The lower bound of hs
  real<lower=0,upper=20> hs_max; //Upper bound of hs
  real<lower=-2,upper=0> alpha_min; //Lower bound of alpha
  real<lower=-2,upper=0> alpha_max; //Upper bound of alpha
  real<lower=0> beta_min; //Lower bound of beta
  real<lower=0> beta_max; //Upper bound of beta
}

transformed data {
  vector<lower=0>[N] weight;
  vector<lower=0>[N] m;
  for (i in 1:N) m[i] <- pow(10,log_m_meas[i]);
  weight <- m/mean(m); // <-- WEIGHTS ARE PROPORTIONAL TO MASSES
}

parameters {
  real<lower=beta_min,upper=beta_max> beta;
  real<lower=hs_min,upper=hs_max> logHs;
  real<lower=alpha_min,upper=alpha_max> alpha;
}

model {
  real Hs;
  real mmin;

  Hs <- pow(10,logHs);
  mmin <- pow(10,min(log_m_meas));

  // ENSURE WE PASS ALPHA+1, NOT JUST ALPHA
  m ~ full(Hs, logHs, alpha+1, beta, mmin, weight);
}

This works (ie. it converges to parameters that are pretty close to what I started with in a generative model). However, in reality, the masses are not known perfectly. They are measured with some uncertainty (for now, assume its standard log-normal with constant $\sigma$). To model this, I need to introduce a hierarchical model.

Assuming I didn't have the weighting problem, I would use the following:

data {
  int<lower=0> N; // number of halos
  vector<lower=0>[N] log_m_meas; // *measured* mass of halos
  real err; //uncertainty in halo masses

  // CONTROLS FOR PARAMETER BOUNDS
  real<lower=0> hs_min; //The lower bound of hs
  real<lower=0,upper=20> hs_max; //Upper bound of hs
  real<lower=-2,upper=0> alpha_min; //Lower bound of alpha
  real<lower=-2,upper=0> alpha_max; //Upper bound of alpha
  real<lower=0> beta_min; //Lower bound of beta
  real<lower=0> beta_max; //Upper bound of beta

  //ADDITIONAL PARAMETER BOUNDS    
  real<lower=0,upper=20>log_mtrue_max; // Upper bound of true masses
  real<lower=0> log_mmin_min; //Lower bound of log_mmin
  real<lower=0> log_mmin_max; //Upper bound of log_mmin

}

parameters {
  real<lower=beta_min,upper=beta_max> beta;
  real<lower=hs_min,upper=hs_max> logHs;
  real<lower=alpha_min,upper=alpha_max> alpha;
  real<lower=log_mmin_min,upper=log_mmin_max> log_mmin;
  vector<lower=log_mmin,upper=log_mtrue_max>[N] log_mtrue; //True mass estimates
}

transformed parameters {
  real Hs;
  real mmin;
  vector[N] m;
  Hs <- pow(10,logHs);
  mmin <- pow(10,log_mmin);
  for (i in 1:N) m[i] <- pow(10,log_mtrue[i]);
}

model {
  m ~ truncated_GGD(Hs, logHs, alpha, beta, mmin); // <-- USES A TRUNCATED_GGD definition with no weights

  //ALSO USE PRIORS ON LOG_MTRUE
  log_mtrue ~ normal(log_m_meas,err);
}

This also works, given I simulate data with an intrinsic $\alpha>-1$. However, to add the weighting, so I can use $\alpha<-1$, I have to add the following transformed data block:

transformed data {
  vector<lower=0>[N] weight;
  vector<lower=0>[N] m;
  vector<lower=0>[N] werr;
  for (i in 1:N) m[i] <- pow(10,log_m_meas[i]);
  weight <- m/mean(m); // <-- WEIGHTS ARE PROPORTIONAL TO MASSES
  for (i in 1:N) werr[i] <- err/sqrt(weight[i]);
}

and modify the model block to this:

model {
  // USE WEIGHTS AND ALPHA+1
  m ~ truncated_GGD(Hs, logHs, alpha+1, beta, mmin,weight);

  // USE WERR = ERR/SQRT(WEIGHT) instead of just ERR
  log_mtrue ~ normal(log_m_meas,werr);

}

And this does not work. Just in case it was something to do with intrinsically not being able to weight things this way, I tried swapping out the truncated_GGD for a Pareto distribution (which is the same except simpler, not having an exponential cutoff), using the same kind of machinery everywhere else. See the question I mentioned at the top. That works. This doesn't.

Specifically, $\log m_{\rm min}$ runs down into its lower limit, $\beta$ runs up into its upper limit, and $\alpha$ gets pushed more negative than it should be. $H_s$ is estimated significantly higher than it should be, but otherwise is at least constrained.

Interestingly, though $\log m_{\rm min}$ ends up at 11 (its lower limit) when it should be 12, the actual minimum 'true' mass estimated is 12.3. This makes absolutely no mathematical sense, as completely independent of everything else, $m_{\rm min}$ has highest likelihood at the minimum estimated mass.

This has been a long question, but I guess I want to know if there's some inherent conceptual flaw in all of this, or just some silly coding error, or perhaps a bug in stan (I'd say that the probability is decreasing in that order).

Answer 1

I still do not know how to specifically answer the question I asked, but have found a way around the situation, which I thought may be useful to anyone with the same question.

To avoid having to use the weights, the problem can be reformulated in terms of only positive shape parameters by using a recurrence relation: $$ \Gamma(z,x) = \frac{\Gamma(z+1,x)-x^z e^{-x}}{z}. $$ In my case, since $z$ is always expected to be $-2 < z < 0 $, I can use the following expanded relation in general: $$ \Gamma(z,x) = \frac{1}{z}\left(\frac{\Gamma(z+2,x)}{z+1} - x^z e^{-x}\left(\frac{x}{z+1}+1\right)\right), $$ which only ever needs to calculate $\Gamma$ with $z+2$, and therefore will always work in Stan (and any library that includes $\Gamma$ for $z>0$).

Furthermore, for those who are interested, if more arbitrary negative values of $z$ are required, a basic way to implement it in terms of only positive shape is (see eg. http://dlmf.nist.gov/8.8#E9) $$ \Gamma(z,x) = \Gamma(z)\left(Q(z+n,x) - x^z e^{-x} \sum_{k=0}^{n-1} \frac{x^k}{\Gamma(z+k+1)}\right), $$ where $Q$ is the upper regularized gamma. Though this involves (complete) gamma functions with negative parameters, this is usually implemented (and certainly this is the case for Stan).

Obviously, we just set $n$ to be $\lceil -z \rceil$. The further $z$ is from 0, the longer the calculation will take, and the more residual error there will be.

Specifically, in Stan, my final functions block looks like this (I only implemented the form for $z>-2$, since it is faster for my application):

functions {       
/**
* g() is the shape part of the generalised gamma distribution (GGD)
* @param vector m, masses 
* @param real Hs, turn-over mass
* @param real alpha, power-law slope
* @param real beta, cut-off parameter

* @returns the un-normalised GGD
*/
vector g(vector m, real Hs, real alpha, real beta){
    vector[num_elements(m)] y;
    vector[num_elements(m)] x;

    y <- m/Hs;
    for (i in 1:num_elements(m)) x[i] <- pow(y[i],beta);
    return log(beta) + alpha*log(y) - x;
}

/**
* gammainc() is the upper incomplete gamma function (not regularized)
* @param real a, shape parameter, must be >0
* @param real x, position, >0
* @returns the non-regularized incomplete gamma
*/
real gammainc(real a, real x){
     return gamma_q(a,x) * tgamma(a);
}

/**
* gammainc_neg() is the upper incomplete gamma function for a > -2
* @param real a, shape parameter, must be >-2
* @param real x, position, >0
* @returns the non-regularized incomplete gamma
*
* NOTES: this routine has been written for speed for this application.
*        More general functions for negative a can be gotten by recursion
*        formulae etc.
*/
real gammainc_neg(real a, real x){
 return ((gammainc(a+2,x)-pow(x,(a+1))*exp(-x))/(a+1)-exp(-x)*pow(x,a))/a;
}

/**
* q() is the normalisation of the truncated generalised gamma distribution
* @param real mmin, the truncation mass
* @param real Hs, the turnover mass
* @param real logHs, log10 of the turnover mass
* @param real alpha, the power-law slope
* @param real beta, the cut-off parameter
* @returns the normalisation of the truncated GGD
*/ 
real q(real mmin, real Hs, real logHs, real alpha, real beta){
    real z;

z <- (alpha+1)/beta;
    return log(10)*logHs + log(gammainc_neg(z,pow((mmin/Hs),beta)));
}


/**
* truncated_GGD_log gives the log PDF of the lower-truncated generalised gamma distribution.
* @param vector m, masses
* @param real Hs, turnover mass
* @param real logHs, log10 of turnover mass
* @param real alpha, power-law slope
* @param real beta, cut-off parameter
* @param real mmin, truncation mass
* @param vector weight, weights of masses
*/
real truncated_GGD_log(vector m, real Hs, real logHs, real alpha, real beta, real mmin){
  return sum(g(m,Hs,alpha,beta)-q(mmin,Hs, logHs, alpha,beta));
}
}