### Name: vbmp
### Title: Variational Bayesian Multinomial Probit Regression with Gaussian
###   Process Priors.
### Aliases: vbmp
### Keywords: models

### ** Examples


## -----------------------------------------------------------------------------
## EXAMPLE 1 - Theta estimate with synthetic data
## -----------------------------------------------------------------------------
## Samples of 2-D data points drawn from three nonlinearly separable
## classes which take the form of two annular rings and one zero-centered
## Gaussian are used in this little illustrative example. 
genSample <- function(n, noiseVar=0) {
    ## class 1 and 2 (x ~ U(0,1))
    u <- 4. * matrix(runif(2*n), nrow=n, ncol=2) - 2.;
    i <- which(((u[, 1]^2 + u[, 2]^2) > .1) & ((u[, 1]^2 + u[, 2]^2) < .5) );
    j <- which(((u[, 1]^2 + u[, 2]^2) > .6) & ((u[, 1]^2 + u[, 2]^2) < 1) );
    X <- u[c(i, j),];
    t.class <- c(rep(1, length(i)),rep(2, length(j)));
    ## class 3 (x ~ N(0,1))
    x <- 0.1 * matrix(rnorm(2*length(i)), ncol=2, nrow=length(i) );
    k <- which((x[, 1]^2 + x[, 2]^2) < 0.1);
    X <- rbind(X, x[k, ]);
    t.class <- c(t.class, rep(3, length(k)));
    ## add random coloumns
    if (noiseVar>0) X <- cbind(X, matrix(rnorm(noiseVar*nrow(X)), ncol=noiseVar, nrow=nrow(X)));
    structure( list( t.class=t.class, X=X), class="MultiNoisyData");
}

set.seed(100); ## Init random number generator

## Generate training and test samples as an independent 
## test set to assess out-of-sample prediction error 
## and predictive likelihoods.
nNoisyInputs <- 0;       ## number of additional noisy input parameters
Ntest <- Ntrain <- 500;  ## sample sizes
dataXt.train <- genSample(Ntrain, nNoisyInputs);
dataXt.test  <- genSample(Ntest,  nNoisyInputs);

## Not run: 
##D  
##D  theta <- runif(ncol(dataXt.train$X));
##D  res <- vbmp( dataXt.train$X, dataXt.train$t.class,
##D         dataXt.test$X, dataXt.test$t.class, theta, 
##D          control=list(bThetaEstimate = T, 
##D          bPlotFitting=T, maxIts=50));
## End(Not run)

## set theta params (previously estimated) 
theta <- c(0.09488309, 0.16141604);   
## Fit the vbmp
res <- vbmp( dataXt.train$X, dataXt.train$t.class,
        dataXt.test$X, dataXt.test$t.class, theta, 
        control=list(maxIts=5));
## print out-of-sample error estimate
predError(res);

## Not run: 
##D ## ----------------------------------------------------------
##D ## EXAMPLE 2 - BRCA12 genomic data
##D ## ----------------------------------------------------------
##D ## Leave-one-out (LOO) cross-validation prediction error of the probabilistic 
##D ## Gaussian process classifier used in Zsofia Kote-Jarai et al. 
##D ## Clin Cancer Res 2006;12(13);3896-3901
##D 
##D   if(any(installed.packages()[,1]=="Biobase")) {
##D     library("Biobase");
##D     data("BRCA12");
##D     brca.y <- BRCA12$Target.class;
##D     brca.x <- t(exprs(BRCA12));
##D   } else {
##D     print("Deprecated.....");
##D     load(url("http://www.dcs.gla.ac.uk/people/personal/girolami/pubs_2005/VBGP/BRCA12.RData"));
##D     brca.y <- as.numeric(BRCA12$y);
##D     brca.x <- as.matrix(BRCA12[,-1]);
##D   }
##D   
##D   sKernelType <- "iprod";  ## Covariance function type
##D   Thresh <- 1e-8;  ## Iteration threshold
##D   InfoLevel <- 1;
##D   theta <- rep(1.0, ncol(brca.x));
##D   ITER.THETA <- 24;
##D   n     <- nrow(brca.x) ;
##D   Kfold <- n; # number of folds , if equal to n then LOO
##D   samps <- sample(rep(1:Kfold, length=n), n, replace=FALSE); 
##D   res   <- rep(NA, n);
##D   print(paste("LOO crossvalidation started...... (",n,"steps)"));
##D   for (x in 1:Kfold) {
##D       cat(paste(x,", ",sep="")); flush.console();
##D       resX <- vbmp( brca.x[samps!=x,], brca.y[samps!=x], 
##D                     brca.x[samps==x,], brca.y[samps==x], 
##D                     theta,  control=list(bThetaEstimate=F, 
##D                     bPlotFitting=F, maxIts=ITER.THETA, 
##D                     sKernelType=sKernelType, Thresh=Thresh));    
##D       res[samps==x] <- predClass(resX); 
##D   }
##D   print("(end)");
##D   print(paste("Crossvalidated error rate", round(sum(res!=brca.y)/n,2)));
## End(Not run)   




