function []=ProtocolA(u1,u0,v10,v11,v00);
%estimate the parameters with data from Protocol A (Double failues) with no augmented data (i.e. no remeasurement of
%double failures) follow-up on all one out of two passes with gold std.

%Arguments: 
%u1 is the number of components that pass initially and are nonconforming,   
%u0 is the number that pass initially and are conforming,   
%v11 is the number that pass the second inspection and are nonconforming,   
%v10 is the number that pass the second inspection and are conforming and finally,
%v00 is the number that fail on both inspections

%data from example in paper:  u1=1892; u0=23;  v10=26;  v11=256;  v00=253   %(total sample size =2450)
%example call:  ProtocolA(1892,23,26,256,253);

%translate data to other parameterization
g=u1;  g2=v11;  up=u0+u1;  upf=v10+v11;  m=u1+u0+v10+v11+v00;  
%e.g. data from the example in the paper up=1915;  upf=282;  g=1892;   g2=256; 

%constraints for the optimizer for fixed effects model
Aeq=[]; Beq=[];
Afix=[1,1,0]; bfix=[1]*0.9999;   %mua+mub<1
lbfix=[1,1,1]*0.0001;  ubfix=[.5,.5,1]*0.9999; %alpha and beta between 0 and 0.5, pic between 0 and 1
   
options = optimset('Display','off'); % Turn off Display, uncomment next line if display is desired
%options=[];
warning off

%set initial guesses for parameters - these values are arbitrary
alphahat=0.05;  betahat=0.10;  pichat=0.8;  %initial guesses
x0fix=[alphahat,betahat,pichat];  %initial guess

%fit FE model
xfix = fmincon(@(para) followup2likelihood(para,m,0,0,up,upf,0,g,g2),x0fix,Afix,bfix,Aeq,Beq,lbfix,ubfix,[],options);
%store the resulting estimates
resfix(1:3)=xfix;    
%also determine the estimates for Pbadship and Pgoodnoship (these are functions of the others)
resfix(4)=(1 - (1 - xfix(1)) ^ 2) * (1 - xfix(3)) / (1 - (1 - xfix(1)) ^ 2 * (1 - xfix(3)) - xfix(2) ^ 2 * xfix(3));
resfix(5)=xfix(2) ^ 2 * xfix(3) / ((1 - xfix(1)) ^ 2 * (1 - xfix(3)) + xfix(2) ^ 2 * xfix(3));
%find asymptotics for fixed effect model
[sealpha,sebeta,sepic,sePbad,sePgood]=Einfofollowup2FE(0,0,m,resfix(1),resfix(2),resfix(3));

disp('Parmeter Estimates and Approximate Standard Errors in brackets:')
disp(['alpha = ',num2str(resfix(1)),' (',num2str(sealpha),')'])
disp(['beta = ',num2str(resfix(2)),' (',num2str(sebeta),')'])
disp(['pi = ',num2str(resfix(3)),' (',num2str(sepic),')'])
disp(['theta0 = ',num2str(resfix(4)),' (',num2str(sePbad),')'])
disp(['theta1 = ',num2str(resfix(5)),' (',num2str(sePgood),')'])

%subroutines

function [loglik]=followup2likelihood(para,m,n,r,u,u2,npart,g,g2)
%returns the log-likelihood for fitting the fixed effects model to data
%use conditional sampling plan with baseline and followup
%RIM procedure: any unit that passes first time or second time is sent to customer where status is verified
%m=baseline size, u=#passes in first time baseline, u2=#passes in second time baseline
%n is number of parts in repeated meas study, r is the number of repeated measurements
%npart gives the number of times we have s=0, 1, ..., r passes in the repeated measurements part
%number of followup equals u, g gives the number of good parts in followup
%para = (alpha,beta,pic)
%assumes all parts that pass either first or second time are in followup with GS, g and g2 give number of good parts
%select sample for repeated measurements from twice rejects in baseline

alpha=para(1);  beta=para(2);  pic=para(3);
%pip=(1-beta)*pic+alpha*(1-pic);  %Pr(pass)

