function [x,fval,exitFlag] = sub_Opt1(x0,params,data) % explore the behaviour of the likelihood under the null and the alternative. % 5 params under the null. U1=U2 % 6 params under the alt. U1 0); xMin = xGTZero(end); machinePrecision = 2.220446049250313e-016; %% fmincon % sets up the linear constraints (define bounds Ax < b) A = [-1 0 0; 1 0 0; 0 -1 0; 0 1 0; 0 0 -1; 0 0 1]; b = [-params.clb; params.cub; -params.klb; params.kub; -params.Ulb(params.TS)-.0001;params.Uub]; % changes tolerance settings etc. possibly need to tweak these options = optimset('display','off','maxFunEvals',1000*prod(size(x0)),'Tolfun', 1e-6,'tolcon',1e-6,'tolx',1e-6); % call optimization routine under null [x,fval,exitFlag,output] = fmincon(@likelihood,x0,A,b,[],[],[],[],[],options); % value of the likelihood under the null fval = -fval; function like = likelihood(x0) ic = ic + 1; c=x0(1); k=x0(2); U=x0(3); % 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:numTimes % always 1 % 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 < U(ilo)) & (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(ind).*log(probi) + (params.nReplicates(ind)-data(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