# Gamma process - is this feasible?

**URL:** <https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544>\
**Category:** Modeling\
**Tags:** bioinformatics, techniques\
**Created:** [May 15, 2017, 2:41am UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544 "2017-05-15T02:41:46Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![stemangiola](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stemangiola/32/260_2.png) [@stemangiola](https://discourse.mc-stan.org/u/stemangiola)\
**Post date:** [May 15, 2017, 2:41am UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/1 "2017-05-15T02:41:46Z")

</div>

Hello,

I am trying to developed a gamma generating process

negative\_binomial -\> gamma -\> hyper\_prior

and I found this proposed solution. Has someone approached this process before? Do you think there is the foundation for creating (myself) a custom density function in stan

![](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/1X/cbd32c8ea7d7dfa50733a5431792707662961039.png)

Few useful links:

- [http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.157.5540&rep=rep1&type=pdf](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.157.5540&rep=rep1&type=pdf)
- [https://stats.stackexchange.com/questions/2635/conjugate-prior-for-a-gamma-distribution](https://stats.stackexchange.com/questions/2635/conjugate-prior-for-a-gamma-distribution)

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [May 15, 2017, 7:18pm UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/2 "2017-05-15T19:18:27Z")

</div>

All those functions are built into Stan, so you can just write your own `_lpdf` as a Stan function. It’d be more efficient to do it in C++ with vectorization and analytical gradients, but it shouldn’t be bad encoded directly in Stan. You might be able to code it up in terms of sufficient statistics to make things faster.

No, I don’t have experience with this distribution.

---

<div class="post-metadata">

**Author:** ![stemangiola](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stemangiola/32/260_2.png) [@stemangiola](https://discourse.mc-stan.org/u/stemangiola)\
**Post date:** [May 16, 2017, 12:18am UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/3 "2017-05-16T00:18:26Z")

</div>

Thanks Bob,

I don’t need to add the normalisation constant right?

I will implement in Stan. @sakrejda do you have any general design suggestion before I dive into that?

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [May 16, 2017, 12:28am UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/4 "2017-05-16T00:28:31Z")

</div>

Nope, I haven’t tried something like this so my only suggestion is to make sure to log and simplify on paper before you try to code it (sort of a duh, but I thought I’d say it, never write a density that has a tgamma in it) but it looks like fun stuff. What kind of data are you using this for?

---

<div class="post-metadata">

**Author:** ![stemangiola](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stemangiola/32/260_2.png) [@stemangiola](https://discourse.mc-stan.org/u/stemangiola)\
**Post date:** [May 16, 2017, 12:35am UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/5 "2017-05-16T00:35:27Z")

</div>

It is gene expression (counts) data. And the process follows nb -\> gamma -\>  
hyperprior with shape and rate unknown. A nightmare so far.

Do you have anything to say about sufficient statistics looking at the book  
image i posted (page 25). I’m quite blank on that.

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [May 16, 2017, 12:41am UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/6 "2017-05-16T00:41:36Z")

</div>

Can you be clearer about which part is confusing? It seems pretty clear to me but maybe I’m missing the point (?) The sufficient statistics S and P are defined on p. 25, once you log you get:

```
(a-1) * sum_i=1^n (log(x_i)) - log b - log (sum_i=1^n x_i) - n (lgamma(a) - a * log(b)) 

```

That makes it pretty efficient.

---

<div class="post-metadata">

**Author:** ![stemangiola](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stemangiola/32/260_2.png) [@stemangiola](https://discourse.mc-stan.org/u/stemangiola)\
**Post date:** [May 16, 2017, 1:24am UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/7 "2017-05-16T01:24:23Z")

</div>

(shouldn’t it be `[b*sum(x)]` instead of `[log b - log (sum_i=1^n x_i)]` ? )

Just to make sure I understand:

the equation 96 is should be used in my case instead of gamma as prior for the mean of the negative\_binomial\_2, and the equation 97 should be used as prior of alphas and betas from the equation 96

negative\_binomial\_2[n,i] -\> 96[i] -\> 97

OR

negative\_binomial\_2[n,i] -\> gamma[i] -\> 96

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [May 16, 2017, 2:02am UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/8 "2017-05-16T02:02:15Z")

</div>

(Yes you’re right)

> [@stemangiola](#):
>
> the equation 96 is should be used in my case instead of gamma as prior for the mean of the negative\_binomial\_2, and the equation 97 should be used as prior of alphas and betas from the equation 96

yes

> [@stemangiola](#):
>
> negative\_binomial\_2[n,i] -\> 96[i] -\> 97
> 
> OR
> 
> negative\_binomial\_2[n,i] -\> gamma[i] -\> 96

I don’t understand what you’re asking. (?)

---

<div class="post-metadata">

**Author:** ![stemangiola](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stemangiola/32/260_2.png) [@stemangiola](https://discourse.mc-stan.org/u/stemangiola)\
**Post date:** [May 16, 2017, 2:05am UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/9 "2017-05-16T02:05:24Z")

</div>

Thanks,

simply this

`negative_binomial_2[n,i] -> 96[i] -> 97`

is same as

_“the equation 96 is should be used in my case instead of gamma as prior for the mean of the negative\_binomial\_2, and the equation 97 should be used as prior of alphas and betas from the equation 96”_

The alternative is

negative\_binomial\_2[n,i] -\> gamma[i] -\> 96

being a gamma used as prior for NB and 96 used for prior for gamma, but I realise now that this doesn’t make sense. :)

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [May 16, 2017, 5:41pm UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/10 "2017-05-16T17:41:17Z")

</div>

You need the normalization constant if it depends on parameters. You can drop additive terms in the log density that only depend on data and literals (that is, they don’t depend on parameters).

---

<div class="post-metadata">

**Author:** ![stemangiola](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stemangiola/32/260_2.png) [@stemangiola](https://discourse.mc-stan.org/u/stemangiola)\
**Post date:** [May 16, 2017, 11:44pm UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/11 "2017-05-16T23:44:25Z")

</div>

Thanks a lot.

That was a good call

---

<div class="post-metadata">

**Author:** ![stemangiola](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/stemangiola/32/260_2.png) [@stemangiola](https://discourse.mc-stan.org/u/stemangiola)\
**Post date:** [May 17, 2017, 12:42am UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/12 "2017-05-17T00:42:48Z")

</div>

@Bob_Carpenter do you have any suggestion on how to approach this normalising part? Integral of a function including gamma is a bit beyond my expertise…

![](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/1X/5b85cd94e070b146d8cb537706f573e7e62fdeee.png)

---

<div class="post-metadata">

**Author:** ![wds15](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/wds15/32/908_2.png) [@wds15](https://discourse.mc-stan.org/u/wds15)\
**Post date:** [May 17, 2017, 11:57am UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/13 "2017-05-17T11:57:21Z")

</div>

You could try to throw this integral onto Mathematica.

Last option is to use the ode integrator in Stan to calculate the constant using numerical integration - won’t be fast but it can be done.

---

<div class="post-metadata">

**Author:** ![sakrejda](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sakrejda/32/845_2.png) [@sakrejda](https://discourse.mc-stan.org/u/sakrejda)\
**Post date:** [May 17, 2017, 12:52pm UTC](https://discourse.mc-stan.org/t/gamma-process-is-this-feasible/544/14 "2017-05-17T12:52:39Z")

</div>

Good point about Mathematica, I recently did some similar ones so you’re likely to end up with some incomplete gamma functions. If you do, those are gamma\_p in Stan and you likely will need a bugfix that’s in-progress.