inst/doc/vbmp.R

In vbmp: Variational Bayesian Multinomial Probit Regression

### R code from vignette source 'vbmp.Rnw'

###################################################
### code chunk number 1: example1
###################################################
## Init random number generator
set.seed(123); 


###################################################
### code chunk number 2: example11
###################################################
genSample <- function(n, noiseVar=0) {
    KF <- 10; # magnify factor
    nt <- round(n/3,0)#
    ## class 1 and 2 (x ~ U(0,1))
    u <- 4. * matrix(runif(2*KF*n), nrow=KF*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(sample(i,nt), sample(j,nt)),];
    t.class <- c(rep(1, nt),rep(2, nt));
    ## class 3 (x ~ N(0,1))
    x <- 0.1 * matrix(rnorm(2*KF*length(i)), ncol=2, nrow=length(i)*KF );
    k <- which((x[, 1]^2 + x[, 2]^2) < 0.1);
    nt <- n - 2*nt;
    X <- rbind(X, x[sample(k,nt), ]);
    t.class <- c(t.class, rep(3, nt));
    ## limit number
    #n <- min(n, nrow(X)); i <- sample(1:nrow(X),n);X <- X[i,]; t.class <- t.class[i];
    ## add random coloumns
    if (noiseVar>0) X <- cbind(X, matrix(rnorm(noiseVar*n), ncol=noiseVar, nrow=nrow(X)));
    structure( list( t.class=t.class, X=X), class="MultiNoisyData");
}


###################################################
### code chunk number 3: example12
###################################################
nNoisyInputs <- 4;       


###################################################
### code chunk number 4: example13
###################################################
Ntrain <- 100;  
Ntest  <- Ntrain * 5;


###################################################
### code chunk number 5: example14
###################################################
dataXtTrain <- genSample(Ntrain, nNoisyInputs);
dataXtTest  <- genSample(Ntest,  nNoisyInputs);


###################################################
### code chunk number 6: example15
###################################################
theta <- rep(1., ncol(dataXtTrain$X));
max.train.iter <- 12;


###################################################
### code chunk number 7: example16
###################################################
library(vbmp);
res <- vbmp( dataXtTrain$X, dataXtTrain$t.class,
             dataXtTest$X, dataXtTest$t.class, theta, 
             control=list(bThetaEstimate=T, bMonitor=T, maxIts=max.train.iter)); 


###################################################
### code chunk number 8: example17
###################################################
predError(res);   


###################################################
### code chunk number 9: vbmp.Rnw:130-137
###################################################
    i.t1 <- which(dataXtTest$t.class == 1);
    i.t2 <- which(dataXtTest$t.class == 2);
    i.t3 <- which(dataXtTest$t.class == 3);
    plot(dataXtTest$X[, 1], dataXtTest$X[,2], type="n", xlab="X1", ylab="X2");
    points(dataXtTest$X[i.t1, 1],dataXtTest$X[i.t1, 2], type="p", col="blue");
    points(dataXtTest$X[i.t2, 1],dataXtTest$X[i.t2, 2], type="p", col="red");
    points(dataXtTest$X[i.t3, 1],dataXtTest$X[i.t3, 2], type="p", col="green");


###################################################
### code chunk number 10: vbmp.Rnw:146-148
###################################################
## plot convergence diagnostics (same as setting bPlotFitting=T)
plotDiagnostics(res);       


###################################################
### code chunk number 11: example2
###################################################
  library("Biobase");
  data("BRCA12");
  brca.y <- BRCA12$Target.class;
  brca.x <- t(exprs(BRCA12));


###################################################
### code chunk number 12: example21
###################################################
  predictVBMP <- function(x) {
      sKernelType <- "iprod";  
      Thresh      <- 1e-8; 
      theta       <- rep(1.0, ncol(brca.x));
      max.train.iter <- 24;
      resX <- vbmp( brca.x[!x,], brca.y[!x], 
                    brca.x[ x,], brca.y[ x], 
                    theta,  control=list(bThetaEstimate=F, 
                    bPlotFitting=F, maxIts=max.train.iter, 
                    sKernelType=sKernelType, Thresh=Thresh));      
      predClass(resX);
  }


###################################################
### code chunk number 13: example212
###################################################
  # fake predictVBMP to speed-up the RCMD TEST
  predictVBMP <- function(x) {
      rep(0,sum(x));
  }


###################################################
### code chunk number 14: example22
###################################################
  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("Crossvalidation started...... (",n,"steps )"));
  for (x in 1:Kfold) {
      cat(paste(x,", ",ifelse(x%%10==0,"\n",""),sep="")); 
      flush.console();
      res[samps==x] <- predictVBMP(samps==x); 
  }


###################################################
### code chunk number 15: example23
###################################################
    CVerrorRate <- round(sum(res!=brca.y)/n,2);
  # don't print out the results, owing to the fake predictVBMP