Parameter estimation for the sum of negative binomial distributions

Modeling
10

First, a big thank you to @martinmodrak for your insight. I read your blogpost and tried to follow your lead for my model. Here’s the code for reference, and hopefully it’ll help anyone else working on a similar problem:

functions {
  /* Using the method of moments approximations for computing mu and phi overall.
   * This approach follows that in Martin Modrak's blog (martinmodrak.cz).
   */
  real negBinomialSumMoments_lpmf(int[] sumY, vector mu, vector phi) {
    real muOverall = sum(mu);
    real phiOverall = square(muOverall) / sum(square(mu) ./ phi);
    
    return neg_binomial_2_lpmf(sumY | muOverall, phiOverall);
  }
}

data {
  // These are the input data sizes, followed by the actual observed data.
  int<lower=0> NX;
  int<lower=0> NY;
  int<lower=0> NZ;
  
  int yX[NX];
  int yY[NY];
  int yZ[NZ];
}

parameters {
  positive_ordered[3] mu;
  vector<lower=0>[3] cv; // coefficient of variation
}

transformed parameters {
  vector[3] phi;
  vector[3] variances;
  
  phi[1] = inv(square(cv[1]));
  phi[2] = inv(square(cv[2]));
  phi[3] = inv(square(cv[3]));
  
  variances = mu + square(mu) ./ phi;
}

model {
  cv[1] ~ exponential(1);
  cv[2] ~ exponential(1);
  cv[3] ~ exponential(1);
  
  mu[1] ~ gamma(0.01, 0.01);
  mu[2] ~ gamma(0.01, 0.01);
  mu[3] ~ gamma(5, 10);
  
  yX ~ neg_binomial_2(mu[1], phi[1]);
  yY ~ negBinomialSumMoments(mu[:2], phi[:2]);
  yZ ~ negBinomialSumMoments(mu, phi);
}

I used 20,000 iterations and 5 chains. The results are not really what I hoped for though. The main issue is that the means are pretty much identical and it seems that the \phi_i capture the differences among the three NB distributions. My assumption is that the mean of y_X should be small, y_Y moderate, and y_Z highest. Hence, the gamma priors I chose. It’s very possible my implementation isn’t what you had in mind, but please let me know if you (@martinmodrak) or anyone else, has any suggestions or corrections.

I appreciate all the help as always.

1 Like