top of page

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)

â€‹================================================

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

bottom of page