Likelihood of (bivariate?) Bernoulli distribution In 50 individuals we measure the presence of 2 SNP's ("single nucleotide polymorphisms") and 1 trait. We assume that both SNP's and the trait are binary. We assume the model $Y_{ijk} = \text{Bernoulli}(p_{ij})$, where $$Y_{ijk} = 1\text{ if person }k\text{ with SNP 1 }= i \text{ and SNP 2 }= j\quad (i,j\in \{0,1\})\text{ has trait 1}$$ and $Y_{ijk} = 0$ otherwise. 
I want to write down the likelihood of the data but I'm not quite sure what it looks like. Is it true that $P(Y_{ijk} = y) = p_{ij}^y(1-p_{ij})^y$, $y =0,1$, or is this wrong? On the internet I find only one dimensional distributions so I just want to check if the density is correct.    
 A: The likelihood is close, but not quite correct. Your Bernoulli assumption at the top establishes
\begin{gather}
\Pr[Y_{ijk}=1] = p_{ij}
\end{gather}
so it follows that
\begin{gather}
\Pr[Y_{ijk}=0] = 1-p_{ij}
\end{gather}
Combining these two, the likelihood of a single observation can be written as
\begin{gather}
\rho(y_{ijk}) = p_{ij}^{y_{ijk}}(1-p_{ij})^{(1-y_{ijk})}
\end{gather}
Note that the fact that $i$ and $j$ are coded as binary variables does not matter here -- they simply index observations. (It does, however, matter that the $y$'s are binary). 
Finally, if the observations are independent, the likelihood of all the $y_{ijk}$'s is
\begin{gather}
p(y) = \prod_{i=0}^1\prod_{j=0}^1\prod_k \rho(y_{ijk})
\end{gather}
You can also sum up all $y_{ijk}$'s for a given pair $(i,j)$---the distribution of the sum is binomial---leaving you with the product of four binomial distributions:
\begin{gather}
y_{ij} = \sum_{k=1}^N y_{ijk}
\end{gather}
\begin{gather}
p(y) = \prod_{i=0}^1\prod_{j=0}^1 Binom( y_{ij} | N, p_{ij} )
\end{gather}
