function [res]=baselineMLEresults(k,m,n,b,model,SPdata,BASEdata)
%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
%BASEdata is a column vector also ordered first by parts and then by operators

%Sample Call:
%[res]=baselineMLEresults(6,3,2,24,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]);

typeflag=1; %ensure we're dealing with baseline (type A augmentation) case
kab=b;
APdata=BASEdata;
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;