\HeaderA{vbmp}{Variational Bayesian Multinomial Probit Regression with Gaussian Process Priors.}{vbmp}
\keyword{models}{vbmp}
\begin{Description}\relax
Used to fit a Multinomial Probit Regression model, specified by giving the 
matrix design \code{X}, the associated response variables \code{t.class}, kernel type and covariate 
scaling parameters. Covariance paramters can be inferred from the data.
\end{Description}
\begin{Usage}
\begin{verbatim}
vbmp(X, t.class, X.TEST, t.class.TEST, theta, control = list())
\end{verbatim}
\end{Usage}
\begin{Arguments}
\begin{ldescription}
\item[\code{X}] Feature matrix for parameter 'estimation'
\item[\code{t.class}] Target values, integer number used for class labels. 
\item[\code{X.TEST}] Feature matrix to compute out-of-sample (test) prediction errors and likelihoods 
\item[\code{t.class.TEST}] Target values for test data 
\item[\code{theta}] The covariance function parameters (e.g. scaling coefficients for each dimension) 
\item[\code{control}] A list of control parameters. See Details
\end{ldescription}
\end{Arguments}
\begin{Details}\relax
In this implementation a single covariance function is shared across all classes.
Compute the predictive posteriors on the test set and the associated likelihood 
and test errors at each iteration.

The control argument is a list that can supply any of the following components:
\item[InfoLevel] 0 to suppress tracing ( > 0  to print different levels
of monitoring information) 
\item[sFILE.TRACE] File name where to redirect output (default NULL) 
\item[bThetaEstimate] if covariance parameter estimation switched on. Defaults to FALSE (switched off) 
\item[sKernelType] Kernel function used in training and predicting. 
Currently implemented kernels are Gaussian (\code{"gauss"}), 
Cauchy (\code{"cauchy"}), Laplace (\code{"laplace"}),
Polynomial (\code{"poly"}), Homogeneous polynomial (\code{"hpoly"}),
'Thin-plate' spline (\code{"tps"}), 'linear' spline (\code{"lsp"}) and
Inner product(\code{"iprod"}).
Defaults to \code{"gauss"}. 
\item[maxIts] Maximum number of variational EM steps to take. 
Defaults to 50.
\item[Thresh] Convergence threshold on marginal likelihood lowerbound. 
Defaults to 1e-4.
\item[method] Integral computation method: "quadrature" (Gaussian quadrature)
or "classic"(simple sampler).  Defaults to "quadrature".     
\item[nNodesQuad] Number of nodes used for quadrature. Defaults to 49.     
\item[nSampsTG] Number of samples used in obtaining mean of truncated 
Gaussian. Defaults to 1000.
\item[nSampsIS] Number of samples used in the importance sampler. 
Defaults to 1000.
\item[nSmallNo] Small number used to prevent numerical problems 
(ill-conditioned covariance matrix). Defaults to 1e-10.
\item[parGammaTau,parGammaSigma] The location and scale parameters of the 
Gamma prior over covaraince params. Default to 1e-6.
\item[bMonitor] TRUE to collect monitor convergence diagnostics at each
iteration. Defaults to FALSE.
\item[bPlotFitting] TRUE to plot test performance results at each iteration 
during model estimation (if TRUE it forces bMonitor to TRUE). Defaults to FALSE.
\end{Details}
\begin{Value}
\code{vbmp} returns an object of class "VBMP.obj". 
An object of class "VBMP.obj" is a list containing at least the following components: 
\begin{ldescription}
\item[\code{Kc}] Number of classes
\item[\code{Ptest}] Matrix of multinomial class predictive posterior probabilities for the test data 
\item[\code{X}] Feature matrix
\item[\code{invPHI}] Inverse of the Kernel matrix
\item[\code{Y}] Matrix of auxiliary variables
\item[\code{M}] Matrix of GP random variables
\item[\code{theta}] covariance kernel hyperparameters (estimates computed during 
model fitting, if inferred
\item[\code{sKernelType }] Kernel function used in training and predicting
\item[\code{Test.Err}] Out-of-Sample Percent Prediction error estimates computed 
during model fitting (0-1 error loss).
\item[\code{PL}] Predictive Likelihood estimates computed during model fitting
\item[\code{LOWER.BOUND}] Lower bound estimates computed during model fitting
\end{ldescription}
\end{Value}
\begin{Author}\relax
N Lama \email{nicola.lama@unina2.it}, MA Girolami \email{girolami@dcs.gla.ac.uk}
\end{Author}
\begin{References}\relax
Girolami M, Rogers S, \emph{Variational Bayesian Multinomial Probit Regression with Gaussian Process Priors}, Neural Computation 18, 1790-1817 (2006).
Lama N, Girolami M \emph{vbmp: Variational Bayesian Multinomial Probit Regression for multi-class classification in R}, Bioinformatics 24(1):135-136 (2008). 
\url{http://bioinformatics.oxfordjournals.org/cgi/content/short/btm535v1}
\end{References}
\begin{SeeAlso}\relax
See Also as \code{\LinkA{predictCPP}{predictCPP}}, \code{\LinkA{covParams}{covParams}}, 
\code{\LinkA{lowerBound}{lowerBound}}, \code{\LinkA{predError}{predError}}, 
\code{\LinkA{predLik}{predLik}}, \code{\LinkA{predClass}{predClass}}
\end{SeeAlso}
\begin{Examples}
\begin{ExampleCode}

## -----------------------------------------------------------------------------
## 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(123); ## 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: 
 
 theta <- runif(ncol(dataXt.train$X));
 res <- vbmp( dataXt.train$X, dataXt.train$t.class,
        dataXt.test$X, dataXt.test$t.class, theta, 
         control=list(bThetaEstimate = T, 
         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: 
## ----------------------------------------------------------
## EXAMPLE 2 - BRCA12 genomic data
## ----------------------------------------------------------
## Leave-one-out (LOO) cross-validation prediction error of the probabilistic 
## Gaussian process classifier used in Zsofia Kote-Jarai et al. 
## Clin Cancer Res 2006;12(13);3896-3901

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

\end{ExampleCode}
\end{Examples}

