[x, fval] = curvefit_gumble
function [x,fval] = curvefit_gumble
x_dat = [12;13;14.5;15];
n_dat = [2;2;5;5];
oldopts = optimset('fminsearch');
options = optimset(oldopts,'MaxFunEvals',250,'MaxIter',250,'Display','Off');
[x,fval] = fminsearch(@myfunc,[16, 12, 5, 5, 0.5],options);
function L = myfunc(x)
mu1 = x(1);
mu2 = x(2);
sigma1 = x(3);
sigma2 = x(4);
prop2 = x(5);
F1 = @(x,mu1,sigma1)(exp(-exp(-(x-mu1)/sigma1)));
F2 = @(x,mu2,sigma2)(exp(-exp(-(x-mu2)/sigma2)));
f1 = @(x,mu1,sigma1) ((1/sigma1)*exp(-((x-mu1)/sigma1 + exp(-(x-mu1)/sigma1)))).* (exp(-exp(-(x-mu1)/sigma1)));
f2 = @(x,mu2,sigma2) ((1/sigma2)*exp(-((x-mu2)/sigma2 + exp(-(x-mu2)/sigma2)))).* (exp(-exp(-(x-mu2)/sigma2)));
f_mix = @(x,mu1,mu2,sigma1,sigma2,prop2)(1-prop2)*f1(x,mu1,sigma1)+ prop2*f2(x,mu2,sigma2);
f = @(x,n,mu1,mu2,sigma1,sigma2,prop2)(n.*((1-prop2).*F1(x,mu1,sigma1) + prop2.*F2(x,mu2,sigma2)).^(n-1)).*f_mix(x,mu1,mu2,sigma1,sigma2,prop2);
L = sum(log(f(x_dat,n_dat,mu1,mu2,sigma1,sigma2,prop2)));
end
end
Best Answer