# 02.04.2008, Thomas Stibor library(mvtnorm); # to install mvtnorm library # start R and type install.packages("mvtnorm"); ############################################################ data.gaussian <- function(N,mean.C1,cov.C1,mean.C2,cov.C2) { ############################################################ x <- rmvnorm(N,mean=mean.C1,sigma=cov.C1); y <- rmvnorm(N,mean=mean.C2,sigma=cov.C2); typ <- c(rep("red",N),rep("blue",N)); return(data.frame(rbind(x,y),typ)); } ############################################## gaussian.density <- function(x,mean.C,cov.C) { ############################################## return(dmvnorm(x,mean=mean.C,sigma=cov.C)); } ########################## g.sigmoid <- function(a) { ########################## return (1/(1+exp(-a))); } ############################################### decision.line <- function(w,x.min,x.max) { ############################################### x <- c(x.min,x.max); y <- -w[2]/w[3] * x - w[1]/w[3]; return (matrix(c(x,y),nrow=2,ncol=2)); } ############################################### # main ############################################### N <- 200; mean.C1 <- c(3,3); cov.C1 <- matrix(c(1,0,0,1),ncol=2,nrow=2); mean.C2 <- c(1,1); cov.C2 <- cov.C1; data <- data.gaussian(N,mean.C1,cov.C1,mean.C2,cov.C2); names(data)[1] <- "x"; names(data)[2] <- "y"; # inverse covariance matrix cov.inv <- solve(cov.C1); # calculate weights w <- cov.inv %*% (mean.C1 - mean.C2); w.0 <- -1/2 * t(mean.C1) %*% cov.inv %*% mean.C1 + 1/2 * t(mean.C2) %*% cov.inv %*% mean.C2; # + log((3/4)/(1/4)); # priors are equal => term log(P(C1)/P(C2)) = 0, having distinct priors resulting in parallel shifts of the decision boundary # for observation add eg. log((3/4)/(1/4)) to w.0; # determine min-max range rng <- apply(cbind(data$x,data$y),2,range); # setup parameters for postscript output postscript("logistic.eps", width = 10.0, height = 10.0, horizontal = FALSE, onefile = FALSE, paper = "special"); # setting up coord. system plot(data$x, data$y, type = "n", asp=0, xlab = "x", ylab = "y", xlim=c(rng[1,1],rng[2,1]), ylim=c(rng[1,2],rng[2,2]), cex.axis = 1.2, axes=TRUE); # plot points points(data$x,data$y,col=as.vector(data$typ),pch=19,cex=2); # plot decision boundary lines(decision.line(c(w.0,w),rng[1,1],rng[2,1]),type="l",lwd=4,col="black"); # create tx times ty grid resol <- 200; tx <- seq(rng[1,1],rng[2,1],length=resol); ty <- seq(rng[1,2],rng[2,2],length=resol); pnts <- matrix(nrow=length(tx)*length(ty),ncol=2); k <- 1 for(j in 1:length(ty)){ for(i in 1:length(tx)){ pnts[k,] <- c(tx[i],ty[j]); k <- k+1; } } # plot density contour for both classes z.C1 <- matrix(gaussian.density(pnts,mean.C1,cov.C1),nrow=length(tx),ncol=length(ty)); contour(tx, ty, z.C1, add = TRUE, col="red"); z.C2 <- matrix(gaussian.density(pnts,mean.C2,cov.C2),nrow=length(tx),ncol=length(ty)); contour(tx, ty, z.C2, add = TRUE, col="blue"); dev.off();