Let’s say that *p* persons are throwing a basketball into a set of *b* baskets which have different radii and the outcome, Spb is binary, indicating whether the basket was scored or not. Let *ap* be person *p*’s throwing accuracy and *rb* the radius of the *b*th basket. The model could be formulated in Stan as such:

```
model{
S[p,b] ~ bernoulli_logit(a[p] + r[b])
}
```

However, what if we wanted to include correlations between trials? Let’s say that scoring a basket gives you a confidence boost which in turn increases your probability of scoring the next hits. Likewise, failing to score a basket makes you more anxious and decreases your chances of subsequent scores. So if Sp1 = 1, the chances that Sp2 or Sp3 will also be 1 are increased.

How could we code these correlations in Stan? If I was dealing with continuous variables, I guess I would use the `multi_normal_lpdf`

function and pass the covariances between trials via the `Sigma`

argument, but I don’t see a way to do that with binary variables since the `bernoulli_logit`

or other related functions don’t have arguments related to the covariance structure.

Furthermore, I don’t have a feeling that simply using the binomial distribution would estimate these correlations. I may be incorrect with this, though. Any help would be appreciated.