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 3.3 Bayesian normal linear model in R.
==================================================
library(MCMCpack)


# Data
nobs = 5000                                                  # number of obs in model
x1 <- runif(nobs)                                          # random uniform variable
beta0 = 2.0                                                   # intercept
beta1 = 3.0                                                   # angular coefficient
xb <- beta0 + beta1*x1                                # linear predictor
y <- rnorm(nobs, xb, sd=1)

 

# Fit
posteriors <- MCMCregress(y ~ x1, thin=1, seed=1056, burnin=1000,
                                               mcmc=10000, verbose=1)

 

# Output
summary(posteriors)

 

# Plot
plot(posteriors)
==================================================

Output on screen:

Iterations = 1001:11000

Thinning interval = 1

Number of chains = 1

Sample size per chain = 10000

 

1. Empirical mean and standard deviation for each variable,

    plus standard error of the mean:

                      Mean                SD         Naive SE         Time-series SE

(Intercept)     1.980         0.02862        0.0002862                0.0002862

x1                  3.034         0.04945        0.0004945                0.0004945

sigma2          1.011         0.02002         0.0002002               0.0002052

2. Quantiles for each variable:

                        2.5%              25%           50%           75%          97.5%

(Intercept)    1.9239          1.9607          1.980          1.999          2.036

x1                 2.9391          2.9997          3.034          3.067          3.133

sigma2          0.9727         0.9973          1.011          1.025           1.051

© 2017 by Emille E. O. Ishida