function [res,out]=MLEresults(k,m,n,kab,typeflag,model,SPdata,APdata)
%find the MLEs given some data
%k = # parts
%m = # operators
%n = # measurements on each part
%kab is ka or kb = # parts measured once, must be a multiple of m
%typeflag = 1 for type A augment and 2 for type B augmented plan
%model = 0 for no interaction model, anything else for interaction model
%SPdata is a column vector ordered first by parts and then by operators
%APdata is a column vector also ordered first by parts and then by operators
%example call - corresponds to provided example
%[res]=MLEresults(6,3,2,24,1,1,[10.4; 9.6; 10.8; 10.8; 10.6; 10.4; 23.8; 22.5; 23.3; 21.7; 23.1; 23.3; 25.0; 25.3; 26.0; 25.4; 27.6; 28.0; 22.5; 22.3; 24.4; 24.1; 23.3; 23.3; 22.9; 22.6; 22.3; 22.2; 22.5; 23.1; 18.1; 17.5; 19.8; 20.6; 17.7; 19.6],[18.8; 12.0; 31.6; 16.3; 18.9; 18.2; 5.1; 18.0; 11.2; 23.8; 18.5; 15.0; 28.3; 26.4; 28.6; 11.9; 25.6; 22.3; 17.2; 27.1; 27.4; 23.4; 24.6; 24.5]);
eparts=kab; %number of parts in augmented part
if kab/m~=round(kab/m), disp('ERROR kab must be a multiple of m'), end;
if typeflag==1,
u=kab/m;
elseif typeflag==2,
u=kab;
else
u=0;
end;
N=m*(k*n+u); %N is the total number of measurements
Plan=[k,m,n]; if u>0, Plan=[Plan,u,typeflag]; end;
switch length(Plan),
case 3
disp(['Standard Plan with ',num2str(Plan(1)),' parts ',num2str(Plan(2)),' operators and ',num2str(Plan(3)),' repeated measurements'])
disp(['total number of measurements ',num2str(Plan(1)*Plan(2)*Plan(3))])
case 5
type=Plan(5);
if type==1,
disp(['Augmented Plan type A with ',num2str(Plan(1)),' parts ',num2str(Plan(2)),' operators and ',num2str(Plan(3)),' repeated measurements, plus ',num2str(Plan(2)*Plan(4)),' additional parts in the augmented plan'])
disp(['total number of measurements = ',num2str(Plan(1)*Plan(2)*Plan(3)+Plan(2)*Plan(4))])
elseif type==2
disp(['Augmented Plan type B with ',num2str(Plan(1)),' parts ',num2str(Plan(2)),' operators and ',num2str(Plan(3)),' repeated measurements, plus ',num2str(Plan(4)),' additional parts in the augmented plan'])
disp(['total number of measurements = ',num2str(Plan(1)*Plan(2)*Plan(3)+Plan(4)*Plan(2))])
else disp('Plan is an invlaid type of augmented plan')
end;
end;
options=optimset('Display','off');
options = optimset(options,'GradObj','on'); %use gradient information
%calculate some data summaries to use as initial guesses for MLE routine
%NOTE: in column vector SPdata operators go 1 to m over and over
%this is very simple and does not give good estimates
if model==0, sigmad2=0;
else sigmad2=0.1*var([SPdata;APdata]); end;
sigmam2=0.1*var([SPdata;APdata]); %partition total variation
sigmap2=0.9*var([SPdata;APdata]);
temp=reshape(SPdata,m*n,k); %each column is now the results from a single part
if n>1, muj=mean(reshape(mean(temp'),n,m))'; %operator mean estimates
else muj=mean(temp')'; %if n=1
end;
%maximize the likelihood directly to find the MLEs for mu, sigmam2, sigmad2 and sigmap2
%[x,fval,exitflag] = fminunc(f,[mujinit,log(sigmap2init),log(sigmam2init)],options); %use data summaries as initial values
c=[typeflag,k,m,n,u,SPdata',APdata']; %fixed arguments and data
if model==0 || (typeflag==0 && n==1), %use no interaction model, i.e. don't estimate sigmad2
f_nosd = @(x)loglik_nosd(x,c); %create anonymous function to allow passing of fixed parameters in c
[x,fval,exitflag] = fminunc(f_nosd,[muj',log(sigmap2),log(sigmam2+sigmad2)],options); %use true values as initial values
else
f = @(x)loglik(x,c); %create anonymous function to allow passing of fixed parameters in c
[x,fval,exitflag] = fminunc(f,[muj',log(sigmap2),log(sigmad2),log(sigmam2)],options); %use true values as initial values
end; %use ln scale in above optimization to avoid need for constraints
tries=1; maxtries=100;
while (exitflag==0) && (tries<maxtries), %didn't converge, try with different initial conditions
disp(['did not converge, try # ', num2str(tries)]), tries=tries+1;
if model==0 || (typeflag==0 && n==1),
[x,fval,exitflag] = fminunc(f_nosd,[muj',log(rand),log(rand)],options); %pick random initial values
else
[x,fval,exitflag] = fminunc(f,[muj',log(rand),log(rand),log(rand)],options); %pick random initial values
end;
end;
if tries==maxtries, disp('MAX TRIES REACHED'), end;
%translate MLEs given by vector x to metrics of interest
mujhat=x(1:m)';
if model==0 || (typeflag==0 && n==1), %use no interaction model, i.e. don't estimate sigmad2
sigmap2hat=exp(x(end-1));
sigmad2hat=0; %assume no interaction model
sigmam2hat=exp(x(end)); %undo log transformation
else
sigmap2hat=exp(x(end-2));
sigmad2hat=exp(x(end-1));
sigmam2hat=exp(x(end)); %undo log transformation
end;
sigmao2hat=mean((mujhat-mean(mujhat)).^2);
gammahat=sqrt((sigmao2hat+sigmad2hat+sigmam2hat)/(sigmap2hat+sigmao2hat+sigmad2hat+sigmam2hat));
%store results
res=[sqrt(sigmap2hat),sqrt(sigmao2hat),sqrt(sigmad2hat),sqrt(sigmam2hat),gammahat];
%determine standard errors using Fisher information matrix with estimates for sp2, etc.
%would perhaps be better to use the observed information instead
if u==0, %SP only
Jmat=J_kmn(k,m,nmax,sigmap2hat,sigmad2hat,sigmam2hat);
elseif typeflag==1, %AP A
Jmat=J_kmn(k,m,n,sigmap2hat,sigmad2hat,sigmam2hat)+J_typeA(m,u,sigmap2hat,sigmad2hat,sigmam2hat); %add the two information matrices for type A AP
else %AP B
Jmat=J_kmn(k,m,n,sigmap2hat,sigmad2hat,sigmam2hat)+J_kmn(u,m,1,sigmap2hat,sigmad2hat,sigmam2hat); %add the two information matrices
end;
varmat=inv(Jmat); %diagonals give variance for sp2, sd2 and sm2
SEmuj=sqrt(diag(varmat(1:m,1:m)));
SEsp=sqrt(varmat(m+1,m+1))/(2*sqrt(sigmap2hat)); %determine the SE for sp (not sp2)
D=cmatrix(mujhat,sigmap2hat,sigmad2hat,sigmam2hat);
temp=D*inv(Jmat)*D';
if sigmad2hat~=0,
SEso=sqrt(temp(2,2));
SEsd=sqrt(temp(3,3));
SEsm=sqrt(temp(4,4));
SEgamma=sqrt(temp(1,1));
else
if m>=2,
SEso=sqrt(temp(2,2));
SEsd=0;
SEsm=sqrt(temp(3,3));
SEgamma=sqrt(temp(1,1));
else
SEso=0;
SEsd=0;
SEsm=sqrt(temp(2,2));
SEgamma=sqrt(temp(1,1));
end;
end;
out=[mujhat',sqrt(sigmap2hat),sqrt(sigmao2hat),sqrt(sigmad2hat),sqrt(sigmam2hat),gammahat; SEmuj',SEsp,SEso,SEsd,SEsm,SEgamma];
%summarize results
disp(' operator means')
disp(['MLEs = ',num2str(mujhat')])
disp(['MLE SEs = ',num2str(SEmuj')])
disp(' sp so sd sm gamma')
disp(['MLEs = ',num2str([sqrt(sigmap2hat),sqrt(sigmao2hat),sqrt(sigmad2hat),sqrt(sigmam2hat),gammahat])])
disp(['MLE SEs = ',num2str([SEsp,SEso,SEsd,SEsm,SEgamma])])
function [ll,G]=loglik(x,c)
%return the negative log likelihood function and the corresponding gradients
typeflag=c(1); k=c(2); m=c(3); n=c(4); u=c(5);
SPdata=c(6:6+n*k*m-1)'; %change back to column vectors
APdata=c(6+n*k*m:end)';
muj=x(1:m)';
sp2=exp(x(end-2)); sd2=exp(x(end-1)); sm2=exp(x(end)); %undo log transform
[ll,G]=likelihood_and_gradients(typeflag,k,m,n,u,SPdata,APdata,muj,sp2,sd2,sm2);
function [ll,G]=loglik_nosd(x,c)
%return the negative log likelihood function and the corresponding
%gradients, no interaction model
typeflag=c(1); k=c(2); m=c(3); n=c(4); u=c(5);
SPdata=c(6:6+n*k*m-1)'; %change back to column vectors
APdata=c(6+n*k*m:end)';
muj=x(1:m)';
sp2=exp(x(end-1)); sd2=0; sm2=exp(x(end)); %undo log transform
[ll,G]=likelihood_and_gradients(typeflag,k,m,n,u,SPdata,APdata,muj,sp2,sd2,sm2);
G(m+2)=[]; %remove gradient for sd2
function [ll,G]=likelihood_and_gradients(typeflag,k,m,n,u,SPdata,APdata,muj,sp2,sd2,sm2)
%summarize the Stage 1 data i.e. calculate the sum of squares in the log-likelihood
%muv=kron(ones(k,1),kron(muj,ones(n,1))) %order by parts then operators too slow, kron function is slow
muv=repmat(kron(muj,ones(n,1)),k,1);
SPdataminusmeans=SPdata-muv; %subtract off operator means
%calculate functions of the data
%A=SPdataminusmeans'*eye(k*m*n)*SPdataminusmeans %A is sum (yijl-muj)
A=sum(SPdataminusmeans.*SPdataminusmeans);
temp=reshape(SPdataminusmeans,n,k*m); %each column corresponds to different part and operator combination
B=sum(sum(temp).^2);
%B=SPdataminusmeans'*kron(eye(k*m),ones(n,n))*SPdataminusmeans
temp=reshape(SPdataminusmeans,m*n,k); %each column corresponds to a different part
C=sum(sum(temp).^2);
%C=SPdataminusmeans'*kron(eye(k),ones(m*n,m*n))*SPdataminusmeans
%likelihood from Stage 1
ll1= -0.1e1 / sm2 * A / 0.2e1 + sd2 / sm2 / (sm2 + n * sd2) * B / 0.2e1 + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) * C / 0.2e1 - k * log((sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m)) / 0.2e1;
%determine gradients from Stage 1
G=zeros(m+3,1);
%temp=reshape(SPdataminusmeans,m,k*n); %each row is now the results from a single operator
A2=zeros(m,1); %m different values one for each operator
for j=1:m, for l=1:n, %find the sum by operator
A2(j)=A2(j)+sum(SPdataminusmeans(l+(j-1)*n:m*n:end));
end; end;
B3=sum(SPdataminusmeans); %the same for all operators - since we have a balanced design
b1=1/sm2;
b2=-sd2 / sm2 / (sm2 + n * sd2);
b3=-sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2);
G(1:m)=b1*A2+b2*A2*n+b3*B3*n; %gradient wrt to muj
dldsm2 = 0.1e1 / sm2 ^ 2 * A / 0.2e1 - sd2 / sm2 ^ 2 / (sm2 + n * sd2) * B / 0.2e1 - sd2 / sm2 / (sm2 + n * sd2) ^ 2 * B / 0.2e1 - sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) * C / 0.2e1 - sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C / 0.2e1 - k * ((sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) / 0.2e1;
dldsd2 = 0.1e1 / sm2 / (sm2 + n * sd2) * B / 0.2e1 - sd2 / sm2 / (sm2 + n * sd2) ^ 2 * B * n / 0.2e1 - sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) * C * n / 0.2e1 - sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * n / 0.2e1 - k * (n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) / 0.2e1;
dldsp2 = 0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) * C / 0.2e1 - sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n / 0.2e1 - k * m * n / (sm2 + n * sd2 + m * n * sp2) / 0.2e1;
%add results for the augmented Stage 2 data
ll2=0; dlAugdmu = zeros(m,1); dlAugdsp2 = 0; dlAugdsd2 = 0; dlAugdsm2 = 0; %assume no augmented data
if (u>0), %measurements in Stage 2
muv=reshape(muj*ones(1,u),u*m,1); %estimated operator effects, one after each other and repeated
APdataminusmeans=APdata-muv;
A=sum(APdataminusmeans.^2); %Sum #1 in AP A Likelihood Expression
if typeflag==1, %AP A
ll2 = -A / (sp2 + sd2 + sm2) / 0.2e1 - u*m * log(sp2 + sd2 + sm2) / 0.2e1;
%find gradients for likelihood of augmented data
temp2=reshape(APdataminusmeans,m,u); %each row is now the results from a single operator
csum=sum(temp2')'; %m different values one for each operator
dlAugdmu = csum/(sp2+sd2+sm2);
dlAugdsp2 = A / (sp2 + sd2 + sm2) ^ 2 / 0.2e1 - u*m / (sp2 + sd2 + sm2) / 0.2e1;
dlAugdsd2 = A / (sp2 + sd2 + sm2) ^ 2 / 0.2e1 - u*m / (sp2 + sd2 + sm2) / 0.2e1;
dlAugdsm2 = A / (sp2 + sd2 + sm2) ^ 2 / 0.2e1 - u*m / (sp2 + sd2 + sm2) / 0.2e1;
elseif typeflag==2, %AP B, extra likelihood a special case of SP likelihood
B=A; %Sum #2 in AP B Likelihood Expression combines with first sum
n=1; k=u;
b2=-sd2 / sm2 / (sm2 + n * sd2); %b1 doesn't depend on n so we use the existing value
b3=-sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2);
C=APdataminusmeans'*kron(eye(k),ones(m*n,m*n))*APdataminusmeans;
ll2=-0.1e1 / sm2 * A / 0.2e1 + sd2 / sm2 / (sm2 + sd2) * B / 0.2e1 + sp2 / (sm2 + sd2) / (sm2 + sd2 + m * sp2) * C / 0.2e1 - u * log((sm2 + sd2 + m * sp2) * (sm2 + sd2) ^ (m - 0.1e1)) / 0.2e1;
%gradients for AP type B is the same format as for SP data with n=1 and k=u and APdata rather than SPdata
A2=zeros(m,1); %m different values one for each operator
for j=1:m, %find the sum by operator
A2(j)=sum(APdataminusmeans(j:m:end));
end;
B3=sum(APdataminusmeans); %the same for all operators - since we have a balanced design
dlAugdmu = b1*A2+b2*A2*n+b3*B3*n; %gradient wrt to muj
dlAugdsm2 = 0.1e1 / sm2 ^ 2 * A / 0.2e1 - sd2 / sm2 ^ 2 / (sm2 + n * sd2) * B / 0.2e1 - sd2 / sm2 / (sm2 + n * sd2) ^ 2 * B / 0.2e1 - sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) * C / 0.2e1 - sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C / 0.2e1 - k * ((sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) / 0.2e1;
dlAugdsd2 = 0.1e1 / sm2 / (sm2 + n * sd2) * B / 0.2e1 - sd2 / sm2 / (sm2 + n * sd2) ^ 2 * B * n / 0.2e1 - sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) * C * n / 0.2e1 - sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * n / 0.2e1 - k * (n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) / 0.2e1;
dlAugdsp2 = 0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) * C / 0.2e1 - sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n / 0.2e1 - k * m * n / (sm2 + n * sd2 + m * n * sp2) / 0.2e1;
end
end
%update the likelihood and gradients to reflect the augmented data
ll=-(ll1+ll2); %negative so that minimization is finding the MLEs
%update the gradients with terms from the (Stage 2) augmented data
G(1:m) = G(1:m) + dlAugdmu;
dldsm2 = dldsm2+dlAugdsm2;
dldsd2 = dldsd2+dlAugdsd2;
dldsp2 = dldsp2+dlAugdsp2;
%need to translate from sigma2 scale to log(sigma2) scale - apply chain rule, also see Maple results
G(m+1)=dldsp2*sp2; %sp2
G(m+2)=dldsd2*sd2; %sd2
G(m+3)=dldsm2*sm2; %sm2
G=-G; %negative so that minimization is finding the MLEs
function [Jmat]=J_kmn(k,m,n,sp2,sd2,sm2)
%find the information matrix for the SP given k, m, n
if sd2 ~=0, Jmat=zeros(m+3,m+3); else Jmat=zeros(m+2,m+2); end;
%Note that the information matrix Jmat is always block diagonal with terms
%related to the mujs and another block for the rest
b1 = 1/sm2;
b2 = -sd2 / (sm2 * (sm2 + n * sd2));
b3 = -sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2);
%setup top bit related to the mus
Jmat(1:m,1:m)=-k*n*n*b3;
for i=1:m, %replace the diagonal
Jmat(i,i)=-k*n*(b1+n*(b2+b3));
end;
%expected vlaues of the various sums of squares
A=k*m*n*(sp2+sd2+sm2);
B=k*m*n*(n*sp2+n*sd2+sm2);
C=k*m*n*(m*n*sp2+n*sd2+sm2);
%now look at bottom right 3x3 part related to sp2, sd2 and sm2 in that order!
if sd2~=0,
Jmat(m+1,m+1) = -0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m ^ 2 * n ^ 2 + k * m ^ 2 * n ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
Jmat(m+1,m+2) = -0.1e1 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) * C * n / 0.2e1 - 0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * n / 0.2e1 + sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n ^ 2 / 0.2e1 + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m * n ^ 2 + k * m * n ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
Jmat(m+2,m+1) = Jmat(m+1,m+2);
Jmat(m+1,m+3) = -0.1e1 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) * C / 0.2e1 - 0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C / 0.2e1 + sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n / 0.2e1 + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m * n + k * m * n / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
Jmat(m+3,m+1) = Jmat(m+1,m+3);
Jmat(m+2,m+2) = -(1 / sm2 / (sm2 + n * sd2) ^ 2 * B * n) + (sd2 / sm2 / (sm2 + n * sd2) ^ 3 * B * n ^ 2) + (sp2 / (sm2 + n * sd2) ^ 3 / (sm2 + n * sd2 + m * n * sp2) * C * n ^ 2) + (sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * n ^ 2) + (sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * n ^ 2) - (k * (2 * n ^ 2 * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) ^ 2 * n ^ 2 / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) * n ^ 2 / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * (n * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) ^ 2 / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * n) / 0.2e1 + (k * (n * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m - 1) * n / (sm2 + n * sd2)) / 0.2e1;
Jmat(m+2,m+3) = -0.1e1 / sm2 ^ 2 / (sm2 + n * sd2) * B / 0.2e1 - 0.1e1 / sm2 / (sm2 + n * sd2) ^ 2 * B / 0.2e1 + sd2 / sm2 ^ 2 / (sm2 + n * sd2) ^ 2 * B * n / 0.2e1 + sd2 / sm2 / (sm2 + n * sd2) ^ 3 * B * n + sp2 / (sm2 + n * sd2) ^ 3 / (sm2 + n * sd2 + m * n * sp2) * C * n + sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * n + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * n - k * (0.2e1 * n * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) ^ 2 / (sm2 + n * sd2) ^ 2 * n * sm2 ^ (m * n - m) - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) / 0.2e1 + k * (n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) ^ 2 / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) / 0.2e1 + k * (n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) * (m - 0.1e1) / (sm2 + n * sd2) / 0.2e1 + k * (n * (sm2 + n * sd2) ^ (m - 0.1e1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 0.1e1) * (m - 0.1e1) * n / (sm2 + n * sd2) * sm2 ^ (m * n - m)) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 0.1e1) / sm2 ^ (m * n - m) * (m * n - m) / sm2 / 0.2e1;
Jmat(m+3,m+2) = Jmat(m+2,m+3);
Jmat(m+3,m+3) = -(1 / sm2 ^ 3 * A) + (sd2 / sm2 ^ 3 / (sm2 + n * sd2) * B) + (sd2 / sm2 ^ 2 / (sm2 + n * sd2) ^ 2 * B) + (sd2 / sm2 / (sm2 + n * sd2) ^ 3 * B) + (sp2 / (sm2 + n * sd2) ^ 3 / (sm2 + n * sd2 + m * n * sp2) * C) + (sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C) + (sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C) - (k * (2 * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + 2 * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) ^ 2 / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) + 2 * (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) ^ 2 / sm2 ^ 2 - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 ^ 2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m - 1) / (sm2 + n * sd2)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m * n - m) / sm2) / 0.2e1;
else %remove bits related to sd2
Jmat(m+1,m+1) = -0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m ^ 2 * n ^ 2 + k * m ^ 2 * n ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
Jmat(m+1,m+2) = -0.1e1 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) * C / 0.2e1 - 0.1e1 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C / 0.2e1 + sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C * m * n / 0.2e1 + sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C * m * n + k * m * n / (sm2 + n * sd2 + m * n * sp2) ^ 2 / 0.2e1;
Jmat(m+2,m+1) = Jmat(m+1,m+2);
Jmat(m+2,m+2) = -(1 / sm2 ^ 3 * A) + (sd2 / sm2 ^ 3 / (sm2 + n * sd2) * B) + (sd2 / sm2 ^ 2 / (sm2 + n * sd2) ^ 2 * B) + (sd2 / sm2 / (sm2 + n * sd2) ^ 3 * B) + (sp2 / (sm2 + n * sd2) ^ 3 / (sm2 + n * sd2 + m * n * sp2) * C) + (sp2 / (sm2 + n * sd2) ^ 2 / (sm2 + n * sd2 + m * n * sp2) ^ 2 * C) + (sp2 / (sm2 + n * sd2) / (sm2 + n * sd2 + m * n * sp2) ^ 3 * C) - (k * (2 * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + 2 * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) ^ 2 / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) ^ 2 * sm2 ^ (m * n - m) + 2 * (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) * (m * n - m) / sm2 + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) ^ 2 / sm2 ^ 2 - (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2 ^ 2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) ^ 2 / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m - 1) / (sm2 + n * sd2)) / 0.2e1 + (k * ((sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * (m - 1) / (sm2 + n * sd2) * sm2 ^ (m * n - m) + (sm2 + n * sd2 + m * n * sp2) * (sm2 + n * sd2) ^ (m - 1) * sm2 ^ (m * n - m) * (m * n - m) / sm2) / (sm2 + n * sd2 + m * n * sp2) / (sm2 + n * sd2) ^ (m - 1) / sm2 ^ (m * n - m) * (m * n - m) / sm2) / 0.2e1;
end;
%change the sign since we want the negative of the expected value instead of the expected value
Jmat=-Jmat;
function [Jmat]=J_typeA(m,u,sp2,sd2,sm2)
%find the special extra information matrix for the type A AP given m and u
%with AP plan A m operators each measure u *different* parts once each
%total number of parts in extra AP bit is m*u
if sd2 ~=0, Jmat=zeros(m+3,m+3);
else Jmat=zeros(m+2,m+2); end;
%setup top bit related to the mus
Jmat(1:m,1:m)=-u/(sp2+sd2+sm2)*eye(m);
%n=1; k=u; %special case
%now look at bottom right 3x3 part related to sp2, sd2, sm2
A=m*u*(sp2+sd2+sm2);
ka=m*u;
if sd2~=0,
Jmat(m+1:m+3,m+1:m+3) = -A / (sp2 + sd2 + sm2) ^ 3 + ka / (sp2 + sd2 + sm2) ^ 2 / 0.2e1; %3x3 matrix all the same
%off diagonal elements also all the same
Jmat(m+1,m+2)=-A / (sp2 + sd2 + sm2) ^ 3 + ka / (sp2 + sd2 + sm2) ^ 2 / 0.2e1;
Jmat(m+1,m+3)=Jmat(m+1,m+2);
Jmat(m+2,m+3)=Jmat(m+1,m+2);
Jmat(m+2,m+1)=Jmat(m+1,m+2);
Jmat(m+3,m+2)=Jmat(m+1,m+2);
Jmat(m+3,m+1)=Jmat(m+1,m+2);
else %remove bits related to sd2
Jmat(m+1:m+2,m+1:m+2) = -A / (sp2 + sd2 + sm2) ^ 3 + ka / (sp2 + sd2 + sm2) ^ 2 / 0.2e1; %3x3 matrix all the same
%off diagonal elements also all the same
Jmat(m+1,m+2)=-A / (sp2 + sd2 + sm2) ^ 3 + ka / (sp2 + sd2 + sm2) ^ 2 / 0.2e1;
Jmat(m+2,m+1)=Jmat(m+1,m+2);
end;
%change the sign since we want the negative of the expected value instead of the expected value
Jmat=-Jmat;
function [mat]=cmatrix(mu,sp2,sd2,sm2)
%returns the change of variables matrix based on the partial derivatives
%go from muis sp2, sd2 and sm2 (note variance scale) to gamma, so, sd, sm in that order
%skip row and column for sd if we assume the no interaction model
%also remove the so column if there is only one operat, i.e. m=1
m=length(mu); %number of operators
if sd2~=0, %interaction model, this doesn't make sense if there is only one operator
mat=zeros(4,m+3); %four rows
if m==1, disp('WARNING, using the interaction model when there is only one operator does not make sense, select sd2=0'), end;
mat(3,m+2)=1/(2*sqrt(sd2)); %only element in third row for sd
mat(4,m+3)=1/(2*sqrt(sm2)); %only element in fourth row for sm
else %no interaction model
if m>=2
mat=zeros(3,m+2); %three rows
mat(3,m+2)=1/(2*sqrt(sm2)); %only element in third row for sm
else %only one operator case
mat=zeros(2,m+2); %two rows, only for gamma and sm
mat(2,m+2)=1/(2*sqrt(sm2)); %only element in second row for sm
end;
end;
muiminusmubar=mu-mean(mu);
so=sqrt(sum(muiminusmubar.^2)/m);
if m>=2, mat(2,1:m)=muiminusmubar/(m*so); end; %for so
%can't include sp in change of variable matrix, otherwise there is a redundancy
%mat(2,m+1)=1/(2*sqrt(sp2)); %only element in second row for sp
%first row for gamma
gam=sqrt((so^2+sd2+sm2)/(sp2+so^2+sd2+sm2));
mat(1,1:m)=muiminusmubar*(1-gam^2)/gam/(sp2+so^2+sd2+sm2)/m;
so2=so^2; sp=sqrt(sp2); sd=sqrt(sd2); sm=sqrt(sm2);
%change of variables on the variance scale
mat(1,m+1) = -((so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2)) ^ (-0.1e1 / 0.2e1) * (so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2) ^ 2 / 0.2e1;
if sd2~=0,
mat(1,m+2) = ((so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2)) ^ (-0.1e1 / 0.2e1) * (0.1e1 / (sp2 + so2 + sd2 + sm2) - (so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2) ^ 2) / 0.2e1;
mat(1,m+3) = ((so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2)) ^ (-0.1e1 / 0.2e1) * (0.1e1 / (sp2 + so2 + sd2 + sm2) - (so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2) ^ 2) / 0.2e1;
else
mat(1,m+2) = ((so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2)) ^ (-0.1e1 / 0.2e1) * (0.1e1 / (sp2 + so2 + sd2 + sm2) - (so2 + sd2 + sm2) / (sp2 + so2 + sd2 + sm2) ^ 2) / 0.2e1;
end;