function [x,fval,exitFlag] = sub_Opt3Alt(x0,params,data) % explore the behaviour of the likelihood under the alternative. % 6 params under the alt. U1 U the rate should be zero but numerically it can be % complex, -ve or positive. if unchecked this causes problems in the calculation of p. % for x < U the rate is only positive and thus 0 < p < 1 and all is fine % use the latutudes that are < U. For values of x>U the likelihood = 0. % if xMin(ilo) < U then select the data upto this latitude ind = (params.xLat params.mu); x = params.xLat(ind); % selects the data for latitudes > U(ilo). % these latitudes will all have rates of zero and probabilities of zero % if data = zero at these latitudes the contribution to the % likelihood is zero. if data exists the contribution to the % likelihood is -infinity. we will use indInfeasible = logical(1-ind); if sum(data(ilo,indInfeasible) > 0) % then there is at least one latitude > U that has non zero % data. its contribution to the likelihood is -inf likelihoodAddition = -inf; else % we add nothing likelihoodAddition = 0; end % rate will always be positive if x>U. ratio = (x-params.mu)/(U(ilo)-params.mu); rate = c(ilo)*power(1-ratio,(1/k(ilo))-1); % rate must be +ve to give 0 < p <1 probi = 1-exp(-rate); % adjust extreme results as they will cause the nans and inf values probi(probi == 0) = machinePrecision; % make as small as possible but not actually zero. % this probability is never really bigger than 1 but will reach % 1 asymptotically. when this happens make it as close to 1 as % possible so the likelihood can be calculated as > -inf probi(probi >= 1) = 1-machinePrecision; % need to maximize likecontrib = (data(ilo,ind).*log(probi) + (params.nReplicates(ind)-data(ilo,ind)).*log(1-probi)); likecontr = [likecontr,likecontrib,likelihoodAddition]; end % NaNs should never arise. cases where probi = 0 and data = 0; % check likecontr(isnan(likecontr)) = 0; like = -sum(likecontr); if isnan(like) error('nan for likelihood. big trouble!') end end function [c,ceq] = myconAlt(x0) % these are the contraints for U's. % U1 < U2 % the last 3 elements of x0 are the U's c = [ x0(5)-x0(6) ]; %They need to be ordered ceq = 0; % set to zero as we want to skip this part. we have no equality contraints end end