The objective of Bayesian classification is to determine, given an observation $x = (x_1, \dots, x_n)$, the probability for $x$ to belong to classes $C_1 \dots C_K$. These are given by $p(C_k | x_1 \dots x_n)$
Using the Bayes theorem, we can rewrite the previous formula using $$ p(C_k | x) = \frac{p(C_k) p(x | C_k)}{p(x)} $$ The first thing to observe in that new formulation is that the quantity $\underset{1 \leq k \leq K}{\operatorname{argmax}} \; p(C_k | x)$ does not depend on $p(x)$ so this quantity can simply be ignored. Now using multiple times the definition of conditional probability, $$ p(x | C_k) = \prod_{i = 1}^n p(x_i | C_k, x_1, \dots, x_{i-1}) $$ Assuming that each parameter $x_i$ is independent, $p(x_i | C_k, x_j) = p(x_i | C_k)$. The previous equation can simply be rewritten $$ P(C_k | x) = \frac{p(C_k)}{p(x)} \prod_{i=1}^n p(x_i | C_k) $$ This quantity is much simpler to estimate. Let's say that each $x_i$ can take up to $v$ different values, then there exist around $v^n$ combinations for $x$. To estimate $p(x | C_k)$ we must have in our dataset at least $K \cdot v^n$ entries. And if we want a correct estimation we should get at least 1000 times that amount. Let's say that our dataset has $n=15$, $v=1000$ and $K=5$, then the number of entries must already be bigger than the number of atoms in the Moon. The same simple estimation becomes $K \cdot n \cdot v$ which is much more reasonable to estimate.
Finally, $$ \begin{align*} \underset{1 \leq k \leq K}{\operatorname{argmax}} \, p(C_k | x) &= \underset{1 \leq k \leq K}{\operatorname{argmax}} \, p(C_k) \prod_{i=1}^n p(x_i | C_k) \\ &= \underset{1 \leq k \leq K}{\operatorname{argmax}} \, \left[ \log p(C_k) + \sum_{i=1}^n \log p(x_i | C_k) \right] \end{align*} $$
The next step is to find a way to estimate $p(C_k)$ and $p(x_i | C_k)$. Of course this depends on whether $x_i$ is continuous or not.
We note $N_k$ the number of entries in class $C_k$ such that $\sum N_k = N$ and $N_{i, k, v}$ the number of entries $x \in C_k$ such that $x_i = v$. Then we can estimate
$$ \begin{align*} \underset{1 \leq k \leq K}{\operatorname{argmax}} \, p(C_k | x) &\approx \underset{1 \leq k \leq K}{\operatorname{argmax}} \, \left[ \log \frac{N_k}{N} + \sum_{i=1}^n \log \frac{N_{i, k, x_i}}{N_k} \right] \\ &= \underset{1 \leq k \leq K}{\operatorname{argmax}} \, \left[ \log \frac{N_k}{N} - n \log N_k + \sum_{i=1}^n \log N_{i, k, x_i} \right] \\ \end{align*} $$But there is an issue with that estimation: let's say that for $x$ a new entry not in the dataset and $k$ fixed we have $N_{i, k, x_i} = 0$, which can easily happen if there isn't in the dataset any value such that the $i^{\text{th}}$ parameter is $x_i$ and is in the $k^{\text{th}}$ class. Then the whole block's value is $-\infty$ and can never be selected as the class for $x$.
One solution would be to use $\frac{1+N_{i, k, x_i}}{N_k}$ but because this still needs to be a probability we should have the sum equal to 1. Let's note $V_i$ all the possible values for $x_i$ that we assumed to be discrete in this part. Then $$ \sum_{v \in V_i} \frac{1+N_{i, k, v}}{\alpha + N_k} = \frac{1}{\alpha + N_k} \left( |V_i| + \sum_{v \in V_i} N_{i, k, v} \right) = \frac{|V_i| + N_k}{\alpha + N_k} = 1 $$ Because $A_v = \{x \in C_k | x_i = v\}$ is a partition of $C_k$. So $\alpha = |V_i|$. Finally, $$ \begin{align*} \underset{1 \leq k \leq K}{\operatorname{argmax}} \, p(C_k | x) &\approx \underset{1 \leq k \leq K}{\operatorname{argmax}} \, \left[ \log \frac{N_k}{N} - \sum_{i=1}^{n} \log(|V_i| +N_k) + \sum_{i=1}^n \log (1+N_{i, k, x_i}) \right] \\ &= \underset{1 \leq k \leq K}{\operatorname{argmax}} \, \left[ \gamma_k + \sum_{i=1}^n L_{i, k, x_i} \right] \\ \end{align*} $$
Because $\gamma_k$ does not depend on $x$, it can be calculated once and for all at the beginning of the training. Because $x$ can only take a finite number of values, it is also possible to precalculate $L_{i, k, x_i}$ as $n$ two-dimensional matrices of size $V_i \times K$.
for each (x, k) in dataset:
N_k[k] += 1
for i = 1 ... n:
N[i][k][x_i] += 1
for k = 1 ... K:
gamma[k] = log(N_k[k] / N)
for i = 1 ... n:
gamma[k] -= log(|V_i| + N_k[k])
for v in V_i:
L[i][k][v] = log(1 + N[i][k][v])Now let's move on to a simple real-world example. All the code used can be found on my GitHub. I used the Play Tennis dataset, which tries to predict whether a tennis match will take place depending on the weather condition. For a better understanding here are the three first lines of the dataset:
| Outlook | Temperature | Humidity | Wind | Play |
|---|---|---|---|---|
| Sunny | Mild | Normal | Strong | Yes |
| Overcast | Mild | Normal | Weak | Yes |
| Rain | Hot | High | Strong | No |
When $x_i$ is continuous, counting $N_{i, k, x_i}$ no longer makes sense. Instead we'll assume that $p(x_i | C_k)$ is Gaussian, so that $x_i | C_k \sim \mathcal{N}(\mu_{i, k}, \sigma_{i, k}^2)$ where $\mu_{i, k}$ and $\sigma_{i, k}^2$ are simply the empirical mean and variance of the $i^{\text{th}}$ feature of entries in class $C_k$. They can be estimated through $$ \mu_{i, k} = \frac{1}{N_k}\sum_{x \in C_k} x_i \quad \text{and} \quad \sigma_{i, k}^2 = \frac{1}{N_k} \sum_{x \in C_k} (x_i - \mu_{i, k})^2 $$ which gives us the following density: $$ p(x_i | C_k) = \frac{1}{\sqrt{2 \pi \sigma_{i, k}^2}} \exp{\left(- \frac{(x_i - \mu_{i, k})^2}{2 \sigma_{i, k}^2}\right)} $$ Substituting this equality into the general result found earlier, and expanding the square $(x_i - \mu_{i,k})^2 = x_i^2 - 2\mu_{i,k}x_i + \mu_{i,k}^2$, we get $$ \begin{align*} \underset{1 \leq k \leq K}{\operatorname{argmax}} \, p(C_k | x) &= \underset{1 \leq k \leq K}{\operatorname{argmax}} \, \left[ \log p(C_k) + \sum_{i=1}^n \log p(x_i | C_k) \right] \\ &\approx \underset{1 \leq k \leq K}{\operatorname{argmax}} \, \left[ \log \frac{N_k}{N} - \frac{1}{2} \sum_{i=1}^n \log(2 \pi \sigma_{i, k}^2) - \sum_{i=1}^n \frac{(x_i - \mu_{i, k})^2}{2 \sigma_{i, k}^2} \right] \\ &= \underset{1 \leq k \leq K}{\operatorname{argmax}} \, \left[ \underbrace{\log \frac{N_k}{N} - \frac{1}{2} \sum_{i=1}^n \log(2 \pi \sigma_{i, k}^2) - \sum_{i=1}^n \frac{\mu_{i,k}^2}{2\sigma_{i,k}^2}}_{\beta_k} + \sum_{i=1}^n \frac{\mu_{i,k}}{\sigma_{i,k}^2} x_i - \sum_{i=1}^n \frac{x_i^2}{2\sigma_{i,k}^2} \right] \end{align*} $$
Because $\beta_k$ does not depend on the entry $x$, it can precalculated.
function train(dataset):
for each (x, k) in dataset:
N[k] += 1
for i = 1 ... n:
mu[i][k] += x
for k = 1 ... K:
for i = 1 ... n:
mu[i][k] /= N[k]
for each (x, k) in dataset:
for i = 1 ... n:
sigma[i][k] += (x - mu[i][k])^2
for k = 1 ... K:
for i = 1 ... n:
sigma[i][k] /= N[k]
beta[k] = log(N[k] / N)
for i = 1 ... n:
beta[k] -= 1/2 log(2 PI sigma[i][k]^2) + mu[i][k]^2 / (2 sigma[i][k]^2)
return beta, mu, sigmaThe space complexity of the training is $O(nK)$, since $\mu$ and $\sigma$ are simply two $n \times K$ matrices and $\beta$ a vector of size $K$. The time complexity is $O(Nn + Kn)$, i.e. $O(Nn)$ as $N \geq K$ in practice. A prediction is then given in $O(nK)$, exactly like in the discrete case.
Here is a small interactive example with $n=2$ continuous features and $K=5$ classes: five clusters in the plane, each fitted with its own diagonal Gaussian $\mu_k \in \mathbb{R}^2, \sigma_k^2 \in \mathbb{R}^2$, shown as an ellipse at one standard deviation. Drag the black point around to move the query $x$ and see how the score of every class.
In most dataset $\llbracket 0, n \rrbracket = I_{c} \cup I_{d}$ such that for $i \in I_c, x_i$ is continous and for $i \in I_d$, $x_i$ is discrete. The two previous methods can be merged into one single formula that can handle moth the continous and discrete case.
$$ \begin{align*} \underset{1 \leq k \leq K}{\operatorname{argmax}} \, p(C_k | x) &= \underset{1 \leq k \leq K}{\operatorname{argmax}} \, \left[ \log p(C_k) + \sum_{i \in I_d} \log p(x_i | C_k) + \sum_{i \in I_c} \log p(x_i | C_k) \right] \\ &\approx \underset{1 \leq k \leq K}{\operatorname{argmax}} \, \left[ \log \frac{N_k}{N} - \sum_{i \in I_d} \log(|V_i| +N_k) + \sum_{i \in I_d} \log (1+N_{i, k, x_i}) \right. \\ &\left. \quad - \frac{1}{2} \sum_{i \in I_c} \log(2 \pi \sigma_{i, k}^2) - \sum_{i \in I_c} \frac{\mu_{i,k}^2}{2\sigma_{i,k}^2} + \sum_{i \in I_c} \frac{\mu_{i,k}}{\sigma_{i,k}^2} x_i - \sum_{i \in I_c} \frac{x_i^2}{2\sigma_{i,k}^2} \right] \end{align*} $$ Now as usually lets note $\theta_k$ and $L_{i, k, w}$ that can be precalculated during training $$ \begin{align*} \theta_k &= \log \frac{N_k}{N} - \sum_{i \in I_d} \log(|V_i| +N_k) - \frac{1}{2} \sum_{i \in I_c} \log(2 \pi \sigma_{i, k}^2) - \sum_{i \in I_c} \frac{\mu_{i,k}^2}{2\sigma_{i,k}^2} \\ L_{i, k, w} &= \log (1 + |\{x \in C_k : x_i = w\}|) \end{align*} $$ Then we get the much simpler formula $$ \underset{1 \leq k \leq K}{\operatorname{argmax}} \, p(C_k | x) \approx \underset{1 \leq k \leq K}{\operatorname{argmax}} \, \left[ \theta_k + \sum_{i \in I_d} L_{i, k, x_i} + \sum_{i \in I_c} \frac{\mu_{i,k}}{\sigma_{i,k}^2} x_i - \sum_{i \in I_c} \frac{x_i^2}{2\sigma_{i,k}^2} \right] $$ $\theta_k$ and $L_{i, k, w}$ only depend on the dataset, not on the entry $x$ to classify, so just like before they can both be precalculated in a single pass over the training set:for each (x, k) in dataset:
N_k[k] += 1
for i in I_d:
N[i][k][x_i] += 1
for i in I_c:
mu[i][k] += x_i
for k = 1 ... K:
for i in I_c:
mu[i][k] /= N_k[k]
for each (x, k) in dataset:
for i in I_c:
sigma[i][k] += (x_i - mu[i][k])^2
for k = 1 ... K:
theta[k] = log(N_k[k] / N)
for i in I_d:
theta[k] -= log(|V_i| + N_k[k])
for v in V_i:
L[i][k][v] = log(1 + N[i][k][v])
for i in I_c:
sigma[i][k] /= N_k[k]
theta[k] -= 1/2 log(2 PI sigma[i][k]) + mu[i][k]^2 / (2 sigma[i][k])