[R] Exact Maximum Likelihood Package

[Luis David Garcia]

Dear R users,

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

http://www.singular.uni-kl.de

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

http://math.berkeley.edu/~lpachter/papers.html

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