Numerically compute MLE for parameters in a model with latent variable

In order to initialize my parameters for a bayesian MCMC, I need to compute the MLE of beta and gamma, defined as follows. Y_1, ..., Y_N are observed, Y_i takes one of J ordered categories. There exists a latent continuous random variable Z_i distributed as N(x_i^T,beta) and we observe Y_i, where Y_i=j if gamma_(j-1)

How could I find those MLE with Python?