Data augmentation for incomplete multivariate normal data
Data augmentation under a normal-inverted Wishart prior. If no prior is specified by the user, the usual "noninformative" prior for the multivariate normal distribution is used. This function simulates one or more iterations of a single Markov chain. Each iteration consists of a random imputation of the missing data given the observed data and the current parameter value (I-step), followed by a draw from the posterior distribution of the parameter given the observed data and the imputed data (P-step).
da.norm(s, start, prior, steps=1, showits=FALSE, return.ymis=FALSE)
s |
summary list of an incomplete normal data matrix produced by the
function |
start |
starting value of the parameter. This is a parameter vector in packed
storage, such as one created by the function |
prior |
optional prior distribution. This is a list containing the
hyperparameters of a normal-inverted Wishart distribution. In order,
the elements of the list are: tau (a scalar), m (a scalar), mu0 (a
vector of length |
steps |
number of data augmentation iterations to be simulated. |
showits |
if |
return.ymis |
if |
if return.ymis=FALSE
, returns a parameter vector, the result of the last
P-step. If the value of steps
was large enough to guarantee
approximate stationarity, then this parameter can be regarded as a
proper draw from the observed-data posterior, independent of start
.
If return.ymis=TRUE
, then this function returns a list of the following
two components:
parameter |
a parameter vector, the result of the last P-step |
ymis |
a vector of missing values, the result of the last I-step. The length
of this vector is |
Before this function may be used, the random number generator seed
must be initialized with rngseed
at least once in the current S
session.
See Chapter 5 of Schafer (1996).
rngseed
, em.norm
, prelim.norm
, and getparam.norm
.
data(mdata) s <- prelim.norm(mdata) thetahat <- em.norm(s) #find the MLE for a starting value rngseed(1234567) #set random number generator seed theta <- da.norm(s,thetahat,steps=20,showits=TRUE) # take 20 steps getparam.norm(s,theta) # look at result
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.