Simulate the sum of $n$ dependent bernoulli random variables.
I know that the sum of $n$ independent bernoulli radom variables with probability $p$ is distributed as a binomial random variable. I found that there are multiple models that can be usefull to approximate the sum when the bernoulli random variables are dependent with correlation $q$. In this case I'm useing this one: $$f_X(x) = \binom{n}{x}p^x(1-p)^{n-x}e^{-x(n-x)q}$$$
I´m trying to simulate it in R but I don´t know how to mix this two in a random sample. Any suggestions would be great!
Here is my code so far:
sumaBerdep <- function (n,p,q)
{
x <- rbern(n,p)
y <- sample(n, replace = T)
return(sum(x)*exp(-y*(n-y)*q))
}
binom <- replicate(n, sumaBerdep(10, 0.5, 0.3))
p <- table(binom)/10^5
plot(p, type = 'h')
