# Binomial coefficient in JAGS

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

• Thank you! Sorry but I'm pretty new with this – 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])) )