function []=ProtocolB(u1,u0,v0,npart);
%estimate the parameters with data from Protocol B (single failues) with augmented data (i.e. remeasurement of some
%single failures) + follow-up on all passes with gold std.

%Arguments: 
%u1 is the number of components that pass and are nonconforming,   
%u0 is the number that pass and are conforming,   
%v0 is the total number that fail on both inspections
%npart is a vector of length r+1 that gives number of parts that gives the number of times we have s=0, 1, ..., r
%passes in the repeated measurement parts, note sum(npart)=n n should not be greater than v0!

%example call:  ProtocolB(2209,8,233,[7,13])

%translate data to other parameterization
g=u1;  u=u0+u1;  m=u1+u0+v0;  n=sum(npart);  r=length(npart)-1;

%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) followuplikelihood(para,m,n,r,u,npart,g),x0fix,Afix,bfix,Aeq,Beq,lbfix,ubfix,[],options);
%store the resulting estimates
resfix(1:3)=xfix;  
resfix(4)=xfix(1) * (1 - xfix(3)) / (xfix(1) * (1 - xfix(3)) + (1 - xfix(2)) * xfix(3));
resfix(5)=xfix(2) * xfix(3) / (xfix(2) * xfix(3) + (1 - xfix(1)) * (1 - xfix(3)));
    
%find asymptotics for fixed effect model
[sealpha,sebeta,sepic,sepbs,sepgn]=EinfofollowupFE(n,r,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(sepbs),')'])
disp(['theta1 = ',num2str(resfix(5)),' (',num2str(sepgn),')'])

%subroutines

function [loglik]=followuplikelihood(para,m,n,r,u,npart,g)
%returns the log-likelihood for fitting the fixed effects model to data
%in the NO gold standard case with the conditional sampling plan with baseline and followup
%m=baseline size, u=#passes in 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 are in followup with GS
%select sample for repeated measurements from rejects in baseline

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

%followup outcomes, 2 possibilities
lf = g*log((1-beta)*pic)+(u-g)*log(alpha*(1-pic));

%repeated measures log-likelihood for s passes
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));

%other baseline rejects not selected for repeated measurements
lextra=(m-u-n)*log(1-pip);
if (m-u-n)<0,  disp('TROUBLE'), end;

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



function [sealpha,sebeta,sepic,sepbs,sepgn]=EinfofollowupFE(n,r,m,alpha,beta,pic)
%find the asymptoic SEs for the latent class FE model with conditional sampling from first time fails
 
%baseline effect (remaining single failures)
u = m*((1-beta)*pic+alpha*(1-pic));  %expected number of passes parts
g = m*(1-beta)*pic;   %expected number of good parts in stream of passed parts
%expected information matrix for baseline
baselineI = -[-(m - u - n) * (-1 + pic) ^ 2 / (1 - (1 - beta) * pic - alpha * (1 - pic)) ^ 2 -(m - u - n) * (-1 + pic) / (1 - (1 - beta) * pic - alpha * (1 - pic)) ^ 2 * pic (m - u - n) / (1 - (1 - beta) * pic - alpha * (1 - pic)) - (m - u - n) * (-1 + pic) / (1 - (1 - beta) * pic - alpha * (1 - pic)) ^ 2 * (-1 + beta + alpha); -(m - u - n) * (-1 + pic) / (1 - (1 - beta) * pic - alpha * (1 - pic)) ^ 2 * pic -(m - u - n) * pic ^ 2 / (1 - (1 - beta) * pic - alpha * (1 - pic)) ^ 2 (m - u - n) / (1 - (1 - beta) * pic - alpha * (1 - pic)) - (m - u - n) * pic / (1 - (1 - beta) * pic - alpha * (1 - pic)) ^ 2 * (-1 + beta + alpha); (m - u - n) / (1 - (1 - beta) * pic - alpha * (1 - pic)) - (m - u - n) * (-1 + pic) / (1 - (1 - beta) * pic - alpha * (1 - pic)) ^ 2 * (-1 + beta + alpha) (m - u - n) / (1 - (1 - beta) * pic - alpha * (1 - pic)) - (m - u - n) * pic / (1 - (1 - beta) * pic - alpha * (1 - pic)) ^ 2 * (-1 + beta + alpha) -(m - u - n) * (-1 + beta + alpha) ^ 2 / (1 - (1 - beta) * pic - alpha * (1 - pic)) ^ 2;];
%also add followup information
followupI=[(u-g)/(alpha^2),0,0;0,g/((1-beta)^2),0;0,0,g/(pic^2)+(u-g)/((1-pic)^2)];

