function argout = Bootmodel(varargin)
% Bootstrap s.e. on age models, plus PB_score

%INPUTS
% De,model,[sigb err],start
% Where De      is in two columns
%       model   is model name in sinlge quotes:
%                'cam'
%                'mam3'
%                'mam3ul'
%       [sigb err] are the sigma_b and error. The error may be ommitted.
%       start   is the starting estimate for mam3 or mam3ul

% To use the CAM, enter
% argout = Bootmodel(de,'cam')
% any other input terms are ignored if the model is 'cam'
% the output is in the format [central_age error overdispersion error]


% To use the MAM3, enter one of these formats:
% argout = Bootmodel(de,'mam3',0.2)    % for sigma_b of 20%
% argout = Bootmodel(de,'mam3',0.2,1)  % for sigma_b of 20% and start_val of 1
% argout = Bootmodel(de,'mam3',[0.2 0.03],1)  % for sigma_b of 20+-3% and start_val of 1
%
% The output format is [min_age error sigma P]
% see Galbraith 99 for definition of sigma and P


% To use the MAM3ul (with PB score calculation), enter one of the following:
% argout = Bootmodel(de,'mam3ul',0.2)    % for sigma_b of 20%
% argout = Bootmodel(de,'mam3ul',0.2,1)  % for sigma_b of 20% and start_val of 1
% argout = Bootmodel(de,'mam3ul',[0.2 0.03],1)  % for sigma_b of 20+-3% and start_val of 1
%
% The output format is [min_age error sigma P PB_score error]
% see Galbraith 99 for definition of sigma and P

%--------------------------------------------------------------------------
argout=[];
%Extract the input arguments from cell array
if isstruct(varargin{1})
    struc=varargin{1};
    de=struc.de;
    n_samples=length(struc);
else
    de = varargin{1};
    n_samples=1;
end

model=varargin{2};

for j=1:n_samples
n=length(de(:,1));
nboot=200;


%Validation of input arguments
if ischar(model)==0
    error('Model name must be a string, use single quotes')
end

TF = strcmp('CAM',model);
TF2 = strcmp('cam',model);
TF3 = strcmp('MAM3',model);
TF4 = strcmp('mam3',model);
TF5 = strcmp('MAM',model);
TF6 = strcmp('mam',model);
TF7 = strcmp('MAM3UL',model);
TF8 = strcmp('mam3ul',model);
TF9 = strcmp('MAMUL',model);
TF10 = strcmp('mamul',model);
TF11 = strcmp('MAM3ul',model);
TF12= strcmp('MAMul',model);
t1=TF||TF2||TF3||TF4||TF5||TF6||TF7||TF8||TF9||TF10||TF11||TF12;
if t1==0;
    error('Model name is not recognised')
end


% Using CAM

TF = strcmp('CAM',model);
TF2 = strcmp('cam',model);
if or(TF, TF2)==1;
    outp=zeros(nboot,4); %initialise matrix
    for i=1:nboot
        R = unidrnd(n,n,1);
        Rde=de(R,1:2);
        outp(i,:) = getCAM(Rde);
    end
    meanTheta=mean(outp(:,1));
    seTheta=std(outp(:,1));
    meanOD=mean(outp(:,3));
    seOD=std(outp(:,3));
    
    if isstruct(varargin{1})
        struc(j).bootCAM=[meanTheta seTheta meanOD seOD];
        argout=struc;
    else
        argout=[meanTheta seTheta meanOD seOD];
    end
end


% using MAM3
TF = strcmp('MAM3',model);
TF2 = strcmp('mam3',model);
TF3 = strcmp('MAM',model);
TF4 = strcmp('mam',model);


if TF||TF2||TF3||TF4==1;
    outp=zeros(nboot,4); %initialise matrix
    if length(varargin)<3
        error('Not enough input arguments')
    else
        Psigb=varargin{3};
    end
    
    if isstruct(varargin{1}) && isfield(struc,'start')
            start=struc(j).start;
    else
        if length(varargin)==4
            start=varargin{4};
        else
            % generate a reasonable starting value
            if length(de(:,1))>2
                de2=sort(de,1);
                k=ceil(length(de(:,1))/2);
                cam=getCAM(de2(1:k,:));
                start=cam(1);
            else
                cam=getCAM(de);
                start=cam(1);
            end
        end
    end
    
    if length(Psigb)==1
        Psigb=[Psigb 0]; % add zero std on sigb for easy processing
    end
    for i=1:nboot
        R = unidrnd(n,n,1);
        Rde=de(R,1:2);
        sigb=randn*Psigb(2)+Psigb(1);
        outp(i,:)=getmam3(Rde,sigb,start);
    end
    
    if isstruct(varargin{1})
        struc(j).bootMAM3=[mean(outp(:,1)) std(outp(:,1)) mean(outp(:,3)) mean(outp(:,4))];
        argout=struc;
    else
        argout=[mean(outp(:,1)) std(outp(:,1)) mean(outp(:,3)) mean(outp(:,4))];
    end
    