%followup outcomes first time passes, 2 possibilities
lf = g*log((1-beta)*pic)+(u-g)*log(alpha*(1-pic));
%same thing for second time passes
lf2 = g2*log(beta*(1-beta)*pic)+(u2-g2)*log(alpha*(1-alpha)*(1-pic));

%repeated measures log-likelihood for s passes on double rejects
s=0:r;  %all possible values for number of passes
%lsf=log((1-beta).^s.*beta.^(r-s+1)*pic+alpha.^s.*(1-alpha).^(r-s+1).*(1-pic));
lsf=log((1-beta).^s.*beta.^(r-s+2)*pic+alpha.^s.*(1-alpha).^(r-s+2).*(1-pic));

%other baseline double rejects not selected for repeated measurements
lextra=(m-u-u2-n)*log(beta^2*pic+(1-alpha)^2*(1-pic));

lprob=sum(lsf.*npart)+lf+lf2+lextra;   %log-likelihood
loglik=-lprob;  %take negative so we can minimize log-likelihood!


function [sealpha,sebeta,sepic,sepbs,sepgn]=Einfofollowup2FE(n,r,m,alpha,beta,pic)
%find the asymptoic SEs for the latent class FE model with conditional sampling from fails
 

u = m*((1-beta)*pic+alpha*(1-pic));  %expected number of first time passes=d parts
g = m*(1-beta)*pic;   %expected number of good parts in stream of once passed parts
u2 = m*((1-beta)*beta*pic+alpha*(1-alpha)*(1-pic));   %expected number of second time passes
g2 = m*beta*(1-beta)*pic;   %expected number of good parts in second time passes
%%baseline effect (remaining double failures), expected information matrix for baseline
baselineI = -[2 * (m - u - u2 - n) * (1 - pic) / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) - 4 * (m - u - u2 - n) * (1 - alpha) ^ 2 * (1 - pic) ^ 2 / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) ^ 2 4 * (m - u - u2 - n) * (1 - alpha) * (1 - pic) / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) ^ 2 * beta * pic 2 * (m - u - u2 - n) * (1 - alpha) / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) + 2 * (m - u - u2 - n) * (1 - alpha) * (1 - pic) / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) ^ 2 * (beta ^ 2 - (1 - alpha) ^ 2); 4 * (m - u - u2 - n) * (1 - alpha) * (1 - pic) / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) ^ 2 * beta * pic 2 * (m - u - u2 - n) * pic / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) - 4 * (m - u - u2 - n) * beta ^ 2 * pic ^ 2 / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) ^ 2 2 * (m - u - u2 - n) * beta / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) - 2 * (m - u - u2 - n) * beta * pic / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) ^ 2 * (beta ^ 2 - (1 - alpha) ^ 2); 2 * (m - u - u2 - n) * (1 - alpha) / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) + 2 * (m - u - u2 - n) * (1 - alpha) * (1 - pic) / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) ^ 2 * (beta ^ 2 - (1 - alpha) ^ 2) 2 * (m - u - u2 - n) * beta / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) - 2 * (m - u - u2 - n) * beta * pic / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) ^ 2 * (beta ^ 2 - (1 - alpha) ^ 2) -(m - u - u2 - n) * (beta ^ 2 - (1 - alpha) ^ 2) ^ 2 / (beta ^ 2 * pic + (1 - alpha) ^ 2 * (1 - pic)) ^ 2;];

