Exercise 2.2 of The Elements of Statistical Learning The text book first generates some 2-class data via:


which gives:

and then it asks:

I try to solve this by first modelling this with this graphical model:

where $c$ is the label, $h\,(1\le h \le 10)$ is the index of the selected mean $m_h^c$, and $x$ is the data point. This will give
$$
\begin{align*}
\Pr(x\mid m_h^c) =& \mathcal{N}(m_h^c,\mathbf{I}/5)\\
\Pr(m_h^c\mid h,c=\mathrm{blue}) =& \mathcal{N}((1,0)^T,\mathbf{I})\\
\Pr(m_h^c\mid h,c=\mathrm{orange}) =& \mathcal{N}((0,1)^T,\mathbf{I})\\
\Pr(h) =& \frac{1}{10}\\
\Pr(c) =& \frac{1}{2}
\end{align*}
$$
On the other hand, the boundary is $\{x:\Pr(c=\mathrm{blue}\mid x)=\Pr(c=\mathrm{orange}\mid x)\}$. With Bayesian rule, we have
$$
\begin{align*}
\Pr(c\mid x) =& \frac{\Pr(x\mid c)\Pr(c)}{\sum_c\Pr(x\mid c)\Pr(c)}\\
\Pr(x\mid c) =& \sum_h\int_{m_h^c}\Pr(h)\Pr(m_h^c\mid h,c)\Pr(x\mid m_h^c)
\end{align*}
$$
But later I found that the problem setting is symmetric so this may yield $x=y$ as the boundary. If the problem is asking the boundary when $m_h^c$ are conditioned, the equation will include $40$ parameters which I think it's unlikely to be the purpose of the exercise.
So am I misunderstanding anything? Thank you.
 A: I don't think you're supposed to find an analytic expression for the Bayes decision boundary, for a given realisation of the $m_k$'s. Similarly I doubt you're supposed to get the boundary over the distribution of the $m_k$, since that's just $x=y$ by symmetry as you noted.
I think what you need is to show is a program that can compute the decision boundary for a given realisation of the $m_k$'s. This can be done by setting down a grid of $x$ and $y$ values, computing the class-conditional densities, and finding the points where they're equal.
This code is a stab at it. IIRC there's actually code to compute the decision boundary in Modern Applied Statistics with S, but I haven't got that handy right now.
# for dmvnorm/rmvnorm: multivariate normal distribution
library(mvtnorm)

# class-conditional density given mixture centers
f <- function(x, m)
{
    out <- numeric(nrow(x))
    for(i in seq_len(nrow(m)))
        out <- out + dmvnorm(x, m[i, ], diag(0.2, 2))
    out
}

# generate the class mixture centers
m1 <- rmvnorm(10, c(1,0), diag(2))
m2 <- rmvnorm(10, c(0,1), diag(2))
# and plot them
plot(m1, xlim=c(-2, 3), ylim=c(-2, 3), col="blue")
points(m2, col="red")

# display contours of the class-conditional densities
dens <- local({
    x <- y <- seq(-3, 4, len=701)
    f1 <- outer(x, y, function(x, y) f(cbind(x, y), m1))
    f2 <- outer(x, y, function(x, y) f(cbind(x, y), m2))
    list(x=x, y=y, f1=f1, f2=f2)
})

contour(dens$x, dens$y, dens$f1, col="lightblue", lty=2, levels=seq(.3, 3, len=10),
        labels="", add=TRUE)

contour(dens$x, dens$y, dens$f2, col="pink", lty=2, levels=seq(.3, 3, len=10),
        labels="", add=TRUE)

# find which points are on the Bayes decision boundary
eq <- local({
    f1 <- dens$f1
    f2 <- dens$f2
    pts <- seq(-3, 4, len=701)
    eq <- which(abs((dens$f1 - dens$f2)/(dens$f1 + dens$f2)) < 5e-3, arr.ind=TRUE)
    eq[,1] <- pts[eq[,1]]
    eq[,2] <- pts[eq[,2]]
    eq
})
points(eq, pch=16, cex=0.5, col="grey")


Result:


A: Actually, the book does ask to provide an analytical solution to this problem. And yes, you have to condition the boundary, but not on the 40 means: you never get to know them precisely. Instead you have to condition on the 200 data points that you do get to see. So you will need 200 parameters, but due to the use of summation, the answer doesn't look too complicated.
I would never be able to derive this formula, so I only take the credit for realizing that the analytical solution doesn't have to be ugly and then searching for it on google. Luckily, it's provided by the authors some nice people, pages 6-7.
A: Wish i stumbled upon code per above earlier ; just finsihed creating some alternative code per below... for what it is worth 
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")

