R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### * <HEADER>
> ###
> attach(NULL, name = "CheckExEnv")
> assign(".CheckExEnv", as.environment(2), pos = length(search())) # base
> ## add some hooks to label plot pages for base and grid graphics
> setHook("plot.new", ".newplot.hook")
> setHook("persp", ".newplot.hook")
> setHook("grid.newpage", ".gridplot.hook")
> 
> assign("cleanEx",
+        function(env = .GlobalEnv) {
+ 	   rm(list = ls(envir = env, all.names = TRUE), envir = env)
+            RNGkind("default", "default")
+ 	   set.seed(1)
+    	   options(warn = 1)
+ 	   delayedAssign("T", stop("T used instead of TRUE"),
+ 		  assign.env = .CheckExEnv)
+ 	   delayedAssign("F", stop("F used instead of FALSE"),
+ 		  assign.env = .CheckExEnv)
+ 	   sch <- search()
+ 	   newitems <- sch[! sch %in% .oldSearch]
+ 	   for(item in rev(newitems))
+                eval(substitute(detach(item), list(item=item)))
+ 	   missitems <- .oldSearch[! .oldSearch %in% sch]
+ 	   if(length(missitems))
+ 	       warning("items ", paste(missitems, collapse=", "),
+ 		       " have been removed from the search path")
+        },
+        env = .CheckExEnv)
> assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now
> assign("ptime", proc.time(), env = .CheckExEnv)
> grDevices::postscript("norm-Examples.ps")
> assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv)
> options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly"))
> options(warn = 1)    
> library('norm')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "da.norm"
> 
> ### * da.norm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: da.norm
> ### Title: Data augmentation for incomplete multivariate normal data
> ### Aliases: da.norm
> ### Keywords: distribution
> 
> ### ** Examples
> 
> data(mdata)
> s  <-  prelim.norm(mdata)
> thetahat <- em.norm(s)   #find the MLE for a starting value
Iterations of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...
> rngseed(1234567)   #set random number generator seed
> theta <- da.norm(s,thetahat,steps=20,showits=TRUE)  # take 20 steps
Steps of Data Augmentation: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...
> getparam.norm(s,theta) # look at result
$mu
[1]    46.047619    40.786778    11.757231 67293.308547     2.924422

$sigma
             [,1]          [,2]        [,3]         [,4]         [,5]
[1,]    306.83519    214.387810   -9.455950 3.726833e+05    22.284403
[2,]    214.38781    163.730635   -7.409176 2.839411e+05    15.448049
[3,]     -9.45595     -7.409176   20.124269 9.279856e+03    -2.711537
[4,] 372683.31550 283941.066595 9279.855659 8.152364e+08 31184.056069
[5,]     22.28440     15.448049   -2.711537 3.118406e+04     2.547689

> 
> 
> 
> cleanEx(); ..nameEx <- "em.norm"
> 
> ### * em.norm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: em.norm
> ### Title: EM algorithm for incomplete normal data
> ### Aliases: em.norm
> ### Keywords: regression
> 
> ### ** Examples
> 
> data(mdata)
> s <- prelim.norm(mdata)   #do preliminary manipulations
> thetahat <- em.norm(s)   #compute mle
Iterations of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...
> getparam.norm(s,thetahat,corr=TRUE)$r #look at estimated correlations
            [,1]       [,2]        [,3]      [,4]       [,5]
[1,]  1.00000000 0.88753257 -0.09642342 0.4889176  0.7145953
[2,]  0.88753257 1.00000000  0.09345297 0.5197721  0.5429959
[3,] -0.09642342 0.09345297  1.00000000 0.2930083 -0.3598209
[4,]  0.48891759 0.51977207  0.29300826 1.0000000  0.3029066
[5,]  0.71459529 0.54299589 -0.35982094 0.3029066  1.0000000
> 
> 
> 
> cleanEx(); ..nameEx <- "getparam.norm"
> 
> ### * getparam.norm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: getparam.norm
> ### Title: Extract normal parameters from packed storage
> ### Aliases: getparam.norm
> ### Keywords: regression
> 
> ### ** Examples
> 
> data(mdata)
> s <- prelim.norm(mdata)   #do preliminary manipulations
> thetahat <- em.norm(s)   #compute MLE
Iterations of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...
> getparam.norm(s,thetahat,corr=TRUE)$r #look at estimated correlations
            [,1]       [,2]        [,3]      [,4]       [,5]
