Binomial coefficient in JAGS

Question

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)
  }
}

Answer 1

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])) )