%also add followup information on first time passes
followupI=[(u-g)/(alpha^2),0,0;0,g/((1-beta)^2),0;0,0,g/(pic^2)+(u-g)/((1-pic)^2)];
followup2I=-[(u2 - g2) * (-2 + 2 * pic) / alpha / (1 - alpha) / (1 - pic) - (u2 - g2) * ((1 - alpha) * (1 - pic) - alpha * (1 - pic)) / alpha ^ 2 / (1 - alpha) / (1 - pic) + (u2 - g2) * ((1 - alpha) * (1 - pic) - alpha * (1 - pic)) / alpha / (1 - alpha) ^ 2 / (1 - pic) 0 (u2 - g2) * (-1 + 2 * alpha) / alpha / (1 - alpha) / (1 - pic) + (u2 - g2) * ((1 - alpha) * (1 - pic) - alpha * (1 - pic)) / alpha / (1 - alpha) / (1 - pic) ^ 2; 0 -2 * g2 / beta / (1 - beta) - g2 * ((1 - beta) * pic - beta * pic) / beta ^ 2 / (1 - beta) / pic + g2 * ((1 - beta) * pic - beta * pic) / beta / (1 - beta) ^ 2 / pic g2 * (1 - 2 * beta) / beta / (1 - beta) / pic - g2 * ((1 - beta) * pic - beta * pic) / beta / (1 - beta) / pic ^ 2; (u2 - g2) * (-1 + 2 * alpha) / alpha / (1 - alpha) / (1 - pic) + (u2 - g2) * ((1 - alpha) * (1 - pic) - alpha * (1 - pic)) / alpha / (1 - alpha) / (1 - pic) ^ 2 g2 * (1 - 2 * beta) / beta / (1 - beta) / pic - g2 * ((1 - beta) * pic - beta * pic) / beta / (1 - beta) / pic ^ 2 -g2 / pic ^ 2 - (u2 - g2) / (1 - pic) ^ 2;];

%add also information from repeated measurements from double failures
repeatedI=zeros(3,3); 
for s=0:r,   %add up over all possible numbers of passes for the n parts repeatedly measured
  ESf = n*nchoosek(r,s)*(((1-beta)^s)*(beta^(r-s+2))*pic+(alpha^s)*((1-alpha)^(r-s+2))*(1-pic))/(beta^2*pic+(1-alpha)^2*(1-pic));  %expected number of double failures with s passes
  repeatedI = repeatedI - ESf*[(alpha ^ s * s ^ 2 / alpha ^ 2 * (1 - alpha) ^ (r - s + 2) * (1 - pic) - alpha ^ s * s / alpha ^ 2 * (1 - alpha) ^ (r - s + 2) * (1 - pic) - 2 * alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 2) * (r - s + 2) / (1 - alpha) * (1 - pic) + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (r - s + 2) ^ 2 / (1 - alpha) ^ 2 * (1 - pic) - alpha ^ s * (1 - alpha) ^ (r - s + 2) * (r - s + 2) / (1 - alpha) ^ 2 * (1 - pic)) / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) - (alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 2) * (1 - pic) - alpha ^ s * (1 - alpha) ^ (r - s + 2) * (r - s + 2) / (1 - alpha) * (1 - pic)) ^ 2 / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) ^ 2 -(alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 2) * (1 - pic) - alpha ^ s * (1 - alpha) ^ (r - s + 2) * (r - s + 2) / (1 - alpha) * (1 - pic)) / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) ^ 2 * (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 2) * pic + (1 - beta) ^ s * beta ^ (r - s + 2) * (r - s + 2) / beta * pic) (-alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 2) + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (r - s + 2) / (1 - alpha)) / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) - (alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 2) * (1 - pic) - alpha ^ s * (1 - alpha) ^ (r - s + 2) * (r - s + 2) / (1 - alpha) * (1 - pic)) / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) ^ 2 * ((1 - beta) ^ s * beta ^ (r - s + 2) - alpha ^ s * (1 - alpha) ^ (r - s + 2)); -(alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 2) * (1 - pic) - alpha ^ s * (1 - alpha) ^ (r - s + 2) * (r - s + 2) / (1 - alpha) * (1 - pic)) / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) ^ 2 * (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 2) * pic + (1 - beta) ^ s * beta ^ (r - s + 2) * (r - s + 2) / beta * pic) ((1 - beta) ^ s * s ^ 2 / (1 - beta) ^ 2 * beta ^ (r - s + 2) * pic - (1 - beta) ^ s * s / (1 - beta) ^ 2 * beta ^ (r - s + 2) * pic - 2 * (1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 2) * (r - s + 2) / beta * pic + (1 - beta) ^ s * beta ^ (r - s + 2) * (r - s + 2) ^ 2 / beta ^ 2 * pic - (1 - beta) ^ s * beta ^ (r - s + 2) * (r - s + 2) / beta ^ 2 * pic) / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) - (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 2) * pic + (1 - beta) ^ s * beta ^ (r - s + 2) * (r - s + 2) / beta * pic) ^ 2 / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) ^ 2 (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 2) + (1 - beta) ^ s * beta ^ (r - s + 2) * (r - s + 2) / beta) / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) - (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 2) * pic + (1 - beta) ^ s * beta ^ (r - s + 2) * (r - s + 2) / beta * pic) / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) ^ 2 * ((1 - beta) ^ s * beta ^ (r - s + 2) - alpha ^ s * (1 - alpha) ^ (r - s + 2)); (-alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 2) + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (r - s + 2) / (1 - alpha)) / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) - (alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 2) * (1 - pic) - alpha ^ s * (1 - alpha) ^ (r - s + 2) * (r - s + 2) / (1 - alpha) * (1 - pic)) / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) ^ 2 * ((1 - beta) ^ s * beta ^ (r - s + 2) - alpha ^ s * (1 - alpha) ^ (r - s + 2)) (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 2) + (1 - beta) ^ s * beta ^ (r - s + 2) * (r - s + 2) / beta) / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) - (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 2) * pic + (1 - beta) ^ s * beta ^ (r - s + 2) * (r - s + 2) / beta * pic) / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) ^ 2 * ((1 - beta) ^ s * beta ^ (r - s + 2) - alpha ^ s * (1 - alpha) ^ (r - s + 2)) -((1 - beta) ^ s * beta ^ (r - s + 2) - alpha ^ s * (1 - alpha) ^ (r - s + 2)) ^ 2 / ((1 - beta) ^ s * beta ^ (r - s + 2) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 2) * (1 - pic)) ^ 2;];
