Optimal Decision surface with two classes.  Try moving the means in the plot.  Also see what happens as you give different variances to each class?  How about nonzero covariance?

Needs["MultivariateStatistics`"];
Manipulate[
r = 3;
μ={meanClass1,meanClass2};
Σ ={ (({{varClass1U, covClass1}, {covClass1, varClass1V}})),(({{varClass2U, covClass2}, {covClass2, varClass2V}}))};
DensityPlot[
Max[PDF[MultinormalDistribution[μ[[1]],Σ[[1]]],{u,v}],PDF[MultinormalDistribution[μ[[2]],Σ[[2]]],{u,v}]],
{u,-r,r},{v,-r,r}],
{{covClass1,0},-1,1}, {{varClass1U,1},0,3},
{{varClass1V,1},0,3},{{covClass2,0},-1,1}, 
{{varClass2U,1},0,3},{{varClass2V,1},0,3},
{{meanClass1,{-1,1}}, Locator},{{meanClass2,{1,-1}}, Locator}]

The optimal predicted probability of being in class 1

Manipulate[
r = 3;
μ={ma,mb};
Σ ={ (({{varAU, covA}, {covA, varAV}})),(({{varBU, covB}, {covB, varBV}}))};
DensityPlot[
PDF[MultinormalDistribution[μ[[1]],Σ[[1]]],{u,v}]  (((PDF[MultinormalDistribution[μ[[1]],Σ[[1]]],{u,v}] + PDF[MultinormalDistribution[μ[[2]],Σ[[2]]],{u,v}]))^(-1)),
{u,-r,r},{v,-r,r}],
{{covA,0},-1,1}, {{varAU,1},0,3},
{{varAV,1},0,3},{{covB,0},-1,1}, 
{{varBU,1},0,3},{{varBV,1},0,3},
{{ma,{-1,1}}, Locator},{{mb,{1,-1}}, Locator}]

Define a function to create artificial data

makeData = Function[{n,μ,Σ},
Table[
RandomReal[MultinormalDistribution[μ[[i]], Σ[[i]]],n],
{i,1,Length[μ]}
]
];

Generate and Plot some artificial data:

Clear[n,k,m,x,μ,σ];
n = 100;
m={{-1,1},{1,-1}};
k ={ (({{1.2, .3}, {.3, 1}})),(({{1, 0}, {0, .5}}))}; 
data = makeData[n,m,k];
ListPlot[data, PlotStyle->{{PointSize[Large], RGBColor[1,0,0]},{PointSize[Large],RGBColor[0,0,1]}}]

Define a function to implement Gaussian Naive Bayes

gaussianNB = Function[data,
{Map[Mean,data],Map[StandardDeviation, data]}
];

Compute the predicted mean and standard deviation:

{pμ, pσ} = gaussianNB[data];
Print["Class 1: μ=", MatrixForm[pμ[[1]]],", Σ=", MatrixForm[((DiagonalMatrix[pσ[[1]]])^2)]]
Print["Class 2: μ=", MatrixForm[pμ[[2]]],", Σ=", MatrixForm[((DiagonalMatrix[pσ[[2]]])^2)]]

We now define a function to compute P(Y|X)

pYgX=Function[{y, x, μ,σ},
pY = (1/2);
(((∏_(i=1)^(Length[x]))PDF[NormalDistribution[μ[[y,i]],σ[[y,i]]],x[[i]]]) pY)
((((∑_(j=1)^(Length[μ]))pY(∏_(i=1)^(Length[x]))PDF[NormalDistribution[μ[[j,i]],σ[[j,i]]],x[[i]]]) )^(-1))
];

Finally we can plot the predicted probability of being in class 1 for the entire space.

r=3;
Show[
DensityPlot[pYgX[1,{u,v}, pμ, pσ],{u,-r,r},{v,-r,r}],
ListPlot[data, PlotStyle->{{PointSize[Large], RGBColor[1,0,0]},{PointSize[Large],RGBColor[0,0,1]}}]
]