[R] fitdistr, mle's and gamma distribution

Lourens Olivier Walters lwalters at cs.uct.ac.za
Thu Oct 9 10:31:55 CEST 2003


Thanks, again I managed to get convergence, although looking at the
resultant Pareto distribution plotted in R, it doesn't seem that the log
likelihood estimated parameters are more accurate than the methods of
moments estimator parameters - I will have to look into it. For the
record, here is the code I used: 

###########

data.sorted <- sort(data)
alpha.val <- data.sorted[1]
beta.val <- 1/((1/n) * sum(log(data.sorted/alpha.val)))

log.alpha.val <- log(alpha.val)
# Optimize log.del = log(0.0000001), which is the log(difference 
# between alpha and log likelihood estimated alpha), chose small 
# difference to start of with
log.del = -16.11809565

paretoLoglik <- function(params, negative=TRUE){
  alpha <- data.sorted[1]+exp(params[1])
  lglk <- sum(pareto.density(data.sorted, alpha=exp(alpha),
+ beta=exp(params[2]), log=TRUE))
  if(negative)
     return(-lglk)
  else
     return(lglk)
}

optim.list <- optim(c(log.del, log.beta.val), paretoLoglik,
+ method="BFGS")
pareto.param1 <- exp(optim.list$par[1]) + alpha.val
pareto.param2 <- exp(optim.list$par[2])

##########

On Wed, 2003-10-08 at 15:07, Spencer Graves wrote:
>       I'm sorry, but I don't have time to read all your code. 
However, 
> I saw that you tested for x > alpha in your Pareto distribution 
> example.  Have you considered reparameterizing to estimate log.del = 
> log(alpha-min(x))?  Pass log.del as part of the vector of parameters
to 
> estimate, then immediately inside the function, compute
> 
>       alpha <- (min(x)+exp(log.del))
> 
>        I've fixed many convergence problems with simple tricks like
this. 
> 
>       hope this helps.  spencer graves




More information about the R-help mailing list