[R-sig-eco] trouble in fitting a distirbution (matteo dossena)

[Edwards, Andrew]

Dear Matteo,

I haven't looked into the details of the mle2 code, but three potential issues spring to mind.

1. It looks like you're trying to find the maximum likelihood estimate (MLE) for each of the two parameters of the Pareto distribution (for which the probability density function is a power-law above some value, call it a, and zero elsewhere). I'll write the probability density function as

f(x) = C x^(- mu), x >= a, where C is the normalisation constant; C = (mu - 1) a^(mu - 1).

The MLEs can be found analytically, you don't need do it numerically. The MLE for mu is just
 
1 - n / (n log a - sum (log x_i) )

where x_i are the data with sample size n.

The MLE for a is just the minimum of your x_i.

2. Trying to estimate the MLE for a numerically may be what's causing trouble. The likelihood function value won't be have a zero gradient at the MLE (unlike what usually happens), because the gradient is always positive.

To see this, I actually have the log-likelihood function given as equation (4) in a recent paper (Using likelihood to test for Lévy flight search patterns and for general power-law distributions in nature, Journal of Animal Ecology 77:1212-1222. If you take that equation and differentiate with respect to a (which I presume is the same as your 'location' parameter), you just get

n (mu - 1) / a  

where n is the sample size and mu>1 is the power-law exponent term (I think the negative of your shape parameter). This derivative is therefore always positive, so does not have a local maximum, and so the maximum likelihood attainable happens when a is as large as possible. However, since a is, by definition, the minimum value of the distribution, the maximum value it can take is the minimum value of your data (since otherwise you would invalidate an assumption of your distribution).  

So if your numerical procedure is looking for a zero gradient of the likelihood function, it won't find it with respect to a, as the gradient is always positive.

3. I'm wondering if the - sign in your line:

> 	density <- -(location+1)*(location^(shape+1))*x^shape

is correct, as I would have thought density should be positive.

I hope this helps. People across several fields have made fitting Pareto distributions much more complicated than it needs to be, but you're on the right track using likelihood. It's just that you can do it analytically, so there's no need for any optimisation. 

Also useful might be: 
White, E.P., Enquist, B.J., & Green, J.L. (2008) On estimating the exponent of power-law frequency distributions. Ecology, 89, 905-912.

Edwards, A.M. (2011) Overturning conclusions of Lévy flight movement patterns by fishing boats and foraging animals. Ecology, 92(6):1247-1257. 

which both also discuss bounded distributions (which can't be solved analytically).

I hope that this helps,

Andy

-- 
Dr. Andrew Edwards 
Andrew.Edwards at dfo-mpo.gc.ca  **NEW** 
http://www.chebucto.ns.ca/~english 

Pacific Biological Station, 
Fisheries and Oceans Canada, 
3190 Hammond Bay Road, 
Nanaimo, British Columbia, 
V9T 6N7, Canada. 
Tel: +1 250-756-7146 
Fax: +1 250 756-7053 

> ------------------------------
> 
> Message: 3
> Date: Wed, 5 Oct 2011 18:19:46 +0100
> From: matteo dossena <matteo.dossena at gmail.com>
> To: r-sig-ecology at r-project.org
> Subject: [R-sig-eco] trouble in fitting a distirbution
> Message-ID: <608BF773-D9AA-4CFD-A5C1-2C510ED0344A at gmail.com>
> Content-Type: text/plain
> 
> Hi All,
> 
> I'm trying to fit a theoretical distribution to some data, 
> following the procedure on Bolker's Ecological Models and 
> Data in R, but I'm having some problem... could someone 
> please point me out what I'm doing  wrong?
> 
> here is the code
> 
> #create a random variable 
> x <- sort(rpareto(100,1,0.5))
> 
> #the pdf
> pdf.pareto <- function(x,location,shape){
> 	density <- -(location+1)*(location^(shape+1))*x^shape
> 	return(density)}
> 
> #predict the probability according to my pdf
> y <- pdf.pareto(x,1,0.5)
> 
> #check that pdf function does what it should
> plot(x,y,typ="l")
> 
> #the negative log-likelihood function
> NLL.pareto <- function(location,shape){
> 	density <- pdf.pareto(x,location,shape)
> 	-sum(dpareto(density,location=location,shape=shape,log=TRUE))}
>  
> #run the optimisation with mle2 to derive the parameters of the pdf
> fit.pareto <- mle2(NLL.pareto,start=list(location=1,shape=0.5))
> 
> 
> when i run it, I get this message:
> Error in optim(par = c(1, 0.5), fn = function (p)  : 
>   initial value in 'vmmin' is not finite
> 
> I was almost convinced that the reasoning behind was correct, 
> and to show it to myself i did the same thing with a less 
> "exotic" distribution. In that case the code doesn't crash 
> but the estimate i get are pretty weird...I'm starting to 
> think that I'm trying to do something totally wrong...  
> really appreciate any help
> thanks
> matteo
>