How do I use MLE on a shifted gamma distribution?

First the Matlab documentation on using the built in distributions is great. I've successfully used gampdf to fit distributions using MLE. However, I would now like to use a custom distribution and I'm running into problems. Error messages such as the "...function returns negative or zero values" are very common when working through. I tried looking at gampdf.m but there's a lot of code in that function: example what is the scaling factor z=x./b?

here is my function below just to start off, any help getting this to work with MLE (in error free form) would be greatly appreciated. Thanks in advance.

function y=mygampdf(x,alpha,beta,x1)
%
% Y = MYGAMPDF(X,ALPHA,BETA,XI);
% This is a shifted gamma function along the x-axis to the right using the
% term XI.  All other factors are the same as the usual gamma distribution
% in gampdf.m
for i=1:length(x)
  y(i)=0;
end
y=(1/(beta^alpha*gamma(alpha))).*(x-x1).^(alpha-1).*exp(-1.0.*(x-x1)./beta);
end

Answer 1

[Peter Perkins]

Pat, thresholded distributions are typically not easy to fit by maximum likelihood. I am not sure of the details for a 3-param gamma, there may be literature specifically dealing with this, I don't know. But for a 3-param lognormal, the MLE is unbounded and in effect estimates the threshold parameter at the smallest observation and the variance parameter at zero, while for the 3-param Weibull, the MLE is unbounded in some cases and in others there are multiple local optima.

You can try fitting by maximum likelihood, but if you're using the MLE function with a custom PDF function, you at least will need to upper bound the threshold parameter by the smallest observation, and probably that minus a small epsilon. You might consider looking at this

http://www.mathworks.com/products/statistics/examples.html?file=/products/demos/shipping/stats/cdffitdemo.html

which describes another fitting method that doesn't suffer from the problems that ML has in some thresholded distributions.

By the way, it's not clear to me why you need that loop in your PDF. If you want to preallocate y, use y = zeros(size(x)). But given the line that follows the loop, you shouldn't even need to preallocate. Hope this helps.

[Pat]

Thanks Peter! I will look at the link closely. I've been using a double gamma on my fit y=P*gampdf(x,a1,b1)+(1-P)*gampdf(x,a2,b2) and this seems to be doing a good job. The main distribution of data comes with a long tail hence the need for a double gamma. I was going to try a double shifted gamma next, but it sounds like I want to tread very carefully here.

[sepideh]

I am trying to do the same thing and this link apparently should be purchased? What is the fitting method? I am trying to find it from another source if possible.

[Peter Perkins]

Try this link:

http://www.mathworks.com/help/stats/examples/fitting-a-univariate-distribution-using-cumulative-probabilities.html?prodcode=ST&language=en

Answer 2

[Matt Tearle]

Peter is the expert on statistical matters, so I'd take his advice on whether what you're trying to do is a good idea. But here's one way to make something work (whether or not it's a good idea!):

% define the pdf
shiftgampdf = @(x,a,b,xi) max(realmin,gampdf(x-xi,a,b));
% make some data
a = pi;
b = 2/pi;
y = gamrnd(a,b,200,1)+10;
% do the MLE fit
mle(y,'pdf',shiftgampdf,'start',[2,1,5])

Note that I'm defining the PDF as a function handle that simply calls gampdf with a shifted x value. However, you can get values that are 0 to machine precision during the MLE fitting, which is why I've added the max(realmin... part. This ensures that the values returned are always very slightly positive (avoiding the error messages that you get from mle).

[Pat]

Thanks Matt! I appreciate your input. I will give your solution a shot, and I agree with you, Peter is certainly the expert. His statements above make me nervous about proceeding but we'll see where this goes. Thanks again.