MVN_GibbsSampler: Gibbs sampler for MVN distribution

Description

Generating random vectors on the basis of a given MVN distribution, through Gibbs sampling algorithm.

Usage

# Bayesian posteriori as data
# data <- MVN_BayesianPosteriori(dataset1)
# Using Gibbs sampler to generate random vectors based on Bayesian posteriori:
MVN_GibbsSampler(n, data, initial, reject_rate, burn)

Arguments

A positive integer. The numbers of random vectors to be generated.

data

A double level list which contains the mean vector (data$mean) and the covariance matrix (data$var) of a given MVN distribution.

initial

Initial vector where Markov chain starts. Default value use a random vector generated by rmvnorm().

reject_rate

A numeric to control burn-in period by ratio. Default value is 0.2, namely the first 20% generated vectors will be rejected. If this arg was customized, the next arg burn should maintain the default value.

burn

A numeric to control burn-in period by numbers. If this arg was customized, final result will be generated by this manner in which it will drop the first n numbers (n=burn).

Value

return a series random vectors in the basis of given MVN distribution.

Details

There're also some literatures suggest using the mean or mode of priori as initial vector. Users can customize this setting according to their own needs.

Examples

Run this code

# NOT RUN {
library(mvtnorm)

# Get parameters of Bayesian posteriori multivariate normal distribution
BPos <- MVN_BayesianPosteriori(dataset1)
BPos

# Using previous result (BPos) to generate random vectors through Gibbs
# sampling: 7000 observations, start from c(1,1,2), use 0.3 burning rate
BPos_Gibbs <- MVN_GibbsSampler(7000, BPos, initial=c(1,1,2), 0.3)
tail(BPos_Gibbs)

# Check for convergence of Markov chain
BPos$mean
colMeans(BPos_Gibbs)
BPos$var
var(BPos_Gibbs)
# }
# NOT RUN {
# 3d Visulization:
library(rgl)
fac1 <- BPos_Gibbs[,1]
fac2 <- BPos_Gibbs[,2]
fac3 <- BPos_Gibbs[,3]
plot3d(x=fac1, y=fac2, z=fac3, col="red", size=2)
# }