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

 A: 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:


*

*Using the logfact function

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

