Modeling mixture of two gamma knowing that the one with the highest value is the output

Modeling
21

Right, think I’ve found something mildly usable. First, theory. Recall that we have Y = \max(Y_1, Y_2), with Y_i \sim \text{gamma}(\text{shape}= \alpha_i, \text{rate}= 1). Under independence of Y_1 and Y_2 we have Pr(Y \leq y) = \prod_i Pr(Y_i \leq y) = \Pi_1(y)\Pi_2(y). Differentiating with respect to y we get the probability density function f_Y(y) = \pi_1(y)\Pi_2(y) + \pi_2(y)\Pi_1(y), i.e., what you would get with your initial approach and log_mix if you fixed \theta = 1/2.

The following R code produces usable results for a good range of values for \alpha_1 and \alpha_2. It “fails” if you do something like \alpha_1 = 18 and \alpha_2 = 20, although it will correctly recover the value for \alpha_2. Following @betanalpha’s recommendation, I have imposed an ordering constraint on the shape parameters, so \alpha_2 will now always be the largest.

#### Notice: now Y2 is *necessarily* the largest of Y1, Y2. This is not how the original model was formulated.
N <- 500
alpha1 <- 3
alpha2 <- 6
y1 <- rgamma(N, alpha1, rate = 1)
y2 <- rgamma(N, alpha2, rate = 1)
z <- y2 > y1
mean(z)
( true.P1 <- extraDistr::pbetapr(1, alpha1, alpha2) )
gammaMax.data <- list(N=N, y = pmax(y1, y2))

hist(gammaMax.data$y, probability = TRUE)
curve(dgamma(x, shape = alpha1), min(gammaMax.data$y), max(gammaMax.data$y), col = 2, lwd = 2 , lty = 2, add = TRUE)
curve(dgamma(x, shape = alpha2), min(gammaMax.data$y), max(gammaMax.data$y), col = 3, lwd = 2 , lty = 2, add = TRUE)

library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = 4)
gamma.max.model <- stan_model("gamma_max.stan")

opt <- optimizing(gamma.max.model, data = gammaMax.data)
opt$theta_tilde
mcmc <- sampling(gamma.max.model, data = gammaMax.data)
mcmc
stan_trace(mcmc)

with the Stan program

functions {
  real betapr_cdf(real x, real alpha, real beta){
    return inc_beta(alpha, beta, x/(1+x));
  }
}
data {
  int N;
  real<lower=0> y[N];
}
parameters {
  ordered[2] alpha;
}
model {
  alpha ~ normal(0.5, 2);
  for (n in 1:N){
    target += log_sum_exp(gamma_lpdf(y[n]|alpha[1], 1) + gamma_lcdf(y[n] | alpha[2], 1),
    gamma_lpdf(y[n]|alpha[2], 1) + gamma_lcdf(y[n] | alpha[1], 1));
  }
}
generated quantities{
   real <lower= 0, upper=1> P1 = betapr_cdf(1.0, alpha[1], alpha[2]);
}
1 Like