day20150119_P1   day20150119_P1.txt  code:

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

@                                                                              @

@       Example GAUSS Code for Maximum Likelihood Estimation of Poisson Model  @

@                                                                              @                                                                          

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

@                                                                              @

@       Note: MAXLIK is a special routine and you have to tell GAUSS that your @

@             will be using it                                                 @

@                                                                              @

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

new;

cls;

load path= D:Pro_gauss;

library maxlik;

maxset;

_max_Algorithm = 4;

output  file = poisson.out reset;

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

@                                                                              @

@       The first step in this program is to generate the pseudo-dataset       @

@                                                                              @

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

nobs    = 200;

x       = (rndp(nobs,1,6) - 6) .^2;

c       = ones(nobs,1);

beta1   = 1;

beta2   = -0.25;

lambda  = exp(beta1 + beta2.*x);

y       = zeros(nobs,1);

for     i (1,nobs,1);

       y[i]     = rndp(1,1,lambda[i,1]);

endfor  ;

dta     = y~c~x;

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

@                                                                              @

@       Procedure to compute loglikelihood value for poisson model             @

@                                                                              @

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

local   yvar, xvar, qvar, llik;

yvar    = dta[.,1];

xvar    = dta[.,2:3];

llik    = yvar .* (xvar * b) - exp(xvar * b) - lnfact(yvar);/* Log-Likelihood */

retp    (llik);

endp;

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

@                                                                              @

@       The next step is to invoke the MAXLIK routine                          @

@       MAXPRT is an auxiliary routine that formats the output of MAXLIK       @

@                                                                              @

@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

stval   = {1,0};                          /*Starting values */

{b,logl,g,vc,ret}=maxprt(maxlik(dta,0,&pmod,stval));

How can I write a maximum likelihood routine in R?

Create a function (find.mle in this example) that takes a vector of data and calculates the MLE based on it, and then use apply to apply that to the columns of data:

library("maxLik")

data <- replicate(20, rnorm(100))

find.mle = function(d) {

   logLikFun <- function(param) {

       mu <- param[1]

       sigma <- param[2]

       sum(dnorm(d, mean = mu, sd = sigma, log = TRUE))

   }

   maxLik(logLik = logLikFun, start = c(mu = 0, sigma = 1))$estimate

}

mles = apply(data, 2, find.mle)

This will give you a 2x20 matrix with your estimates:

> mles

           [,1]      [,2]        [,3]       [,4]       [,5]        [,6]

mu    0.03675611 0.1129927 -0.06499549 0.04651673 0.06593217 -0.08753828

sigma 0.93497523 0.9817961  0.84734600 0.93139761 1.01083924  1.04114752

          [,7]       [,8]      [,9]       [,10]      [,11]       [,12]

mu    0.1629807 0.01665411 0.2306688 -0.02147982 0.07723695 0.009476477

sigma 1.0428713 1.01658241 1.0073277  0.99781761 0.99327722 0.983356049

          [,13]      [,14]      [,15]      [,16]     [,17]     [,18]

mu    0.06524147 0.02442983 -0.1305258 -0.1050299 0.1449996 0.1172218

sigma 1.04004799 0.89963009  0.9979824  1.0227063 0.9319562 0.9916734

          [,19]       [,20]

mu    -0.1288296 -0.05769467

sigma  0.9975368  0.89506586

thanks for that, i am actually looking for the mu and sigma for all 20 samples. the summary(mle) already gives me the overall mu and sigma but i want it disaggregated for each of hte 20 different samples. Thanks –  YesSure Jul 10 '12 at 19:06

Ah! In that case, my edited answer is what you're looking for. –  David Robinson Jul 10 '12 at 19:21

very nice, i had just discovered a couple of these terms and was working towards this but would have taken me a lot longer. Interestingly when i call mle.estimates, R is telling me it cant be found. –  YesSure Jul 10 '12 at 19:49

any ideas on my query above? –  YesSure Jul 10 '12 at 21:01

Sorry- I meant to put mles (the name of the variable I'd created in the previous code, I had changed it mid posting) –  David Robinson Jul 10 '12 at 21:52

I really think that there is not need to write any function in order to obtain the maximum likelihood (ML hereafter) estimators for the mean and sd. If X is a normal random variable then the ML estimators for the population mean and sd are the sample mean and the sample sd, and we know that the sample mean is an unbiased estimator for the population mean but the ML estimator for the variance is biased (downward), since the denominator for the variance is n instead of being n-1.

So R calculates the sample quasi-variance (corrected for the degrees of freedom) and this is the unbiased estimator, so it is not the ML estimator, but we can obtain the ML estimator from the R estimation, simply we only have to multiply it by (n-1)(1/n) and the result will be the ML estimation of the variance, then apply squared root and voilá you'll have the ML estimation for the sd, but I like easy stuff so, just multiply the sd by (n-1)(1/n) and this is your answer. For a detailed explanation see Population variance and sample variance on http://en.wikipedia.org/wiki/Variance

I want to run some maximum likelihood code on the data sample I have created. This is what I have so far:

library("maxLik")

data <- replicate(20, rnorm(100))

logLikFun <- function(param) {

mu <- param[1]

sigma <- param[2]

sum(dnorm(data, mean = mu, sd = sigma, log = TRUE))

}

mle <- maxLik(logLik = logLikFun, start = c(mu = 0, sigma = 1))

summary(mle)

I am having some problems extracting the mean and standard deviation for each sample of the 20, I amended the apply function to try to suit this but nothing has worked yet. Any ideas?