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("ensembleBMA-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('ensembleBMA')
> 
> assign(".oldSearch", search(), env = .CheckExEnv)
> assign(".oldNS", loadedNamespaces(), env = .CheckExEnv)
> cleanEx(); ..nameEx <- "CRPS"
> 
> ### * CRPS
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: CRPS
> ### Title: Calculate the Continuous Ranked Probability Score for an
> ###   ensemble forecast.
> ### Aliases: CRPS
> ### Keywords: file
> 
> ### ** Examples
> 
> #read in the sea-level pressure data and calculate BMA estimates
> #for forecasting on the 35th day in the data set
> data(slp)
> unique.dates <- unique(slp$date)
> date.list <- NULL
> 
> for(i in 1:length(unique.dates))
+ {
+   date.list[slp$date==unique.dates[i]] <- i
+ }
> 
> X <- cbind(slp$F1,slp$F2,slp$F3,slp$F4,slp$F5)
> Y <- slp$Y
> 
> EMresult <- EM.for.date(date = 35,date.list = date.list,X = X,Y = Y )
> 
> #now calculate the CRPS over the training period (observations 43 through 161)
> CRPS(a = EMresult$a,b = EMresult$b, sigma = EMresult$sigma, w =  EMresult$w, X=X[43:161,], Y=Y[43:161])
[1] 2.081511
> 
> #calculate the BMA estimates without CRPS minimization, and compare the new CRPS score
> EMresult.without <- EM.for.date(date = 35,date.list = date.list,X = X,Y = Y, min.CRPS=FALSE )
> CRPS(a = EMresult.without$a,b = EMresult.without$b, sigma = EMresult.without$sigma, w =  EMresult.without$w, X=X[43:161,], Y=Y[43:161])
[1] 2.081943
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "EM.for.date"
> 
> ### * EM.for.date
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: EM.for.date
> ### Title: EM.normals wrapper
> ### Aliases: EM.for.date
> ### Keywords: file
> 
> ### ** Examples
> 
> #read in the sea-level pressure data and calculate BMA estimates
> #for forecasting on the 35th day in the data set
> data(slp)
> unique.dates <- unique(slp$date)
> date.list <- NULL
> 
> for(i in 1:length(unique.dates))
+ {
+   date.list[slp$date==unique.dates[i]] <- i
+ }
> 
> X <- cbind(slp$F1,slp$F2,slp$F3,slp$F4,slp$F5)
> Y <- slp$Y
> 
> EMresult <- EM.for.date(date = 35,date.list = date.list,X = X,Y = Y )
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "EM.normals"
> 
> ### * EM.normals
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: EM.normals
> ### Title: EM algorithm for a BMA mixture of normals
> ### Aliases: EM.normals
> ### Keywords: file
> 
> ### ** Examples
> 
> 
> #create a simulated dataset with equal weights, no bias,
> #and standard deviation of 1 in each component
> x <- matrix(rnorm(1000,0,2),nrow = 200, ncol = 5)
> 
> y.latent <- floor(runif(200,1,6))
> y.means <- NULL
> for(i in 1:200)
+ {
+   y.means[i] <- x[i,y.latent[i]]
+ }
> y <- rnorm(200,y.means, sd = 1)
> 
> #calculate the BMA estimates of the parameters
> EMresult <- EM.normals(x, y, reg.adjust=FALSE, min.CRPS=FALSE)
> 
> #read in the sea-level pressure data and calculate BMA estimates
> #for forecasting on the 35th day in the data set
> #(this requires training on observations 43 through 161)
> data(slp)
> unique.dates <- unique(slp$date)
> date.list <- NULL
> 
> for(i in 1:length(unique.dates))
+ {
+   date.list[slp$date==unique.dates[i]] <- i
+ }
> 
> X <- cbind(slp$F1,slp$F2,slp$F3,slp$F4,slp$F5)
> Y <- slp$Y
> 
> #calculate the BMA estimates of the parameters
> EMresult <- EM.normals(X = X[43:161,],Y = Y[43:161] )
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "bmacdf"
> 
> ### * bmacdf
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bmacdf
> ### Title: Find the BMA cdf at x for the K-component BMA mixture of normals
> ### Aliases: bmacdf
> ### Keywords: file
> 
> ### ** Examples
> 
> #create a simulated dataset with equal weights, no bias,
> #and standard deviation of 1 in each component
> x <- matrix(rnorm(1000,0,2),nrow = 200, ncol = 5)
> 
> y.latent <- floor(runif(200,1,6))
> y.means <- NULL
> for(i in 1:200)
+ {
+   y.means[i] <- x[i,y.latent[i]]
+ }
> y <- rnorm(200,y.means, sd = 1)
> 
> #calculate the BMA estimates of the parameters
> EMresult <- EM.normals(x, y, reg.adjust=FALSE, min.CRPS=FALSE)
> 
> #evaluate and plot the BMA cdf for the first observation
> index <- seq(-5,5,by=.1)
> cdf.at.index=NULL
> for(i in 1:length(index))
+ {
+   cdf.at.index[i]=bmacdf(a = EMresult$a,b = EMresult$b, sigma = EMresult$sigma,
+   w =  EMresult$w, x[1,],index[i])
+ }
> plot(index, cdf.at.index, type="l")
> 
> #read in the sea-level pressure data and calculate BMA estimates
> #for forecasting on the 35th day in the data set
> data(slp)
> unique.dates <- unique(slp$date)
> date.list <- NULL
> 
> for(i in 1:length(unique.dates))
+ {
+   date.list[slp$date==unique.dates[i]] <- i
+ }
> 
> X <- cbind(slp$F1,slp$F2,slp$F3,slp$F4,slp$F5)
> Y <- slp$Y
> 
> EMresult <- EM.for.date(date = 35,date.list = date.list,X = X,Y = Y )
> 
> #evaluate and plot the BMA cdf for the first observation
> index <- seq(1000,1025,by=.1)
> cdf.at.index=NULL
> for(i in 1:length(index))
+ {
+   cdf.at.index[i]=bmacdf(a = EMresult$a,b = EMresult$b, sigma = EMresult$sigma,
+   w =  EMresult$w, X[1,],index[i])
+ }
> plot(index, cdf.at.index, type="l")
> 
> 
> 
> 
> cleanEx(); ..nameEx <- "bmaquant"
> 
> ### * bmaquant
> 
> flush(stderr()); flush(stdout())
> 
> ### Name: bmaquant
> ### Title: Find a specific quantile of a BMA mixture of normal
> ###   distributions
> ### Aliases: bmaquant
> ### Keywords: file
> 
> ### ** Examples
> 
> #create a simulated dataset with equal weights, no bias,
> #and standard deviation of 1 in each component
> x <- matrix(rnorm(1000,0,2),nrow = 200, ncol = 5)
> 
> y.latent <- floor(runif(200,1,6))
> y.means <- NULL
> for(i in 1:200)
+ {
+   y.means[i] <- x[i,y.latent[i]]
+ }
> y <- rnorm(200,y.means, sd = 1)
> 
> #calculate the BMA estimates of the parameters
> EMresult <- EM.normals(x, y, reg.adjust=FALSE, min.CRPS=FALSE)
> 
> # 95th percentile
> bmaquant(a = EMresult$a,b = EMresult$b, sigma = EMresult$sigma,
+ w =  EMresult$w, alpha = .95, x[1,])
[1] 3.103819
> 
> # 5th percentile
> bmaquant(a = EMresult$a,b = EMresult$b, sigma = EMresult$sigma,
+ w =  EMresult$w, alpha = .05 ,x[1,])
[1] -3.126049
> 
> #read in the sea-level pressure data and calculate BMA estimates
> #for forecasting on the 35th day in the data set
> data(slp)
> unique.dates <- unique(slp$date)
> date.list <- NULL
> 
> for(i in 1:length(unique.dates))
+ {
+   date.list[slp$date==unique.dates[i]] <- i
+ }
> 
> X <- cbind(slp$F1,slp$F2,slp$F3,slp$F4,slp$F5)
> Y <- slp$Y
> 
> EMresult <- EM.for.date(date = 35,date.list = date.list,X = X,Y = Y )
> 
> # 5th and 95th percentiles
> bmaquant(a = EMresult$a,b = EMresult$b, sigma = EMresult$sigma,
+ w =  EMresult$w, alpha = c(.05,.95) ,X[1,])
[1] 1007.636 1018.190
> 
> 
> 
> 
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", env = .CheckExEnv),"\n")
Time elapsed:  16.95 0.06 17.06 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')