Fitting a Gamma distribution with threshold parameter: error message?

Hey,

I'm trying to fit a Gamma distribution to some raw data which includes a lower threshold parameter, and i'm getting the following error message. Could somebody help here please?

CODE %%%%%%%%%%%%%%%%%

c = min(data)/4
k = mean(data)^2 
theta = var(data)/mean(data)^2

x = sort(data);
n = length(x);
pEmp = ((1:n)-0.5)' ./ n;
wgt = 1 ./ sqrt(pEmp.*(1-pEmp));
gammaObj = @(params) sum(wgt.*(gamcdf((x-params(1)),exp(params(2)),exp(params(3))) - pEmp).^2)
paramHat = fminsearch(gammaObj,[c,log(k),log(theta)])
paramHat1 = paramHat(1);
paramHat2 = exp(paramHat(2));
paramHat3 = exp(paramHat(3));

ERROR %%%%%%%%%%%%%%%%%

|??? Error using ==> minus Matrix dimensions must agree.

Error in ==> @(params)sum(wgt.*(gamcdf((x-params(1)),exp(params(2)),exp(params(3)))-pEmp).^2)

Error in ==> fminsearch at 205 fv(:,1) = funfcn(x,varargin{:});

Error in ==> gammaFit at 582 paramHat = fminsearch(gammaObj,[c,log(k),log(theta)])

| Thanks a bundle !

Answer 1

[Nirmal Gunaseelan]

The error points to the MINUS operator. It appears that the variables x and n (and hence pEmp?) have the same dimensions and params(1) is different from this. What are the dimensions of these variables?

[Paul]

they should all be scalers!

Here is the weird thing. Running this test code everything works well :

n = 50;

sim_data = gamrnd(2,1,n,1)+2;

x = sort(sim_data);

pEmp = ((1:n)-0.5)' ./ n;

a1 = mean(x)^2;

b1 = var(x)/mean(x);

c1 = min(x)*0.90;

wgt = 1 ./ sqrt(pEmp.*(1-pEmp));

gammaObj = @(params) sum(wgt.*(gamcdf((x-exp(params(3))),exp(params(1)),exp(params(2)))-pEmp).^2);

paramHat = fminsearch(gammaObj,[log(a1),log(b1),log(c1)]);

paramHat = exp(paramHat)

Well, the fit may not be great but that's another story. Replacing sim_data with some real data I get the error above?

[Nirmal Gunaseelan]

Good. So the difference now is between the dimensions of sim_data and the real_data that you are using. You could use this:

size(sim_data)

to verify.

[Paul]

very simple:

real_data was not properly vectorized !

Thanks a lot.

Answer 2

[haneen jawa]

(R .* (C.^3)/2) .*(3/K) .*gamma(3/k); WHAT IS THE WRONG