%add also information from repeated measurements
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+1))*pic+(alpha^s)*((1-alpha)^(r-s+1))*(1-pic))/((1-(1-beta))*pic+(1-alpha)*(1-pic));  %expected number of first time failures with s passes
  repeatedI = repeatedI - ESf*[(alpha ^ s * s ^ 2 / alpha ^ 2 * (1 - alpha) ^ (r - s + 1) * (1 - pic) - alpha ^ s * s / alpha ^ 2 * (1 - alpha) ^ (r - s + 1) * (1 - pic) - 2 * alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 1) * (r - s + 1) / (1 - alpha) * (1 - pic) + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (r - s + 1) ^ 2 / (1 - alpha) ^ 2 * (1 - pic) - alpha ^ s * (1 - alpha) ^ (r - s + 1) * (r - s + 1) / (1 - alpha) ^ 2 * (1 - pic)) / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) - (alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 1) * (1 - pic) - alpha ^ s * (1 - alpha) ^ (r - s + 1) * (r - s + 1) / (1 - alpha) * (1 - pic)) ^ 2 / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) ^ 2 -(alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 1) * (1 - pic) - alpha ^ s * (1 - alpha) ^ (r - s + 1) * (r - s + 1) / (1 - alpha) * (1 - pic)) / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) ^ 2 * (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 1) * pic + (1 - beta) ^ s * beta ^ (r - s + 1) * (r - s + 1) / beta * pic) (-alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 1) + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (r - s + 1) / (1 - alpha)) / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) - (alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 1) * (1 - pic) - alpha ^ s * (1 - alpha) ^ (r - s + 1) * (r - s + 1) / (1 - alpha) * (1 - pic)) / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) ^ 2 * ((1 - beta) ^ s * beta ^ (r - s + 1) - alpha ^ s * (1 - alpha) ^ (r - s + 1)); -(alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 1) * (1 - pic) - alpha ^ s * (1 - alpha) ^ (r - s + 1) * (r - s + 1) / (1 - alpha) * (1 - pic)) / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) ^ 2 * (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 1) * pic + (1 - beta) ^ s * beta ^ (r - s + 1) * (r - s + 1) / beta * pic) ((1 - beta) ^ s * s ^ 2 / (1 - beta) ^ 2 * beta ^ (r - s + 1) * pic - (1 - beta) ^ s * s / (1 - beta) ^ 2 * beta ^ (r - s + 1) * pic - 2 * (1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 1) * (r - s + 1) / beta * pic + (1 - beta) ^ s * beta ^ (r - s + 1) * (r - s + 1) ^ 2 / beta ^ 2 * pic - (1 - beta) ^ s * beta ^ (r - s + 1) * (r - s + 1) / beta ^ 2 * pic) / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) - (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 1) * pic + (1 - beta) ^ s * beta ^ (r - s + 1) * (r - s + 1) / beta * pic) ^ 2 / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) ^ 2 (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 1) + (1 - beta) ^ s * beta ^ (r - s + 1) * (r - s + 1) / beta) / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) - (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 1) * pic + (1 - beta) ^ s * beta ^ (r - s + 1) * (r - s + 1) / beta * pic) / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) ^ 2 * ((1 - beta) ^ s * beta ^ (r - s + 1) - alpha ^ s * (1 - alpha) ^ (r - s + 1)); (-alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 1) + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (r - s + 1) / (1 - alpha)) / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) - (alpha ^ s * s / alpha * (1 - alpha) ^ (r - s + 1) * (1 - pic) - alpha ^ s * (1 - alpha) ^ (r - s + 1) * (r - s + 1) / (1 - alpha) * (1 - pic)) / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) ^ 2 * ((1 - beta) ^ s * beta ^ (r - s + 1) - alpha ^ s * (1 - alpha) ^ (r - s + 1)) (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 1) + (1 - beta) ^ s * beta ^ (r - s + 1) * (r - s + 1) / beta) / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) - (-(1 - beta) ^ s * s / (1 - beta) * beta ^ (r - s + 1) * pic + (1 - beta) ^ s * beta ^ (r - s + 1) * (r - s + 1) / beta * pic) / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) ^ 2 * ((1 - beta) ^ s * beta ^ (r - s + 1) - alpha ^ s * (1 - alpha) ^ (r - s + 1)) -((1 - beta) ^ s * beta ^ (r - s + 1) - alpha ^ s * (1 - alpha) ^ (r - s + 1)) ^ 2 / ((1 - beta) ^ s * beta ^ (r - s + 1) * pic + alpha ^ s * (1 - alpha) ^ (r - s + 1) * (1 - pic)) ^ 2;];
end;
J=baselineI+repeatedI+followupI; 

V = inv(J);   %find the inverse

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

%sealpha=.0137; sebeta=.0090;  sepic=.0071;

%also determine the asymptotic standard errors for Pbadship and Pgoodnoship

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

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