naive bayes classifier for non-binary feature values
\begin{align}
&amp;P(x^{(1)}, \ldots, x^{(m)}, y^{(1)}, \ldots, y^{(m)})\\
&amp;=\prod_{i=1}^m P(x^{(i)},y^{(i)})\\
&amp;= \prod_{i=1}^m P(y^{(i)})P(x^{(i)}|y^{(i)})\\
&amp;= \prod_{i=1}^m \left[\textbf{1}(y^{(i)}=1)(1-\theta_0)\prod_{j=1}^nP(x^{(i)}_j|y^{(i)}=1)+\textbf{1}(y^{(i)}=0)\theta_0\prod_{j=1}^nP(x^{(i)}_j|y^{(i)}=1)\right]\\
&amp;=(1-\theta_0)^{\sum_{i=1}^m\textbf1(y^{(i)}=1)}\theta_0^{\sum_{i=1}^m\textbf1(y^{(i)}=0)}\prod_{j=1}^n\prod_{k=1}^s \theta_{j,k|y=1}^{\sum_{i=1}^m\textbf1(x_j^{(i)}=k, y^{(i)}=1)}\theta_{j,k|y=0}^{\sum_{i=1}^m\textbf1(x_j^{(i)}=k, y^{(i)}=0)}
\end{align}

Let's take the logarithm of the expression above and we want to maximize the following 

\begin{align}&amp;\sum_{i=1}^m \textbf1(y^{(i)}=1)\log(1-\theta_0)+\sum_{i=1}^m \textbf1(y^{(i)}=0)\log(\theta_0)+\\&amp; \sum_{j=1}^n \sum_{k=1}^s \sum_{i=1}^m\left[ \textbf1(x_j^{(i)}=k, y^{(i)}=1) \log \theta_{j,k|y=1} + \textbf1(x_j^{(i)}=k, y^{(i)}=0) \log \theta_{j,k|y=0}\right]\end{align}

subject to the condition that 

$$\sum_{k=1}^s \theta_{j,k|y=1}=1$$
$$\sum_{k=1}^s \theta_{j,k|y=0}=1$$

as well as the nonnegative constraints.

Let the Langrange multiplier for the sum to $1$ constraint be $\lambda_{j,1}$ and $\lambda_{j,0}$.

The Langrangian that we want to optimize is 

\begin{align}&amp;\sum_{i=1}^m \textbf1(y^{(i)}=1)\log(1-\theta_0)+\sum_{i=1}^m \textbf1(y^{(i)}=0)\log(\theta_0)+\\&amp; \sum_{j=1}^n \sum_{k=1}^s \sum_{i=1}^m\left[ \textbf1(x_j^{(i)}=k, y^{(i)}=1) \log \theta_{j,k|y=1} + \textbf1(x_j^{(i)}=k, y^{(i)}=0) \log \theta_{j,k|y=0}\right]\\
&amp;-\sum_{j=1}^n \lambda_{j,1}\left(\sum_{k=1}^s \theta_{j,k|y=1}-1\right)-\sum_{j=1}^n \lambda_{j,0}\left(\sum_{k=1}^s \theta_{j,k|y=0}-1\right)\end{align}

Differentiate with respect to $\theta_0$ and equate it to zero, we get 

$$-\frac{\sum_{i=1}^m\textbf1(y^{(i)}=1)}{1-\theta_0}+\frac{\sum_{i=1}^m\textbf1(y^{(i)}=0)}{\theta_0}=0$$

Solving it which is just solving a linear equation, we can check that $\theta_0=\frac{\sum_{i=1}^m\textbf1(y^{(i)}=0)}{m}$.

We can also differentiate with respect to $\theta_{j,k|y=1}$,

$$ \frac{\sum_{i=1}^m\textbf{1}(x_j^{(i)}=k, y^{(i)}=1)}{\theta_{j,k|y=1}}=\lambda_{j,1}$$

That is fixing $j$,$\theta_{j,k|y=1}$ is proportional to   $ \sum_{i=1}^m\textbf{1}(x_j^{(i)}=k, y^{(i)}=1)$. With the sum to $1$ constraint, we have 

$$\theta_{j,k|y=1}=\frac{\sum_{i=1}^m\textbf{1}(x_j^{(i)}=k, y^{(i)}=1)}{\sum_{k=1}^s\sum_{i=1}^m\textbf{1}(x_j^{(i)}=k, y^{(i)}=1)}$$

I will leave the last part as an exercise.