[1,]  1.00000000 0.88753257 -0.09642342 0.4889176  0.7145953
[2,]  0.88753257 1.00000000  0.09345297 0.5197721  0.5429959
[3,] -0.09642342 0.09345297  1.00000000 0.2930083 -0.3598209
[4,]  0.48891759 0.51977207  0.29300826 1.0000000  0.3029066
[5,]  0.71459529 0.54299589 -0.35982094 0.3029066  1.0000000
> 
> 
> 
> cleanEx(); ..nameEx <- "imp.norm"
> 
> ### * imp.norm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: imp.norm
> ### Title: Impute missing multivariate normal data
> ### Aliases: imp.norm
> ### Keywords: regression
> 
> ### ** Examples
> 
> data(mdata)
> s <- prelim.norm(mdata)   #do preliminary manipulations
> thetahat <- em.norm(s)   #find the mle
Iterations of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...
> rngseed(1234567)   #set random number generator seed
> ximp <- imp.norm(s,thetahat,mdata)  #impute missing data under the MLE
> 
> 
> 
> cleanEx(); ..nameEx <- "loglik.norm"
> 
> ### * loglik.norm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: loglik.norm
> ### Title: Observed-data loglikelihood for normal data
> ### Aliases: loglik.norm
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> data(mdata)
> s <- prelim.norm(mdata)   #do preliminary manipulations
> thetahat <- em.norm(s)   #compute MLE
Iterations of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...
> loglik.norm(s,thetahat)  #loglikelihood at the MLE
[1] -21.62696
> 
> 
> 
> cleanEx(); ..nameEx <- "logpost.norm"
> 
> ### * logpost.norm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: logpost.norm
> ### Title: Observed-data log-posterior for normal data
> ### Aliases: logpost.norm
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> data(mdata)
> s <- prelim.norm(mdata)   #do preliminary manipulations
> prior <- list(0,.5,rep(0,ncol(mdata)),
+    .5*diag(rep(1,ncol(mdata)))) #ridge prior with .5 df
> thetahat <- em.norm(s,prior=prior)   #compute posterior mode
Iterations of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...
> logpost.norm(s,thetahat,prior)  #log-posterior at mode
[1] -12.51476
> 
> 
> 
> cleanEx(); ..nameEx <- "makeparam.norm"
> 
> ### * makeparam.norm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: makeparam.norm
> ### Title: Convert normal parameters to packed storage
> ### Aliases: makeparam.norm
> ### Keywords: regression
> 
> ### ** Examples
> 
> data(mdata)
> s <- prelim.norm(mdata)   #do preliminary manipulations
> thetahat <- em.norm(s)   #compute mle
Iterations of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...
> thetahat <- getparam.norm(s,thetahat,corr=TRUE)   #extract parameters
> thetahat$r   #look at mle correlations
            [,1]       [,2]        [,3]      [,4]       [,5]
[1,]  1.00000000 0.88753257 -0.09642342 0.4889176  0.7145953
[2,]  0.88753257 1.00000000  0.09345297 0.5197721  0.5429959
[3,] -0.09642342 0.09345297  1.00000000 0.2930083 -0.3598209
[4,]  0.48891759 0.51977207  0.29300826 1.0000000  0.3029066
[5,]  0.71459529 0.54299589 -0.35982094 0.3029066  1.0000000
> thetahat$r[1,2] <- .5   #tweak a parameter
> thetahat <- makeparam.norm(s,thetahat)   #convert to packed storage
> thetahat <- em.norm(s,thetahat) #run EM again from new starting value
Iterations of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...
> 
> 
> 
> cleanEx(); ..nameEx <- "mda.norm"
> 
> ### * mda.norm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mda.norm
> ### Title: Monotone data augmentation for incomplete multivariate normal
> ###   data
> ### Aliases: mda.norm
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> data(mdata)
> s <- prelim.norm(mdata)
> thetahat <- em.norm(s)   #find the MLE for a starting value
Iterations of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...
> rngseed(1234567)   #set random number generator seed
> theta <- mda.norm(s,thetahat,steps=20,showits=TRUE)  # take 20 steps
Steps of Monotone Data Augmentation: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...
> getparam.norm(s,theta) # look at result
$mu
[1]    42.241085    36.478323    13.558893 54054.686286     2.383784

$sigma
             [,1]         [,2]         [,3]          [,4]         [,5]
[1,] 1.065414e+02    56.499139    12.558988    135664.856     3.634417
[2,] 5.649914e+01    37.038244     8.996055     90189.082     1.494470
[3,] 1.255899e+01     8.996055    26.496029     61020.445    -2.598673
[4,] 1.356649e+05 90189.081764 61020.444609 653204090.398 -4411.142211
[5,] 3.634417e+00     1.494470    -2.598673     -4411.142     1.350514

> 
> 
> 
> cleanEx(); ..nameEx <- "ninvwish"
> 
> ### * ninvwish
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: ninvwish
> ### Title: Random normal-inverted Wishart variate
> ### Aliases: ninvwish
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> data(mdata)
> s <- prelim.norm(mdata)   #do preliminary manipulations
> params <- list(1,.5,rep(0,ncol(mdata)), .5*diag(rep(1,ncol(mdata)))) # gives widely dispersed values
> rngseed(1234567)
> start <- ninvwish(s,params)  # draw a variate
> thetahat <- em.norm(s,start=start) # run EM from this starting value
Iterations of EM: 
1...2...3...4...5...6...7...8...9...10...11...12...13...14...15...16...17...18...19...20...21...
> 
> 
> 
> cleanEx(); ..nameEx <- "prelim.norm"
> 
> ### * prelim.norm
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: prelim.norm
> ### Title: Preliminary manipulations for a matrix of incomplete continuous
> ###   data.
> ### Aliases: prelim.norm
> ### Keywords: multivariate
> 
> ### ** Examples
> 
> data(mdata)
> s <- prelim.norm(mdata)  #do preliminary manipulations 
> s$nmis[s$co] #look at nmis 
named numeric(0)
> s$r #look at missing data patterns
   [,1] [,2] [,3] [,4] [,5]
14    1    1    1    1    1
 4    1    0    1    1    1
 1    1    1    0    1    1
 1    1    0    0    1    1
 3    1    1    1    0    1
 1    1    1    0    0    1
 1    1    1    1    0    0
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  0.36 0.05 0.41 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')