[R] Exact Maximum Likelihood Package
spencer.graves at pdf.com
Sat Jul 10 16:07:48 CEST 2004
Beyond this, I believe it may be possible to write an interface
between R and Singular; see "Writing R Extensions", available from
www.r-project.org -> Documentation: Manuals or via help.start() from
within and R session. There already exists an interface between Excel
and R; see www.r-project.org -> CRAN -> (select a local CRAN Mirror,
e.g., http://cran.stat.ucla.edu/) -> Software: Other -> R-(D)COM
I have not done this, but I suspect it should be possible to
interface R and Singular, and there my even be folks in your local
Statistics Department who have done crudely similar things.
hope this helps. spencer graves
p.s. Gabor mentioned both "grid.expand" and "expand.grid";
"expand.grid works, while "grid.expand" seems to be a typo.
Gabor Grothendieck wrote:
>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.
>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
>: 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
>: 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
>: 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
>: four roots (p,s,t) and the corresponding Likelihoood value:
>: (0.99998692268562, 0.666609253585517, 5.056808689815264)
>: (0.0000130773153387599, 5.056808689815264, 0.666609254644253)
>: (0.312819438453599, 0.307857785707541, 0.830004221502637)
>: (0.687180561546401, 0.830004221502637, 0.307857785707541)
>: 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 at stat.math.ethz.ch mailing list
>PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
More information about the R-help