end;
J=baselineI+repeatedI+followupI+followup2I; 

V = inv(J);   %find the inverse

sealpha = sqrt(V(1,1));
sebeta = sqrt(V(2,2));
sepic = sqrt(V(3,3));

%also determine the asymptotic standard errors for Pbadship and Pgoodnoship

GPbs=[(2 - 2 * alpha) * (1 - pic) / (1 - (1 - alpha) ^ 2 * (1 - pic) - beta ^ 2 * pic) - 2 * (1 - (1 - alpha) ^ 2) * (1 - pic) ^ 2 / (1 - (1 - alpha) ^ 2 * (1 - pic) - beta ^ 2 * pic) ^ 2 * (1 - alpha) 2 * (1 - (1 - alpha) ^ 2) * (1 - pic) / (1 - (1 - alpha) ^ 2 * (1 - pic) - beta ^ 2 * pic) ^ 2 * beta * pic -(1 - (1 - alpha) ^ 2) / (1 - (1 - alpha) ^ 2 * (1 - pic) - beta ^ 2 * pic) - (1 - (1 - alpha) ^ 2) * (1 - pic) / (1 - (1 - alpha) ^ 2 * (1 - pic) - beta ^ 2 * pic) ^ 2 * ((1 - alpha) ^ 2 - beta ^ 2)];
GPgn=[2 * beta ^ 2 * pic / ((1 - alpha) ^ 2 * (1 - pic) + beta ^ 2 * pic) ^ 2 * (1 - alpha) * (1 - pic) 2 * beta * pic / ((1 - alpha) ^ 2 * (1 - pic) + beta ^ 2 * pic) - 2 * beta ^ 3 * pic ^ 2 / ((1 - alpha) ^ 2 * (1 - pic) + beta ^ 2 * pic) ^ 2 beta ^ 2 / ((1 - alpha) ^ 2 * (1 - pic) + beta ^ 2 * pic) - beta ^ 2 * pic / ((1 - alpha) ^ 2 * (1 - pic) + beta ^ 2 * pic) ^ 2 * (-(1 - alpha) ^ 2 + beta ^ 2)];

sepbs=sqrt(GPbs*V*GPbs');
sepgn=sqrt(GPgn*V*GPgn');