[R] Exact Maximum Likelihood Package

From: Luis David Garcia <lgarcia_at_Math.Berkeley.EDU>
Date: Sat 10 Jul 2004 - 12:04:38 EST

I am a mathematics postdoc at UC Berkeley. I have written a package in a Computational Algebra System named Singular

to compute the Maximum Likelihood of a given probability distribution over several discrete random variables. This package gives exact answers to the problem. But more importantly, it gives All MLE solutions.

My understanding is that current algorithms (like the one in R) only find one solution at a time, and there is no way to decide if that local max is really a global one. That is the power that computer algebra brings to the table.

I have two goals in mind:

1. I would like to create a link between Singular (or any other CAS) and R. The problem of MLE computation is just one instance of the benefits that symbolic computations has to offer to statisticians.

For a more detail account of this interaction, you can check the articles written by Pachter and Sturmfels that will be published in Science

In those articles it is explained how computational algebra can be used to help in many problems of Computational Biology, like Sequence Analysis.

MY PROBLEM IS THAT I JUST STARTED LEARNING R. I really don't know what to do to create such a link, and I would be very interested in having some help/advice on this direction.

2. This is just an example of my ignoRance. I have tried to use R to solve the following MLE problem. But I cannot figure out how to do it. Your help would also be appreciated.

Consider the mixture of a pair of four-times repeated Bernoulli trials. Let s and t be the Bernoulli parameters and p the mixing parameter. There are 5 possible outcomes

```f0 = p*(1-s)^4 + (1-p)*(1-t)^4;
f1 = 4*p*s*(1-s)^3 + 4*(1-p)*t*(1-t)^3;
f2 = 6*p*s^2*(1-s)^2 + 6*(1-p)*t^2*(1-t)^2;
f3 = 4*p*s^3*(1-s) + 4*(1-p)*t^3*(1-t);
f4 = p*s^4 + (1-p)*t^4;

```

The polynomial f_i represents the probability of seeing i successes. Suppose we repeat this experiment 1000 times, and u_i is the number of times we saw i successes.

The likelihood of this event is f_0^u_0*f_1^u_1*f_2^u_2*f_3^u_3*f_4^u_4, and we seek to find those parameter values for s,t,p which maximize the likelihood.

My Singular package has as input the 5 polynomials and a data vector u. For the particular example of u = (3,5,7,11,13). The output are the following
four roots (p,s,t) and the corresponding Likelihoood value:

(0.99998692268562, 0.666609253585517, 5.056808689815264)

0.118420271386274e-27

(0.0000130773153387599, 5.056808689815264, 0.666609254644253)

0.118420271386274e-27

(0.312819438453599, 0.307857785707541, 0.830004221502637)

0.508592337077813e-25

(0.687180561546401, 0.830004221502637, 0.307857785707541)

0.508592337077813e-25

Can anyone tell me how to do the same thing in R. I tried using nlm, but the help on this function is not enough for me, and I cannot get it right.

Thank you very much,

Luis David Garcia

R-help@stat.math.ethz.ch mailing list