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("MNP-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('MNP')
Loading required package: MASS
MNP: R Package for Fitting the Multinomial Probit Model
Version2.3-3
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "mnp"
> 
> ### * mnp
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: mnp
> ### Title: Fitting the Multinomial Probit Model via Markov chain Monte
> ###   Carlo
> ### Aliases: mnp MNP
> ### Keywords: models
> 
> ### ** Examples
> 
> ###
> ### NOTE: this example is not fully analyzed. In particular, the
> ### convergence has not been assessed. A full analysis of these data
> ### sets appear in Imai and van Dyk (2005b).
> ###
> 
> ## load the detergent data
> data(detergent)
> ## run the standard multinomial probit model with intercepts and the price
> res1 <- mnp(choice ~ 1, choiceX = list(Surf=SurfPrice, Tide=TidePrice,
+                                        Wisk=WiskPrice, EraPlus=EraPlusPrice,
+                                        Solo=SoloPrice, All=AllPrice),
+             cXnames = "price", data = detergent, n.draws = 500, burnin = 100,
+             thin = 3, verbose = TRUE)

The base category is `All'.

The total number of alternatives is 6.

The choice-specific variables of the base category are subtracted from the corresponding variables of the non-base categories.

The dimension of beta is 6.

The number of observations is 2657.

Improper prior will be used for beta.

Starting Gibbs sampler...
 10 percent done.
 20 percent done.
 30 percent done.
 40 percent done.
 50 percent done.
 60 percent done.
 70 percent done.
 80 percent done.
 90 percent done.
100 percent done.
> ## summarize the results
> summary(res1)

Call:
mnp(formula = choice ~ 1, data = detergent, choiceX = list(Surf = SurfPrice, 
    Tide = TidePrice, Wisk = WiskPrice, EraPlus = EraPlusPrice, 
    Solo = SoloPrice, All = AllPrice), cXnames = "price", n.draws = 500, 
    burnin = 100, thin = 3, verbose = TRUE)


Coefficients:
                        mean std.dev.     2.5%   97.5%
(Intercept):EraPlus   2.1966   0.2192   1.8160   2.545
(Intercept):Solo      1.3831   0.2242   0.9208   1.706
(Intercept):Surf      1.2647   0.1012   1.0387   1.446
(Intercept):Tide      2.3833   0.1845   2.0634   2.677
(Intercept):Wisk      1.3192   0.1279   1.1138   1.565
price               -79.1934   6.9255 -88.1821 -65.282

Covariances:
                    mean std.dev.     2.5% 97.5%
EraPlus:EraPlus  1.00000  0.00000  1.00000 1.000
EraPlus:Solo     0.34087  0.14374  0.06304 0.587
EraPlus:Surf    -0.02658  0.12259 -0.22693 0.205
EraPlus:Tide     0.12939  0.09563 -0.02110 0.290
EraPlus:Wisk     0.62405  0.16749  0.30813 0.907
Solo:Solo        1.55901  0.32987  0.93771 2.181
Solo:Surf       -0.08547  0.15307 -0.30159 0.236
Solo:Tide       -0.20211  0.18039 -0.50049 0.076
Solo:Wisk        0.69290  0.15359  0.40745 1.045
Surf:Surf        1.31737  0.29845  0.87349 1.777
Surf:Tide        0.54326  0.16867  0.32564 0.896
Surf:Wisk        0.52049  0.17661  0.14852 0.811
Tide:Tide        1.02872  0.24572  0.64421 1.439
Tide:Wisk        0.66747  0.23213  0.28968 1.094
Wisk:Wisk        1.94255  0.37056  1.22170 2.511

Base category: All
Number of alternatives: 6
Number of observations: 2657
Number of stored MCMC draws: 100

> ## calculate the predicted probabilities for the first 5 observations
> predict(res1, newdata = detergent[1:3,], type="prob", verbose = TRUE)
There are 3 observations to predict. Please wait...
 10 percent done.
 20 percent done.
 30 percent done.
 40 percent done.
 50 percent done.
 60 percent done.
 70 percent done.
 80 percent done.
 90 percent done.
 100 percent done.
$p
      All EraPlus Solo Surf Tide Wisk
[1,] 0.01    0.25 0.12 0.46 0.11 0.05
[2,] 0.04    0.13 0.04 0.03 0.40 0.36
[3,] 0.03    0.14 0.15 0.00 0.36 0.32

> 
> ## load the Japanese election data
> data(japan)
> ## run the multinomial probit model with ordered preferences
> res2 <- mnp(cbind(LDP, NFP, SKG, JCP) ~ gender + education + age, data = japan,
+             verbose = TRUE)

The base category is `JCP'.

The total number of alternatives is 4.

The dimension of beta is 12.

The number of observations is 418.

Improper prior will be used for beta.

Starting Gibbs sampler...
 10 percent done.
 20 percent done.
 30 percent done.
 40 percent done.
 50 percent done.
 60 percent done.
 70 percent done.
 80 percent done.
 90 percent done.
100 percent done.
> ## summarize the results
> summary(res2)

Call:
mnp(formula = cbind(LDP, NFP, SKG, JCP) ~ gender + education + 
    age, data = japan, verbose = TRUE)


Coefficients:
                      mean   std.dev.       2.5% 97.5%
(Intercept):LDP  0.8377753  0.4260994  0.0210850 1.674
(Intercept):NFP  1.1546068  0.4736017  0.2117028 2.085
(Intercept):SKG  0.4269408  0.3770675 -0.3218562 1.149
gendermale:LDP  -0.1105414  0.1601280 -0.4296236 0.201
gendermale:NFP  -0.2268267  0.1745650 -0.5792927 0.102
gendermale:SKG  -0.1401373  0.1462414 -0.4315418 0.147
education:LDP   -0.1034982  0.0762202 -0.2531932 0.046
education:NFP   -0.1046450  0.0857684 -0.2731799 0.060
education:SKG   -0.0014338  0.0670919 -0.1365451 0.132
age:LDP          0.0132518  0.0065138  0.0002626 0.026
age:NFP          0.0065625  0.0070828 -0.0072325 0.020
age:SKG          0.0093534  0.0057698 -0.0019475 0.021

Covariances:
           mean std.dev.    2.5% 97.5%
LDP:LDP 1.00000  0.00000 1.00000 1.000
LDP:NFP 1.05835  0.06166 0.94411 1.169
LDP:SKG 0.71878  0.05978 0.60000 0.830
NFP:NFP 1.43448  0.14022 1.18697 1.739
NFP:SKG 0.76435  0.08359 0.60458 0.927
SKG:SKG 0.71507  0.08839 0.55776 0.892

Base category: JCP
Number of alternatives: 4
Number of observations: 418
Number of stored MCMC draws: 5000

> ## calculate the predicted probabilities for the 10th observation
> predict(res2, newdata = japan[10,], type = "prob")
$p
       JCP   LDP    NFP    SKG
[1,] 0.105 0.357 0.3308 0.2072

> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  41.21 0.21 41.78 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')