I was struggling with the same question; created some code that might be helpful. Didnt manage to create a good function to add contour lines, but per below should show the principle.
set.seed(1)
library(MASS)
#create original 10 center points/means for each class
I.mat=diag(2)
mu1=c(1,0);mu2=c(0,1)
mv.dist1=mvrnorm(n = 10, mu1, I.mat)
mv.dist2=mvrnorm(n = 10, mu2, I.mat)
values1=NULL;values2=NULL
#create 100 observations for each class, after random sampling of a center point, based on an assumed bivariate probability distribution around each center point
for(i in 1:10){
mv.values1=mv.dist1[sample(nrow(mv.dist1),size=1,replace=TRUE),]
sub.mv.dist1=mvrnorm(n = 10, mv.values1, I.mat/5)
values1=rbind(sub.mv.dist1,values1)
}
values1
#similar as per above, for second class
for(i in 1:10){
mv.values2=mv.dist2[sample(nrow(mv.dist2),size=1,replace=TRUE),]
sub.mv.dist2=mvrnorm(n = 10, mv.values2, I.mat/5)
values2=rbind(sub.mv.dist2,values2)
}
values2
#did not find probability function in MASS, so used mnormt
library(mnormt)
#create grid of points
grid.vector1=seq(-2,2,0.1)
grid.vector2=seq(-2,2,0.1)
length(grid.vector1)*length(grid.vector2)
grid=expand.grid(grid.vector1,grid.vector2)
#calculate density for each point on grid for each of the 100 multivariates distributions
prob.1=matrix(0:0,nrow=1681,ncol=10) #initialize grid
for (i in 1:1681){
for (j in 1:10){
prob.1[i,j]=dmnorm(grid[i,], mv.dist1[j,], I.mat/5)
}
}
prob.1
prob1.max=apply(prob.1,1,max)
#second class - as per above
prob.2=matrix(0:0,nrow=1681,ncol=10) #initialize grid
for (i in 1:1681){
for (j in 1:10){
prob.2[i,j]=dmnorm(grid[i,], mv.dist2[j,], I.mat/5)
}
}
prob.2
prob2.max=apply(prob.2,1,max)
#bind
prob.total=cbind(prob1.max,prob2.max)
class=rep(1,1681)
class[prob1.max<prob2.max]=2
cbind(prob.total,class)
#plot points
plot(grid[,1], grid[,2],pch=".", cex=3,col=ifelse(class==1, "coral", "cornflowerblue"))
points(values1,col="coral")
points(values2,col="cornflowerblue")
#check - original centers
# points(mv.dist1,col="coral")
# points(mv.dist2,col="cornflowerblue")
# bayesian decision boundary; first get conditional probability for class 1
prob.bayes <- matrix(prob1.max/(prob1.max + prob2.max), length(grid.vector1), length(grid.vector2))
contour(grid.vector1, grid.vector2, prob.bayes, levels = 0.5, labels = "", lwd = 2, add = TRUE)