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

(c) 2017,  Joseph M. Hilbe, Rafael S. de Souza and Emille E. O. Ishida  

 

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)

​================================================

 

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