end
        

% using MAM3ul

TF = strcmp('MAM3UL',model);
TF2 = strcmp('mam3ul',model);
TF3 = strcmp('MAMUL',model);
TF4 = strcmp('mamul',model);
TF5 = strcmp('MAM3ul',model);
TF6 = strcmp('MAMul',model);

if TF||TF2||TF3||TF4||TF5||TF6==1;
    outp=zeros(nboot,5); %initialise matrix
    if length(varargin)<3
        error('Not enough input arguments')
    else
        Psigb=varargin{3};
    end
    
    if isstruct(varargin{1}) && isfield(struc,'start')
            start=struc(j).start;
    else
        if length(varargin)==4
            start=varargin{4};
        else
            % generate a reasonable starting value
            if length(de(:,1))>2
                de2=sort(de,1);
                k=ceil(length(de(:,1))/2);
                cam=getCAM(de2(1:k,:));
                start=cam(1);
            else
                cam=getCAM(de);
                start=cam(1);
            end
        end
    end
    
    if length(Psigb)==1
        Psigb=[Psigb 0]; % add zero std on sigb for easy processing
    end
    
    for i=1:nboot
        R = unidrnd(n,n,1);
        Rde=de(R,1:2);
        sigb=randn*Psigb(2)+Psigb(1);
        outp(i,1:4)=getmam3ul(Rde,sigb,start);
        
        %and get PBscore for this bootstrap sample
        ests=get_sigma_ul(Rde,sigb,outp(i,1)); %ests = [sigma p0] with gamma and ACMOD sig_b
        outp(i,5)=(1-ests(2))*sqrt(ests(1)); % The PB score
    end
    
    if isstruct(varargin{1})
        struc(j).bootMAM3ul=[mean(outp(:,1)) std(outp(:,1)) mean(outp(:,3)) mean(outp(:,4)) mean(outp(:,5)) std(outp(:,5))];
        struc(j).PBscore=[mean(outp(:,5)) std(outp(:,5))];
        argout=struc;
    else
        argout=[mean(outp(:,1)) std(outp(:,1)) mean(outp(:,3)) mean(outp(:,4)) mean(outp(:,5)) std(outp(:,5))];
    end
end

end

%-------------------- Auxilary functions----------------------------------

function ests=get_sigma_ul(de,sigma_b,gamma)
% used in calculation of PB_score

% Starting values
p0=0.5;         % proportion at gamma
%gamma=start;
%mu=gamma;%mean(de(:,1));     % the other normal not used in mam3
sigma=std(de(:,1));  % s.d of the half-normal

% lower and upper bounds [gamma sigma p0]
xlb=[0.01 0.01];
xub=[100 1];

dat=[de(:,1); de(:,2); sigma_b; gamma]; % format data into 1 vector, need for passing to mle

% mle: alpha sets the confidence interval e.g. 0.05 for 95%
[ests]=mle(dat,'nloglf',@neglog_find_sigma_ul,'start',[sigma p0],'lowerbound',xlb,'upperbound',xub);


function negll = neglog_find_sigma_ul(param,dat,cens,freq)
 
%	# recover the data
	n=length(dat)-2;
	zi=dat(1:n/2);
	si=dat(n/2+1:n);
    sigma_b=dat(n+1);
    gamma=dat(n+2);

%	# recover the parameters
	
	mu=gamma;
	sigma=param(1);
	p0=param(2);
    
    
    % NOW SIG_b ONLY FOR PORTION OF DE SINCE MINDOSE
    abOD=gamma*sigma_b;
    %si=abs(zi.*sqrt((si./zi).^2+abOD^2));
    si=sqrt(si.^2+abOD.^2);
    
%	# calculate sigma^2 + seld^2, mu0 and sigma0
	s2=sigma^2 + si.^2;
	sigma0=1./sqrt(1./sigma.^2 + 1./si.^2);
	mu0=(mu./sigma^2 + zi./si.^2)./(1./sigma^2 + 1./si.^2);


%	# calculate the log-likelihood

	logsqrt2pi=0.5*log(2*pi);
	res0=(gamma - mu0)./sigma0;
	res1=(gamma - mu)./sigma;
	lf1i=log(p0) - log(si) - 0.5*((zi-gamma)./si).^2   - logsqrt2pi;
	lf2i=log(1-p0) - 0.5*log(s2) - 0.5*(zi-mu).^2./s2  - logsqrt2pi;
	lf2i=lf2i + log(1-normcdf(res0)) - log(1-normcdf(res1));
	llik=log(exp(lf1i) + exp(lf2i));
    
%	# log-likelihood for censored data
%	#
%	#resc<- (mu - zc)/sigma
%	#llikc<- lnnprob.f(resc)

	negll=-sum(llik);
    if negll==-inf
        negll=0;
    else if isnan(negll)
            negll=0;
        end
    end
