function [x,fval,exitFlag] = sub_Opt3Null(x0,params,data) % explore the behaviour of the likelihood under the null and the alternative. % 5 params under the null. U1=U2 kub = 1; %upper bound on k machinePrecision = 2.220446049250313e-016; %% fmincon % sets up the linear constraints (define bounds Ax < b) A = [-1 0 0 0 0; 1 0 0 0 0 ;0 -1 0 0 0 ;0 1 0 0 0;]; % c's A = [A;0 0 -1 0 0 ; 0 0 1 0 0 ;0 0 0 -1 0; 0 0 0 1 0 ;]; %k's A = [A; 0 0 0 0 -1; 0 0 0 0 1]; %single U b = [-params.clb; params.cub;-params.clb; params.cub; -params.klb; kub; -params.klb;kub; -params.UMax-.0001;params.Uub]; % changes tolerance settings etc. options = optimset('display','off','maxFunEvals',1000*prod(size(x0)),'Tolfun', 1e-6,'tolcon',1e-10,'tolx',1e-10); % call optimization routine under null [x,fval,exitFlag,output] = fmincon(@likelihoodNull,x0,A,b,[],[],[],[],[],options); % value of the likelihood under the null fval = -fval; function like = likelihoodNull(x0) c = x0(1:2); k = x0(3:4); U = x0(5); % U1=U2 % not that if we start this routine using x0 in the unfeasible % region (specifically U < x) then this will not converge. likecontr = []; for ilo = 1:params.numTimes % fmincon will evaluate this function in both the feasible and % infeasible solution space. Once in the feasible solution space it % can still select params in the infeasible solution space. % When in the infeasible space complications arise in the numerics. % specifically when x > 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 -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 end