In the example you have below, two of the
solutions should be rejected since the parameters
are outside of the valid range and the remaining
two can be derived from each other due to the
symmetry of the problem so actually you only need
one solution in this case.
At any rate, you can create a coarse grid using
grid.expand and apply the likelihood function using apply
to evaluate your likelihood on each point of the
grid. Try running optim a few times starting at
each of the top few values found in your grid search.
You will want to give optim derivatives and for this
R has deriv to calculate symbolic derivatives.
See:
?expand.grid
?apply
?deriv
?optim
Luis David Garcia <lgarcia <at> math.berkeley.edu> writes:
:
: 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