function MAM3ul=getmam3ul(de,sigma_b,start)
% This is an unlogged version of the MAM3. calls neglogMAM3, which was
% translated from Arnold's s-plus version.

% OUTPUT is four numbers 
% 1) minimum age
% 2) standard error (CI estimation with free gamma only)
% 3) s.d. of the half-normal distribution of poorly bleached aliquots
% 4) proportion of aliquots at the minimum age
% ENSURE THAT THESE NUMBERS MAKE SENSE, OTHERWISE CHECK UPPER/LOWER BOUNDS

% INPUTS 
% de        % Name of the data in the workspace - Requires 2 cols [de se]
% sigma_b   % Specified overdispersion (e.g. 0.10 for 10%)
% start     % Starting value for minimum age

% ----------------------------------------------------------------------

% 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=[-2 0.01 0.01];
xub=[100 100 1];

%combine errors and format data
%comb_errs=sqrt(de(:,2).^2+sigma_b^2); %this is what Arnold has
comb_errs=abs(de(:,1).*sqrt((de(:,2)./de(:,1)).^2+sigma_b^2)); % This seems more likely

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

% mle: alpha sets the confidence interval e.g. 0.05 for 95%
[ests ci68]=mle(dat,'nloglf',@neglogMAM3,'start',[gamma sigma p0],'alpha',0.32,'lowerbound',xlb,'upperbound',xub);
%pal=ests(1);
p0end=ests(3);
sig=ests(2);
% uncomment this to show the negative log likelihood
%negloglike=neglogMAM3(ests,dat)

% Matlab's ci estimates, too big because too many free parameters?
%err=mean([ests(1)-ci68(1,1),ci68(2,1)-ests(1)]);
%paleo_se_sig_p0=[ests(1) err sig p0end]

% Estimate error by fixing all parameters other than gamma
dat2=[dat; sig; p0end]; % attatch other params as data
[Fests Fci68]=mle(dat2,'nloglf',@neglogMAM1,'start',ests(1),'alpha',0.32); % call 1 comp mam
err=mean([ests(1)-Fci68(1,1),Fci68(2,1)-ests(1)]);
MAM3ul=[Fests err sig p0end];

function negll = neglogMAM3(param,dat,cens,freq)

%   # this calculates the negative of the log likelihood of the  
%	# data (dat,rep) for a given set of parameters (param)
%	#
%	# param is the vector of model parameters: 
%	#        (gamma, mu, sigma, p0)
%	# datmat is a nx2 matrix of data: ld, seld (including sigma_b)
	
       
%	# recover the data
	n=length(dat);
	zi=dat(1:n/2);
	si=dat(n/2+1:n);


%	# recover the parameters
	gamma=param(1);
	mu=gamma;
	sigma=param(2);
	p0=param(3);

%	# 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

    function negll = neglogMAM1(param,dat,cens,freq)

%   # this calculates the negative of the log likelihood of the  
%	# data (dat,rep) for a given set of parameters (param)
%	#
%	# param is the vector of model parameters: 
%	#        (gamma, mu, sigma, p0)
%	# datmat is a nx2 matrix of data: ld, seld (including sigma_b)
	
       
%	# recover the data
	n=length(dat)-2;
	zi=dat(1:n/2);
	si=dat(n/2+1:n);


%	# recover the parameters
	gamma=param(1);
	mu=gamma;
	sigma=dat(n+1);
	p0=dat(n+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
