I'm trying to write such likelihood using JAGS in R (in order to estimate parameters $\xi$ and $\pi$, being $m$ fixed equal 7) using the "ones trick" and got stuck in defining the binomial coefficient (term3). How would you define the binomial coefficient? $$ P(Y=y)=\left\{ \pi\binom{m-1}{y-1} (1-\xi)^{y-1}\xi^{m-r}+(1-\pi){1\over m}\right\} $$

This is my model:

model {
  for (i in 1:N) {    
    term1[i] <- pai[i]*pow(1-csi[i],(y[i]-1))*pow(csi[i],(m-y[i]))
    term2[i] <- (1-pai[i])*(1/m)
    term3[i] <- logfact(m-1)/logfact(m-y[i]-1)*logfact(y[i])

    L[i] <- term3[i]*term1[i]+term2[i]
    phi[i] <- L[i]/C
    ones[i] ~ dbern(phi[i]) 

    csi[i] ~ dbeta (alfa, beta)
    pai[i] ~ dunif (0,1)
  • $\begingroup$ Thank you! Sorry but I'm pretty new with this $\endgroup$ – user27115 Jun 23 '13 at 23:06

The logarithm of the binomial coefficient can be implemented in WinBUGS/JAGS by using the function logfact which is the logarithm of the factorial: $\ln(x!)$. Alternatively, we could use the function of the logarithm of the gamma function $\ln(\Gamma(x))$ loggam.

The binomial coefficient is defined as: $$ \binom{n}{k} = \frac{n!}{k!(n-k)!} $$ Using the logarithmic identities: $$ \log_{b}(x\cdot y)=\log_{b}(x) + \log_{b}(y) $$ and $$ \log_{b}\left(\frac{x}{y}\right)=\log_{b}(x) - \log_{b}(y) $$ we can write the binomial coefficient $\binom{m-1}{y-1}$ as

$$ \begin{align} \binom{m-1}{y-1} &=\frac{(m-1)!}{(y-1)!(m-y)!} \\ &=\exp \left\{ \ln\left[(m-1)!\right]-\left(\ln\left[(y-1)!\right]+\ln\left[(m-y)!\right] \right) \right\} \\ &=\frac{\Gamma(m)}{\Gamma(1+m-y)\Gamma(y)}\\ &=\exp \left\{ \ln\left[\Gamma(m)\right]-\left(\ln\left[\Gamma(1+m-y)\right]+\ln\left[\Gamma(y)\right] \right) \right\} \\ \end{align} $$

So we have two possibilities to implement the binomial coefficient:

  1. Using the logfact function
  2. Using the loggam function (i.e. the logarithm of the gamma function, $\ln(\Gamma(x))$

The first implementation in JAGS/WinBUGS is:

term3 <- exp( logfact(m-1) - (logfact(y[i]-1) + logfact(m-y[i])) )

and the second using the log-gamma function:

term3 <- exp( loggam(m) - (loggam(1+m-y[i]) + loggam(y[i])) )

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.