2
$\begingroup$

I am applying a MAP log-likelihood approach in order to fit a Markov mixture model, where objective function to be maximized is given by the formula:

$$ L(X|\Theta _K)=\sum_{i=1}^{n}f(X_i|\Theta_K)+\sum_{j=1}^{K}\sum_{n=0}^{M}\log p(\theta_n^{j}|a_n^{j}) $$

where the second argument is a sum of Dirichlet priors ( I am using the formula given by Wikipedia) and the first argument is the sum of log-likelihood across all sequences and all components.

At this point I have achieved a lot in implementing the algorithm in R, thanks to answers to my previously posted questions related to topic.

At this stage my question is - after performing one step of expectation-maximization algorith ( with 2 components to start with), my value of $L(X|\Theta_K)$ became much smaller than it was ( from -2200 to -8000). I believe that my code is correct and do not understand why this could be happening ( the next 2 steps show steady increase). Can there be fluctuation of the algorithm in the beginning?

There are 2 possible issues, however I cannot pinpoint if they are indeed the causing this: underflow problem(some of the values in my transition matrix and the resulting likelihood values, can be so negligibly small that in certain calculations R, in which I am working, rounds them up and thus I end up with incorrect calculations; or meaningless Dirichlet priors ( e.g. sometimes the priors are larger than 1).

Additionally: the resulting from M-step multinomial distributions for each row of the transition matrix and start probabilities vector sum to 1 for each of the components. The posterior conditional probabilities of the hidden variables ( components ) also seem to make sense.

The example of conditional posteriors of hidden variables for first ([[1]]) and second [[2]] component for each of 50 sequences are below ( results of first iteration E-step):

    [[1]]
     [1] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00   1.520433e-10 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
    [11] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00  1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
    [21] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
    [31] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 2.726330e-02 1.000000e+00 1.000000e+00
    [41] 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.979822e-01 1.000000e+00 1.000000e+00

     [[2]]
     [1]  7.324479e-65  2.462187e-97  4.146568e-32  3.498317e-97 2.135013e-274  1.000000e+00  7.731884e-47 2.068553e-264 8.497501e-260
     [10] 4.356689e-271 2.983088e-213 7.485688e-110 4.556750e-287 1.360219e- 173  2.340220e-45  2.609916e-59  6.057344e-59 1.286382e-185
     [19]  3.879706e-80 8.488843e-188  1.881308e-14 3.098226e-290 1.681928e-290 1.168018e-211 2.123491e-292 5.767748e-177 9.232827e-198
    [28] 1.120970e-159 1.397181e-257  4.078388e-48 1.524531e-247  1.000000e+00 1.833904e-252 1.452165e-263 1.878481e-111 3.379251e-178
    [37] 1.823290e-247  9.727367e-01 8.553065e-238  2.748773e-32 4.602824e-138  1.212533e-93 1.744806e-271 2.677587e-292 7.822883e-131
     [46] 2.504779e-111  1.775144e-16  8.020178e-01  1.301381e-63 4.437099e-114
$\endgroup$
  • $\begingroup$ If you continue the iterations, do you achieve convergence? $\endgroup$ – gregory_britten Sep 9 '13 at 14:42
  • $\begingroup$ I am at the 4th iteration now (my algorithm is running very slowly at the moment), and after initial dip, log likelihood has been increasing so far ( and one of the components is getting updated). I notice that for my first mixture component, whose transition probability matrix was initialized by summarizing all transitions in the collection and normalizing by the total number of transitions, has stagnated after the first iteration ( as if there is no more room to improve) - the other component is still changing ( and correspondingly prior probabilities of components rebalance). $\endgroup$ – zima Sep 9 '13 at 14:50
  • $\begingroup$ @gregory_britten, after the 4th iteration my MAP log-likelihood is again smaller. Am I right in my understanding that m-step takes first and second derivative of the function and thus the scenario I am describing ( fluctuation ) is impossible? $\endgroup$ – zima Sep 9 '13 at 15:42
  • $\begingroup$ Did you write your own code? If you did, this is probably a bug. $\endgroup$ – StasK Sep 9 '13 at 15:49
  • 2
    $\begingroup$ It appears something is wrong. As far as I know, the EM is guarented to improve the likelihood, and never decrease it. This is reiterated in the following document which I found quickly online here. I suspect there is a bug in your code. $\endgroup$ – gregory_britten Sep 9 '13 at 15:56
2
$\begingroup$

In applications of the expectation-maximization algorithm, the likelihood of the data should monotonically increase.

As discussed by R. Neal and G. Hinton the EM algorithm can be seen as (effectively) a gradient ascent algorithm in where the data likelihood is the objective function expressed as a function of the model parameters..

The main case I can envision where this type of problem would arise (aside from implementation bugs) is when one is using an approximation technique to solve for the values of the parameters that yield the expectation and/or optimization. For example, if one has to, due to the structure of the distributions in question, use approximations for the values of the parameters that maximize the likelihood; then these approximations may allow one to "jump across" a local maximum, just like in any other application of gradient-following algorithms.

For the specific case of apply the EM algorithm for estimating a two-component discrete Markov model, you should be able to evaluate the expectations and do the maximizations exactly, so you shouldn't be having this problem.

$\endgroup$
  • $\begingroup$ thank you @Dave for the thorough answer. I am trying to understand what is the problem.. I also notice that the second term of my MAP-function jumps around with every iteration.. Should the second term also always increase? When the model parameters include extremely small numbers, and alpha concentration parameters are also very small, I get sometimes results larger than 1. I also have to approximate the results of estimation ( because there are numbers like 1e-60 which sometimes get converted to 0 by R. I updated my question with an example of E-step values. Could you perhaps take a look? $\endgroup$ – zima Sep 9 '13 at 17:32
  • $\begingroup$ I also have to mention that I did approximate 0's in my model parameters with 1e-17 ( which later got updated with even smaller values). $\endgroup$ – zima Sep 9 '13 at 17:38
  • $\begingroup$ Three possibilities come to mind: (1) you are not normalizing some of your computations, and thus driving the prior of the second component to 0, (2) you may be trying to fit a mixture when a a single-component model is better , (3) in these situations I like to express the prior as a set of observations, and then treat them just like the actual data. $\endgroup$ – Dave Sep 9 '13 at 18:10
  • $\begingroup$ ,is it possible to elaborate a bit on (3)? Trying to understand what this conceptually means. As far as (1) and (2) - I have verified the code multiple times for normalization.. Still looking to hopefully find a bug,so far to no avail. For (2) - second component had a lot of very small parameter values initially because it was plainly computed from one of the sequences and thus didn't contain a lot of the states in the model. But with 2nd and 3rd iteration this changed ( e.g. the transition probabilities were more even across the states). The sums of rows are always 1 for each iter-n. $\endgroup$ – zima Sep 9 '13 at 18:19
  • $\begingroup$ I have had a thought yesterday regarding the approximations.. I replace 0's in any calculations with the value of 1e-17, so 1e-17 is informal for 0 in my model. But, sometimes I see as a result of an iteration, values much lower than 1e-17 start appearing. So perhaps, I can assume that everything that is =< 1e-17 is also 0? Would this increase any validity to my results? Is it at all feasible to work with numbers like 1e-300 while modeling in R? The program itself replaces such small numbers with 0's in some calculations ( like addition with a much larger number) $\endgroup$ – zima Sep 10 '13 at 9:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.