Probit model using R

From: Bayesian Models for Astrophysical Data, Cambridge Univ. Press

you are kindly asked to include the complete citation if you used this material in a publication

Code 5.20 - Synthetic probit data and model generated in R

================================================

set.seed(135)
nobs <- 1:2000

x1 <- runif(nobs)
x2 <- 2*runif(nobs)

xb <- 2 + 0.75 * x1 - 1.25 * x2

exb <- pnorm(xb) # probit inverse link

py <- rbinom(nobs, size=1, prob=exb)

probdata <- data.frame(py, x1, x2)

================================================

Code 5.21 Probit model using R

================================================

library(MCMCpack)

mypL <- MCMCprobit(py ~ x1 + x2,
burnin = 5000,
mcmc = 100000,
data = probdata)

summary(mypL)
plot(mypL)

================================================

Code 5.21

Output on screen:

Iterations = 5001:105000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 1e+05

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE Time-series SE

(Intercept) 1.9391 0.10794 0.0003413 0.0008675

x1 0.7185 0.12486 0.0003948 0.0008250

x2 -1.2496 0.07079 0.0002239 0.0006034

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

(Intercept) 1.7302 1.8662 1.9381 2.0119 2.1535

x1 0.4756 0.6339 0.7185 0.8022 0.9645

x2 -1.3900 -1.2970 -1.2489 -1.2017 -1.1120

HSI