# Example 1
(myrho <- rhobit(2, inverse = TRUE))
ymat = rbinom2.rho(nn <- 2000, mu1 = 0.8, rho = myrho, exch = TRUE)
(mytab = table(ymat[,1], ymat[,2], dnn = c("Y1","Y2")))
fit = vglm(ymat ~ 1, binom2.rho(exch = TRUE))
coef(fit, matrix = TRUE)
# Example 2
bdata = data.frame(x = sort(runif(nn)))
bdata = transform(bdata, mu1 = probit(-2+4*x, inverse = TRUE),
mu2 = probit(-1+3*x, inverse = TRUE))
dmat = with(bdata, dbinom2.rho(mu1, mu2, myrho))
ymat = with(bdata, rbinom2.rho(nn, mu1, mu2, myrho))
fit2 = vglm(ymat ~ x, binom2.rho, bdata)
coef(fit2, matrix = TRUE)
matplot(with(bdata, x), dmat, lty = 1:4, col = 1:4,
type = "l", main = "Joint probabilities",
ylim = 0:1, lwd = 2, ylab = "Probability")
legend(x = 0.25, y = 0.9, lty = 1:4, col = 1:4, lwd = 2,
legend = c("1 = (y1=0, y2=0)", "2 = (y1=0, y2=1)",
"3 = (y1=1, y2=0)", "4 = (y1=1, y2=1)"))
Run the code above in your browser using DataLab