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

**URL:** <https://discourse.mc-stan.org/t/modeling-mixture-of-two-gamma-knowing-that-the-one-with-the-highest-value-is-the-output/9396>\
**Category:** Modeling\
**Created:** [June 25, 2019, 12:45pm UTC](https://discourse.mc-stan.org/t/modeling-mixture-of-two-gamma-knowing-that-the-one-with-the-highest-value-is-the-output/9396 "2019-06-25T12:45:31Z")\
**Posts on this page:** 1\
**Page:** 2

<div class="post-metadata">

**Author:** ![maxbiostat](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/maxbiostat/32/1326_2.png) [@maxbiostat](https://discourse.mc-stan.org/u/maxbiostat)\
**Post date:** [July 4, 2019, 3:43pm UTC](https://discourse.mc-stan.org/t/modeling-mixture-of-two-gamma-knowing-that-the-one-with-the-highest-value-is-the-output/9396/21 "2019-07-04T15:43:24Z")

</div>

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

```stan
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]);
}

```

[Previous page](https://discourse.mc-stan.org/t/modeling-mixture-of-two-gamma-knowing-that-the-one-with-the-highest-value-is-the-output/9396.md?page=1)