alpha <- (min(x)+exp(log.del))
I've fixed many convergence problems with simple tricks like this.
hope this helps. spencer graves
Lourens Olivier Walters wrote:
Thanks, your advice worked. I don't have much experience with maths, and therefore tried to stay away from dealing with optimization, but going down to this level opens a lot of possibilities. For the record, the code I used, as you suggested:
############### shape <- mean(data)^2/var(data) scale <- var(data)/mean(data) gamma.param1 <- shape gamma.param2 <- scale log.gamma.param1 <- log(gamma.param1) log.gamma.param2 <- log(gamma.param2) gammaLoglik <- function(params, negative=TRUE){ lglk <- sum(dgamma(data, shape=exp(params[1]), scale=exp(params[2]), log=TRUE)) if(negative) return(-lglk) else return(lglk) }
optim.list <- optim(c(log.gamma.param1, log.gamma.param2), gammaLoglik) gamma.param1 <- exp(optim.list$par[1]) gamma.param2 <- exp(optim.list$par[2]) ################
optim converges and the parameters provide a much better fit to the data
than the method of moments parameters do.
Armed with this new knowledge, I attempted to write a log(likelihood)
optimization method for the Pareto distribution. My ignorance of math
however shows, as the code does not work.
Here is what I did:
################ pareto.density <- function(x, alpha, beta, log=FALSE){ # Pareto distribution not defined for x < alpha # Test for values x < alpha, taking into account floating point error test.vals <- x-alpha error.vals <- test.vals[test.vals<0] if (length(error.vals)>0) stop("\nERROR: x > alpha in pareto.distr(x, alpha, beta, + log=FALSE)") density <- beta * (alpha^beta) * (x^(-beta-1)) if(log) return(log(density)) else return(density) }
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) log.beta.val <- log(beta.val)
paretoLoglik <- function(params, negative=TRUE){ lglk <- sum(pareto.density(data.sorted, alpha=exp(params[1]), + beta=exp(params[2]), log=TRUE)) if(negative) return(-lglk) else return(lglk) }
optim.list <- optim(c(log.alpha.val, log.beta.val), paretoLoglik, + method="L-BFGS-B", lower=c(log.alpha.val, 0), upper=c(log.alpha.val, + Inf)) pareto.param1 <- exp(optim.list$par[1]) pareto.param2 <- exp(optim.list$par[2]) #################
I fixed the alpha parameter as my Pareto density function produces an error if a datavalue > alpha.
I get the following output:
source("browsesessoffplotfitted.R")
Error in optim(c(log.alpha.val, log.beta.val), paretoLoglik, method = + "L-BFGS-B", : non-finite finite-difference value [0]
Any ideas would be appreciated, otherwise its back to method of moments for the Pareto distribution for me :)
Thanks Lourens
On Tue, 2003-09-30 at 22:07, Spencer Graves wrote:
In my experience, the most likely cause of this problem is that optim may try to test nonpositive values for shape or scale. I avoid this situation by programming the log(likelihood) in terms of log(shape) and log(scale) as follows:
> gammaLoglik <- + function(x, logShape, logScale, negative=TRUE){ + lglk <- sum(dgamma(x, shape=exp(logShape), scale=exp(logScale), + log=TRUE)) + if(negative) return(-lglk) else return(lglk) + } > tst <- rgamma(10, 1) > gammaLoglik(tst, 0, 0) [1] 12.29849
Then I then call optim directly to minimize the negative of the log(likelihood).
If I've guessed correctly, this should fix the problem. If not, let us know.
hope this helps. spencer graves
______________________________________________
[EMAIL PROTECTED] mailing list
https://www.stat.math.ethz.ch/mailman/listinfo/r-help
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help
