function [res] = ProbAgreeAnalysis(n,r,cad,data)
%Function created by Nathaniel T. Stevens
%University of Waterloo
%Copyright 2015

%This function performs the P(agreement) analysis as outlined in the paper
%"Assessing Agreement between two Measurement Systems". It constructs 
%diagnostic plots, the probability of agreement plot, and outputs estimates 
%and standard errors associated with model parameters.

%n = number of parts (subjects) measured by each MS 
%r = number of replicate measurements made on each part (subject) by each MS
%cad is the upper limit of the clinically accpetable difference i.e. CAD=[-cad,cad]
%data is a column vector of data ordered first by system, then by subject,
%and then by replicate measurements

% Note that we assume MS1 is unbiased and all bias inferences are made from
% MS2 relative to MS1. This is necessary to be able to estimate all
% parameters

%Get the data
Y=data;
Y1=Y(1:n*r);
Y2=Y(n*r+1:end);

%Maximize the likelihood
grad=1; %whether or not gradient information is used during optimization

%choose the optimization option desired, comment out ones that are not wanted
options=optimset('Display','off');
%options=optimset('Display','iter');   %shows the optimization steps
if grad==1, 
    options = optimset(options,'GradObj','on');  %use gradient information
    warning on
else
    options = optimset(options,'GradObj','off');  %don't use gradient information
    warning off all   %do not display warnings that say we don't have gradient information
end;

c=[n,r,Y'];   %fixed arguments and data 

%Get some initial guesses for parameters
muinit=mean(Y);
alphainit=mean(Y2)-mean(Y1);
betainit=1;
Y1mat=reshape(Y1,r,n); %matrix whose columns represent parts, and rows are repeated measurements on those parts (MS1)
Y2mat=reshape(Y2,r,n); %matrix whose columns represent parts, and rows are repeated measurements on those parts (MS2)
s1init=sqrt(mean(var(Y1mat))); %repeatability for MS1 (pooled across individuals)
s2init=sqrt(mean(var(Y2mat))); %repeatability for MS2 (pooled across individuals)
spinit=mean([sqrt(var(Y1)-s1init^2) sqrt(var(Y2)-s2init^2)]);

f = @(x)loglik(x,c);   %create anonymous function to allow passing of fixed parameters in c    
[x,fval,exitflag] = fminunc(f,[muinit,alphainit,betainit,log(spinit),log(s1init),log(s2init)],options); %intuitive intial guesses   
%[x,fval,exitflag] = fminunc(f,[rand,rand,rand,log(rand),log(rand),log(rand)],options); %random initial guesses

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;
   [x,fval,exitflag] = fminunc(f,[muinit,alphainit,betainit,log(rand),log(rand),log(rand)],options); %pick random initial values 
end;
if tries==maxtries,  disp('MAX TRIES REACHED'),  end;

%Get MLEs given by vector x
muhat=x(1);  alphahat=x(2);  betahat=x(3);  
sphat=exp(x(4));  s1hat=exp(x(5));   s2hat=exp(x(6));   %undo log transformation

%Get Asymptotic Standard Deviations
J=FishInf(n,r,muhat,[sphat,s1hat,s2hat],[alphahat,betahat]); %fisher information matrix
Jinv=inv(J);%covariance matrix for the MLEstimators
sd=diag(sqrt(Jinv));

%Inferences for theta
if data(1)==1, %we know true values
    cad=alpha+mu*(beta-1)+1.96*sqrt((sp^2)*(beta-1)^2+s1^2+s2^2);
    k1=(-cad-alpha-mu*(beta-1))/sqrt((sp^2)*(beta-1)^2+s1^2+s2^2);
    k2=(cad-alpha-mu*(beta-1))/sqrt((sp^2)*(beta-1)^2+s1^2+s2^2);
    theta=normcdf(k2)-normcdf(k1); %true value of theta
end

%estimated versions
k1hat=(-cad-alphahat-muhat*(betahat-1))/sqrt((sphat^2)*(betahat-1)^2+s1hat^2+s2hat^2);
k2hat=(cad-alphahat-muhat*(betahat-1))/sqrt((sphat^2)*(betahat-1)^2+s1hat^2+s2hat^2);
thetahat=normcdf(k2hat)-normcdf(k1hat); %estimated theta (evaluated at MLEs)

%Delta method for SEtheta
D=DeltaMeth(cad,muhat,[sphat,s1hat,s2hat],[alphahat,betahat]); %gradient vector for theta
SEtheta=sqrt(D'*Jinv*D);

ps=muhat-2*sphat;  pm=muhat;  pl=muhat+2*sphat;

%Use the following for website software:
est=[muhat,alphahat,betahat,sphat,s1hat,s2hat,thetahat]';
sd=[sd;SEtheta];
res=[est sd]';

%check that the gradients evaluated at the MLE are zero
estcheck=[muhat,alphahat,betahat,log(sphat),log(s1hat),log(s2hat)]';
[ll,g]=loglik(estcheck,c);

ll=-ll; %log-likelihood evaluated at MLEs

%res=[ll,est']; %Used when being call by BoxCoxAnalysis

display=1;
if display==1
    ModelCheck(n,r,Y'); %construct diagnostic plots

    %Display results
    disp('        MLE      SE')
    disp(['    mu:   ',num2str(muhat),'    ',num2str(sd(1))])
    disp([' alpha:   ',num2str(alphahat),'    ',num2str(sd(2))])
    disp(['  beta:   ',num2str(betahat),'    ',num2str(sd(3))])
    disp(['sigmap: ',num2str(sphat),'    ',num2str(sd(4))])
    disp(['sigma1: ',num2str(s1hat),'    ',num2str(sd(5))])
    disp(['sigma2: ',num2str(s2hat),'    ',num2str(sd(6))])
    disp([' theta: ',num2str(thetahat),'    ',num2str(SEtheta)])
    
    %Probability of Agreement Plots
    figure
    pmin=floor(muhat-3*sphat); %lower tail of P
    pmax=ceil(muhat+3*sphat); %upper tail of P
    p=[pmin:.1:pmax]; %range of values for P
    [temp,np]=size(p);
    thetap=[];
    SEthetap=[];
    for i=1:np,
        %For thetap
        kp1=(-cad-alphahat-p(i)*(betahat-1))/sqrt(s1hat^2+s2hat^2);
        kp2=(cad-alphahat-p(i)*(betahat-1))/sqrt(s1hat^2+s2hat^2);
        thetap=[thetap;normcdf(kp2)-normcdf(kp1)];
        %Delta method for SEthetap
        Dp=DeltaMethP(cad,p(i),[s1hat,s2hat],[alphahat,betahat]);
        SEthetap=[SEthetap;sqrt(Dp'*inv(J)*Dp)];
    end
    
    %xlswrite('temp.xls',thetap)
    UCL=thetap+1.96*SEthetap; %pointwise upper confidence limits for thetap
    UCL(UCL>1) = 1;   UCL(UCL<0) = 0; %confine CI limits to [0,1]
    LCL=thetap-1.96*SEthetap; %pointwise lower confidence limits for thetap
    LCL(LCL>1) = 1;   LCL(LCL<0) = 0; %confine CI limits to [0,1]
    plot(p,thetap,'b-')
    hold on
    plot(p,UCL,'r--')
    plot(p,LCL,'r--')
    title(['Probability of Agreement Plot with c = ',num2str(cad)])
    %title(['Probability of Agreement Plot with \delta = ',num2str(cad)])
    xlabel('s')
    ylabel('\theta(s)')
    axis([pmin pmax 0 1])
    legend('P(Agreement)', '95% CI','Location','southwest')
    hold off  
%     print -r300 -dtiff PAplot.tiff
end

    
function [ll,g]=loglik(x,c)
%return the negative log likelihood function and the corresponding gradients
n=c(1);  r=c(2);
Y=c(3:end)'; %get the data from c
mu=x(1);  alpha=x(2);  beta=x(3);  
sp=exp(x(4));  s1=exp(x(5));   s2=exp(x(6));   %undo log transform
[ll,g]=likelihood_and_gradients(n,r,Y,mu,alpha,beta,sp,s1,s2);

function [ll,g]=likelihood_and_gradients(n,r,Y,mu,alpha,beta,sp,s1,s2)
%Calculate the sums of squares needed for likelihood and gradients
%Manipulate data vector to ease calculation of SSs
Y1=Y(1:n*r);   Y2=Y(n*r+1:end); %each vetor is ordered by part
mu1=mu*ones(r*n,1);   mu2=(alpha+beta*mu)*ones(r*n,1); %mean vectors for each system
Y1res=Y1-mu1;   Y2res=Y2-mu2; %residual vectors for each system
A=sum(Y1res.^2); %A
C=sum(Y2res.^2); %C
F=sum(Y1res); %F
G=sum(Y2res); %G

Y1resmat=reshape(Y1res,r,n);   Y2resmat=reshape(Y2res,r,n); %matrix columns are repeated measurements on each part
temp1=sum(Y1resmat);
temp1=temp1.^2;
B=sum(temp1); %B

temp2=sum(Y2resmat);
temp2=temp2.^2;
D=sum(temp2); %D

E=sum(Y1resmat)*sum(Y2resmat)'; %E

%log-likelihood function
ll = -r * n * log(0.2e1 * pi) - r * n * (log(s1) + log(s2)) - 0.5e0 * n * log(0.1e1 + r * sp ^ 2 * (0.1e1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) - 0.5e0 / s1 ^ 2 * A + 0.5e0 * sp ^ 2 / (0.1e1 + r * sp ^ 2 * (0.1e1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / s1 ^ 4 * B - 0.5e0 / s2 ^ 2 * C + 0.5e0 * sp ^ 2 / (0.1e1 + r * sp ^ 2 * (0.1e1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 / s2 ^ 4 * D + sp ^ 2 / (0.1e1 + r * sp ^ 2 * (0.1e1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 2 * E;

ll=-ll;  %negative so that minimization is finding the MLEs

%determine the gradients 
g=zeros(6,1);
dldmu = (1 / s1 ^ 2 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / s1 ^ 4 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 / s1 ^ 2 / s2 ^ 2) * F + (beta / s2 ^ 2 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 3 / s2 ^ 4 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 2) * G;
dldalpha = -r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 2 * F + (1 / s2 ^ 2 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 / s2 ^ 4) * G;
dldbeta = -(n * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s2 ^ 2) - (r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta * mu / s1 ^ 2 / s2 ^ 2 * F) - (r * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s2 ^ 2 / s1 ^ 4 * B) + ((mu / s2 ^ 2 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 * mu / s2 ^ 4) * G) + ((-2 * r * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s2 ^ 2 + 2 * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta) * D / s2 ^ 4) / 0.2e1 + ((-2 * r * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 / s2 ^ 2 + sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2))) * E / s1 ^ 2 / s2 ^ 2);
dldsp = -0.10e1 * n * r * sp * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) + 0.10e1 * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / (s1 ^ 4) * B - 0.10e1 * (sp ^ 3) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) / (s1 ^ 4) * B * r * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) + 0.10e1 * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * (beta ^ 2) / (s2 ^ 4) * D - 0.10e1 * (sp ^ 3) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) * (beta ^ 2) / (s2 ^ 4) * D * r * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) + (2 * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 2 * E) - (2 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s1 ^ 2 / s2 ^ 2 * E * r * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2));
dlds1 = -(r * n / s1) + 0.10e1 * n * r * (sp ^ 2) / (s1 ^ 3) / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) + 0.10e1 / (s1 ^ 3) * A + 0.10e1 * (sp ^ 4) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) / (s1 ^ 7) * B * r - 0.20e1 * (sp ^ 2) / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / (s1 ^ 5) * B + 0.10e1 * (sp ^ 4) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) * (beta ^ 2) / (s2 ^ 4) * D * r / (s1 ^ 3) + (2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s1 ^ 5 / s2 ^ 2 * E * r) - (2 * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 3 / s2 ^ 2 * E);
dlds2 = -(r * n / s2) + 0.10e1 * n * r * (sp ^ 2) * (beta ^ 2) / (s2 ^ 3) / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) + 0.10e1 * (sp ^ 4) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) / (s1 ^ 4) * B * r * (beta ^ 2) / (s2 ^ 3) + 0.10e1 / (s2 ^ 3) * C + 0.10e1 * (sp ^ 4) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) * (beta ^ 4) / (s2 ^ 7) * D * r - 0.20e1 * (sp ^ 2) / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * (beta ^ 2) / (s2 ^ 5) * D + (2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s1 ^ 2 / s2 ^ 5 * E * r) - (2 * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 3 * E);

g(1)=dldmu; %mu
g(2)=dldalpha; %alpha
g(3)=dldbeta; %beta
%need to translate from sigma scale to log(sigma) scale - apply chain rule
g(4)=dldsp*sp; %sp
g(5)=dlds1*s1; %s1
g(6)=dlds2*s2; %s2

g=-g;  %negative so that minimization is finding the MLEs

function [Jmat]=FishInf(n,r,mu,sig,bias)
%find the Fisher Information Matrix for data
sp=sig(1);   s1=sig(2); s2=sig(3);   alpha=bias(1);   beta=bias(2); 
% %expected values of the various sums of squares
A=n*r*(sp^2+s1^2);
B=n*r*(r*sp^2+s1^2);
C=n*r*((beta^2)*(sp^2)+s2^2);
D=n*r*(r*(beta^2)*(sp^2)+s2^2);
E=n*(r^2)*beta*sp^2;
F=0;
G=0;

%second partial derivative
dldmu2 = -r * n * (1 / s1 ^ 2 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / s1 ^ 4 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 / s1 ^ 2 / s2 ^ 2) - r * n * beta * (beta / s2 ^ 2 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 3 / s2 ^ 4 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 2);
dldmualpha = -r * n * (beta / s2 ^ 2 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 3 / s2 ^ 4 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 2);
dldmubeta = (2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s2 ^ 2 / s1 ^ 4 - (-2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s2 ^ 2 + 2 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta) / s1 ^ 2 / s2 ^ 2) * F + (1 / s2 ^ 2 - (-2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 4 / s2 ^ 2 + 3 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2) / s2 ^ 4 - (-2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 / s2 ^ 2 + r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2))) / s1 ^ 2 / s2 ^ 2) * G - r * n * mu * (beta / s2 ^ 2 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 3 / s2 ^ 4 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 2);
dldmusp = (-2 * r * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / s1 ^ 4 + 2 * r ^ 2 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 / s1 ^ 4 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) - 2 * r * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 / s1 ^ 2 / s2 ^ 2 + 2 * r ^ 2 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 / s1 ^ 2 / s2 ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * F + (-2 * r * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 3 / s2 ^ 4 + 2 * r ^ 2 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s2 ^ 4 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) - 2 * r * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 2 + 2 * r ^ 2 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s1 ^ 2 / s2 ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * G;
dldmus1 = (-2 / s1 ^ 3 - 2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 / s1 ^ 7 + 4 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / s1 ^ 5 - 2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 / s1 ^ 5 / s2 ^ 2 + 2 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 / s1 ^ 3 / s2 ^ 2) * F + (-2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s2 ^ 4 / s1 ^ 3 - 2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s1 ^ 5 / s2 ^ 2 + 2 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 3 / s2 ^ 2) * G;
dldmus2 = (-2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 / s1 ^ 4 * beta ^ 2 / s2 ^ 3 - 2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 4 / s1 ^ 2 / s2 ^ 5 + 2 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 / s1 ^ 2 / s2 ^ 3) * F + (-2 * beta / s2 ^ 3 - 2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 5 / s2 ^ 7 + 4 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 3 / s2 ^ 5 - 2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s1 ^ 2 / s2 ^ 5 + 2 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 3) * G;
dldalpha2 = -r * n * (1 / s2 ^ 2 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 / s2 ^ 4);
dldalphabeta = -(-2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 / s2 ^ 2 + r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2))) * F / s1 ^ 2 / s2 ^ 2 - (-2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s2 ^ 2 + 2 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta) * G / s2 ^ 4 - r * n * mu * (1 / s2 ^ 2 - r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 / s2 ^ 4);
dldalphasp = -2 * r * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 2 * F + 2 * r ^ 2 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s1 ^ 2 / s2 ^ 2 * F * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) + (-2 * r * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 / s2 ^ 4 + 2 * r ^ 2 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 / s2 ^ 4 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * G;
dldalphas1 = -2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s1 ^ 5 / s2 ^ 2 * F + 2 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 3 / s2 ^ 2 * F - 2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 / s2 ^ 4 / s1 ^ 3 * G;
dldalphas2 = -2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s1 ^ 2 / s2 ^ 5 * F + 2 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 3 * F + (-2 / s2 ^ 3 - 2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 4 / s2 ^ 7 + 4 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 / s2 ^ 5) * G;
dldbeta2 = (2 * n * r ^ 2 * sp ^ 4 * beta ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 4) - (n * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / s2 ^ 2) + (2 * r ^ 2 * sp ^ 4 * beta ^ 2 * mu * F / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 / s1 ^ 2 / s2 ^ 4) - (r * sp ^ 2 * mu * F / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / s1 ^ 2 / s2 ^ 2) + (4 * r ^ 2 * sp ^ 6 * beta ^ 2 * B / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 / s2 ^ 4 / s1 ^ 4) - (r * sp ^ 4 * B / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 / s1 ^ 4) + ((2 * r ^ 2 * sp ^ 4 * beta ^ 3 * mu / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 6 - 2 * r * sp ^ 2 * beta * mu / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / s2 ^ 4) * G) - (n * r * mu * (mu / s2 ^ 2 - r * sp ^ 2 * beta ^ 2 * mu / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / s2 ^ 4)) + ((8 * r ^ 2 * sp ^ 6 * beta ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 / s2 ^ 4 - 10 * r * sp ^ 4 * beta ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + 2 * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2))) * D / s2 ^ 4) / 0.2e1 - (r * mu * (-2 * r * sp ^ 4 * beta ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + 2 * sp ^ 2 * beta / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2))) * G / s2 ^ 4) + ((8 * r ^ 2 * sp ^ 6 * beta ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 / s2 ^ 4 - 6 * r * sp ^ 4 * beta / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2) * E / s1 ^ 2 / s2 ^ 2) - (r * mu * (-2 * r * sp ^ 4 * beta ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 / s2 ^ 2 + sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2))) * F / s1 ^ 2 / s2 ^ 2);
dldbetasp = -(2 * n * r * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s2 ^ 2) + (2 * n * r ^ 2 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s2 ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) - (2 * r * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta * mu / s1 ^ 2 / s2 ^ 2 * F) + (2 * r ^ 2 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta * mu / s1 ^ 2 / s2 ^ 2 * F * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) - (4 * r * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s2 ^ 2 / s1 ^ 4 * B) + (4 * r ^ 2 * sp ^ 5 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta / s2 ^ 2 / s1 ^ 4 * B * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) + ((-2 * r * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 * mu / s2 ^ 4 + 2 * r ^ 2 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * mu / s2 ^ 4 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * G) + ((-8 * r * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s2 ^ 2 + 8 * r ^ 2 * sp ^ 5 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 3 / s2 ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) + 4 * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta - 4 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta * r * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * D / s2 ^ 4) / 0.2e1 + ((-8 * r * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 / s2 ^ 2 + 8 * r ^ 2 * sp ^ 5 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 2 / s2 ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) + 2 * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) - 2 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * r * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * E / s1 ^ 2 / s2 ^ 2);
dldbetas1 = -(2 * n * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s2 ^ 2 / s1 ^ 3) - (2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta * mu / s1 ^ 5 / s2 ^ 2 * F) + (2 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta * mu / s1 ^ 3 / s2 ^ 2 * F) - (4 * r ^ 2 * sp ^ 6 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta / s2 ^ 2 / s1 ^ 7 * B) + (4 * r * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s2 ^ 2 / s1 ^ 5 * B) - (2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 * mu / s2 ^ 4 / s1 ^ 3 * G) + ((-8 * r ^ 2 * sp ^ 6 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 3 / s2 ^ 2 / s1 ^ 3 + 4 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta * r / s1 ^ 3) * D / s2 ^ 4) / 0.2e1 + ((-8 * r ^ 2 * sp ^ 6 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 2 / s2 ^ 2 / s1 ^ 3 + 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * r / s1 ^ 3) * E / s1 ^ 2 / s2 ^ 2) - (2 * (-2 * r * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 / s2 ^ 2 + sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2))) * E / s1 ^ 3 / s2 ^ 2);
dldbetas2 = -(2 * n * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s2 ^ 5) + (2 * n * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s2 ^ 3) - (2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 * mu / s1 ^ 2 / s2 ^ 5 * F) + (2 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta * mu / s1 ^ 2 / s2 ^ 3 * F) - (4 * r ^ 2 * sp ^ 6 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 3 / s2 ^ 5 / s1 ^ 4 * B) + (2 * r * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s2 ^ 3 / s1 ^ 4 * B) + ((-2 * mu / s2 ^ 3 - 2 * r ^ 2 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 4 * mu / s2 ^ 7 + 4 * r * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta ^ 2 * mu / s2 ^ 5) * G) + ((-8 * r ^ 2 * sp ^ 6 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 5 / s2 ^ 5 + 8 * r * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s2 ^ 3) * D / s2 ^ 4) / 0.2e1 - (2 * (-2 * r * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s2 ^ 2 + 2 * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta) * D / s2 ^ 5) + ((-8 * r ^ 2 * sp ^ 6 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 4 / s2 ^ 5 + 6 * r * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 / s2 ^ 3) * E / s1 ^ 2 / s2 ^ 2) - (2 * (-2 * r * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 2 / s2 ^ 2 + sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2))) * E / s1 ^ 2 / s2 ^ 3);
dldsp2 = -0.10e1 * n * r * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) + 0.20e1 * n * (r ^ 2) * (sp ^ 2) * ((1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) ^ 2) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) + 0.10e1 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / (s1 ^ 4) * B - 0.50e1 * (sp ^ 2) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) / (s1 ^ 4) * B * r * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) + 0.40e1 * (sp ^ 4) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3) / (s1 ^ 4) * B * (r ^ 2) * ((1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) ^ 2) + 0.10e1 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * (beta ^ 2) / (s2 ^ 4) * D - 0.50e1 * (sp ^ 2) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) * (beta ^ 2) / (s2 ^ 4) * D * r * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) + 0.40e1 * (sp ^ 4) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3) * (beta ^ 2) / (s2 ^ 4) * D * (r ^ 2) * ((1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) ^ 2) + (2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 2 * E) - (10 * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s1 ^ 2 / s2 ^ 2 * E * r * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) + (8 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta / s1 ^ 2 / s2 ^ 2 * E * r ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) ^ 2);
dldsps1 = 0.20e1 * n * r * sp / (s1 ^ 3) / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) - 0.20e1 * n * (r ^ 2) * (sp ^ 3) * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) / (s1 ^ 3) + 0.40e1 * (sp ^ 3) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) / (s1 ^ 7) * B * r - 0.40e1 * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / (s1 ^ 5) * B - 0.40e1 * (sp ^ 5) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3) / (s1 ^ 7) * B * (r ^ 2) * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) + 0.40e1 * (sp ^ 3) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) / (s1 ^ 5) * B * r * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) + 0.40e1 * (sp ^ 3) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) * (beta ^ 2) / (s2 ^ 4) * D * r / (s1 ^ 3) - 0.40e1 * (sp ^ 5) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3) * (beta ^ 2) / (s2 ^ 4) * D * (r ^ 2) * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) / (s1 ^ 3) + (8 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s1 ^ 5 / s2 ^ 2 * E * r) - (4 * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 3 / s2 ^ 2 * E) - (8 * sp ^ 5 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta / s1 ^ 5 / s2 ^ 2 * E * r ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) + (4 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s1 ^ 3 / s2 ^ 2 * E * r * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2));
dldsps2 = 0.20e1 * n * r * sp * (beta ^ 2) / (s2 ^ 3) / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) - 0.20e1 * n * (r ^ 2) * (sp ^ 3) * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) * (beta ^ 2) / (s2 ^ 3) + 0.40e1 * (sp ^ 3) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) / (s1 ^ 4) * B * r * (beta ^ 2) / (s2 ^ 3) - 0.40e1 * (sp ^ 5) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3) / (s1 ^ 4) * B * (r ^ 2) * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) * (beta ^ 2) / (s2 ^ 3) + 0.40e1 * (sp ^ 3) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) * (beta ^ 4) / (s2 ^ 7) * D * r - 0.40e1 * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * (beta ^ 2) / (s2 ^ 5) * D - 0.40e1 * (sp ^ 5) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3) * (beta ^ 4) / (s2 ^ 7) * D * (r ^ 2) * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) + 0.40e1 * (sp ^ 3) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) * (beta ^ 2) / (s2 ^ 5) * D * r * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2) + (8 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s1 ^ 2 / s2 ^ 5 * E * r) - (4 * sp / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 3 * E) - (8 * sp ^ 5 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 3 / s1 ^ 2 / s2 ^ 5 * E * r ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) + (4 * sp ^ 3 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s1 ^ 2 / s2 ^ 3 * E * r * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2));
dlds12 = (r * n / s1 ^ 2) - 0.30e1 * n * r * (sp ^ 2) / (s1 ^ 4) / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) + 0.20e1 * n * (r ^ 2) * (sp ^ 4) / (s1 ^ 6) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) - 0.30e1 / (s1 ^ 4) * A + 0.40e1 * (sp ^ 6) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3) / (s1 ^ 10) * B * (r ^ 2) - 0.110e2 * (sp ^ 4) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) / (s1 ^ 8) * B * r + 0.100e2 * (sp ^ 2) / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) / (s1 ^ 6) * B + 0.40e1 * (sp ^ 6) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3) * (beta ^ 2) / (s2 ^ 4) * D * (r ^ 2) / (s1 ^ 6) - 0.30e1 * (sp ^ 4) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) * (beta ^ 2) / (s2 ^ 4) * D * r / (s1 ^ 4) + (8 * sp ^ 6 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta / s1 ^ 8 / s2 ^ 2 * E * r ^ 2) - (14 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s1 ^ 6 / s2 ^ 2 * E * r) + (6 * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 4 / s2 ^ 2 * E);
dlds1s2 = 0.20e1 * n * (r ^ 2) * (sp ^ 4) / (s1 ^ 3) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) * (beta ^ 2) / (s2 ^ 3) + 0.40e1 * (sp ^ 6) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3) / (s1 ^ 7) * B * (r ^ 2) * (beta ^ 2) / (s2 ^ 3) - 0.40e1 * (sp ^ 4) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) / (s1 ^ 5) * B * r * (beta ^ 2) / (s2 ^ 3) + 0.40e1 * (sp ^ 6) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3) * (beta ^ 4) / (s2 ^ 7) * D * (r ^ 2) / (s1 ^ 3) - 0.40e1 * (sp ^ 4) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) * (beta ^ 2) / (s2 ^ 5) * D * r / (s1 ^ 3) + (8 * sp ^ 6 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 3 / s1 ^ 5 / s2 ^ 5 * E * r ^ 2) - (4 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta / s1 ^ 5 / s2 ^ 3 * E * r) - (4 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s1 ^ 3 / s2 ^ 5 * E * r) + (4 * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 3 / s2 ^ 3 * E);
dlds22 = (r * n / s2 ^ 2) - 0.30e1 * n * r * (sp ^ 2) * (beta ^ 2) / (s2 ^ 4) / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) + 0.20e1 * n * (r ^ 2) * (sp ^ 4) * (beta ^ 4) / (s2 ^ 6) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) + 0.40e1 * (sp ^ 6) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3) / (s1 ^ 4) * B * (r ^ 2) * (beta ^ 4) / (s2 ^ 6) - 0.30e1 * (sp ^ 4) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) / (s1 ^ 4) * B * r * (beta ^ 2) / (s2 ^ 4) - 0.30e1 / (s2 ^ 4) * C + 0.40e1 * (sp ^ 6) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3) * (beta ^ 6) / (s2 ^ 10) * D * (r ^ 2) - 0.110e2 * (sp ^ 4) / ((1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2) * (beta ^ 4) / (s2 ^ 8) * D * r + 0.100e2 * (sp ^ 2) / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * (beta ^ 2) / (s2 ^ 6) * D + (8 * sp ^ 6 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 3 * beta ^ 5 / s1 ^ 2 / s2 ^ 8 * E * r ^ 2) - (14 * sp ^ 4 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) ^ 2 * beta ^ 3 / s1 ^ 2 / s2 ^ 6 * E * r) + (6 * sp ^ 2 / (1 + r * sp ^ 2 * (1 / s1 ^ 2 + beta ^ 2 / s2 ^ 2)) * beta / s1 ^ 2 / s2 ^ 4 * E);

         
%create the matrix
Jmat=zeros(6,6);
Jmat(1,:)=[dldmu2 dldmualpha dldmubeta dldmusp dldmus1 dldmus2];
Jmat(2,:)=[dldmualpha dldalpha2 dldalphabeta dldalphasp dldalphas1 dldalphas2];
Jmat(3,:)=[dldmubeta dldalphabeta dldbeta2 dldbetasp dldbetas1 dldbetas2];
Jmat(4,:)=[dldmusp dldalphasp dldbetasp dldsp2 dldsps1 dldsps2];
Jmat(5,:)=[dldmus1 dldalphas1 dldbetas1 dldsps1 dlds12 dlds1s2];
Jmat(6,:)=[dldmus2 dldalphas2 dldbetas2 dldsps2 dlds1s2 dlds22];

%change the sign since we want the negative of the expected value instead of the expected value
Jmat=-Jmat;

function [Dvec] = DeltaMeth(c,mu,sig,bias)
%Calculate the gradient vector for theta 
sp=sig(1);   s1=sig(2); s2=sig(3);   alpha=bias(1);   beta=bias(2); 

D1 = -pi ^ (-0.1e1 / 0.2e1) * exp(-((beta * mu + alpha - c - mu) ^ 2 / (beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2)) / 0.2e1) * sqrt(0.2e1) * (beta - 1) * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.1e1 / 0.2e1)) / 0.2e1 + pi ^ (-0.1e1 / 0.2e1) * exp(-((beta * mu + alpha + c - mu) ^ 2 / (beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2)) / 0.2e1) * sqrt(0.2e1) * (beta - 1) * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.1e1 / 0.2e1)) / 0.2e1;
D2 = -pi ^ (-0.1e1 / 0.2e1) * exp(-((beta * mu + alpha - c - mu) ^ 2 / (beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2)) / 0.2e1) * sqrt(0.2e1) * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.1e1 / 0.2e1)) / 0.2e1 + pi ^ (-0.1e1 / 0.2e1) * exp(-((beta * mu + alpha + c - mu) ^ 2 / (beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2)) / 0.2e1) * sqrt(0.2e1) * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.1e1 / 0.2e1)) / 0.2e1;
D3 = -pi ^ (-0.1e1 / 0.2e1) * exp(-((beta * mu + alpha - c - mu) ^ 2 / (beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2)) / 0.2e1) * (sqrt(0.2e1) * mu * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.1e1 / 0.2e1)) / 0.2e1 - sqrt(0.2e1) * (beta * mu + alpha - c - mu) * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.3e1 / 0.2e1)) * (2 * beta * sp ^ 2 - 2 * sp ^ 2) / 0.4e1) + pi ^ (-0.1e1 / 0.2e1) * exp(-((beta * mu + alpha + c - mu) ^ 2 / (beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2)) / 0.2e1) * (sqrt(0.2e1) * mu * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.1e1 / 0.2e1)) / 0.2e1 - sqrt(0.2e1) * (beta * mu + alpha + c - mu) * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.3e1 / 0.2e1)) * (2 * beta * sp ^ 2 - 2 * sp ^ 2) / 0.4e1);
D4 = pi ^ (-0.1e1 / 0.2e1) * exp(-((beta * mu + alpha - c - mu) ^ 2 / (beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2)) / 0.2e1) * sqrt(0.2e1) * (beta * mu + alpha - c - mu) * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.3e1 / 0.2e1)) * (2 * beta ^ 2 * sp - 4 * beta * sp + 2 * sp) / 0.4e1 - pi ^ (-0.1e1 / 0.2e1) * exp(-((beta * mu + alpha + c - mu) ^ 2 / (beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2)) / 0.2e1) * sqrt(0.2e1) * (beta * mu + alpha + c - mu) * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.3e1 / 0.2e1)) * (2 * beta ^ 2 * sp - 4 * beta * sp + 2 * sp) / 0.4e1;
D5 = pi ^ (-0.1e1 / 0.2e1) * exp(-((beta * mu + alpha - c - mu) ^ 2 / (beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2)) / 0.2e1) * sqrt(0.2e1) * (beta * mu + alpha - c - mu) * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.3e1 / 0.2e1)) * s1 / 0.2e1 - pi ^ (-0.1e1 / 0.2e1) * exp(-((beta * mu + alpha + c - mu) ^ 2 / (beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2)) / 0.2e1) * sqrt(0.2e1) * (beta * mu + alpha + c - mu) * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.3e1 / 0.2e1)) * s1 / 0.2e1;
D6 = pi ^ (-0.1e1 / 0.2e1) * exp(-((beta * mu + alpha - c - mu) ^ 2 / (beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2)) / 0.2e1) * sqrt(0.2e1) * (beta * mu + alpha - c - mu) * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.3e1 / 0.2e1)) * s2 / 0.2e1 - pi ^ (-0.1e1 / 0.2e1) * exp(-((beta * mu + alpha + c - mu) ^ 2 / (beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2)) / 0.2e1) * sqrt(0.2e1) * (beta * mu + alpha + c - mu) * ((beta ^ 2 * sp ^ 2 - 2 * beta * sp ^ 2 + s1 ^ 2 + s2 ^ 2 + sp ^ 2) ^ (-0.3e1 / 0.2e1)) * s2 / 0.2e1;

Dvec=[D1;D2;D3;D4;D5;D6];

function [Dvec] = DeltaMethP(c,p,sig,bias)
%Calculate the gradient vector for thetap
s1=sig(1); s2=sig(2);   alpha=bias(1);   beta=bias(2); 

Dp1 = 0;
Dp2= -pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha - c - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 0.2e1) * sqrt(0.2e1) * (s1 ^ 2 + s2 ^ 2) ^ (-0.1e1 / 0.2e1) / 0.2e1 + pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha + c - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 0.2e1) * sqrt(0.2e1) * (s1 ^ 2 + s2 ^ 2) ^ (-0.1e1 / 0.2e1) / 0.2e1;
Dp3 = -pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha - c - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 0.2e1) * sqrt(0.2e1) * p * (s1 ^ 2 + s2 ^ 2) ^ (-0.1e1 / 0.2e1) / 0.2e1 + pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha + c - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 0.2e1) * sqrt(0.2e1) * p * (s1 ^ 2 + s2 ^ 2) ^ (-0.1e1 / 0.2e1) / 0.2e1;
Dp4 = 0;
Dp5 = pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha - c - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 0.2e1) * sqrt(0.2e1) * (p * beta + alpha - c - p) * (s1 ^ 2 + s2 ^ 2) ^ (-0.3e1 / 0.2e1) * s1 / 0.2e1 - pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha + c - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 0.2e1) * sqrt(0.2e1) * (p * beta + alpha + c - p) * (s1 ^ 2 + s2 ^ 2) ^ (-0.3e1 / 0.2e1) * s1 / 0.2e1;
Dp6 = pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha - c - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 0.2e1) * sqrt(0.2e1) * (p * beta + alpha - c - p) * (s1 ^ 2 + s2 ^ 2) ^ (-0.3e1 / 0.2e1) * s2 / 0.2e1 - pi ^ (-0.1e1 / 0.2e1) * exp(-(p * beta + alpha + c - p) ^ 2 / (s1 ^ 2 + s2 ^ 2) / 0.2e1) * sqrt(0.2e1) * (p * beta + alpha + c - p) * (s1 ^ 2 + s2 ^ 2) ^ (-0.3e1 / 0.2e1) * s2 / 0.2e1;

Dvec=[Dp1;Dp2;Dp3;Dp4;Dp5;Dp6];

function [val] = thetafun(cad,p,alpha,beta,s1,s2)
k1=(-cad-alpha-p*(beta-1))/sqrt(s1^2+s2^2);
k2=(cad-alpha-p*(beta-1))/sqrt(s1^2+s2^2);
val=normcdf(k2)-normcdf(k1);

function [res] = ModelCheck(n,r,data)

%This function constructs plots with which we can check the assumptions of
%the model to decide which formulation is best.

%n = number of parts
%r = number of repeated measurements made on each part by each MS
%data' is a (nmr x 1) vector of measured values ordered by repeated
%measurements on specific parts (MS1 in top half, MS2 in bottom half)

res=[];

close all

Y=data'; %make data a column vector
Y1=Y(1:n*r); %get MS1 measurements
Y2=Y(n*r+1:end); %get MS2 measurements
part=kron([1:1:n]',ones(r,1)); %vector of repeating 1:n, used in individuals plot

Y1mat=reshape(Y1,r,n); %matrix whose columns represent parts, and rows are repeated measurements on those parts (MS1)
Y2mat=reshape(Y2,r,n); %matrix whose columns represent parts, and rows are repeated measurements on those parts (MS2)

Y1partavg=mean(Y1mat); %average of repeated MS1 measurements on parts
Y2partavg=mean(Y2mat); %average of repeated MS2 measurements on parts

partavg=mean([Y1partavg;Y2partavg]);  %average of all repeated (MS1 & MS2) measurements on parts

[temp,n]=size(Y1partavg);

%QQ-plots of part-averages
figure
subplot(2,2,1)
oldfordqqplot(n,Y1partavg)
ylabel('Quantiles of MS1 subject averages')
title('QQ-Plot of MS1 subject averages vs Standard Normal')
% print -r300 -dtiff QQplot1.tiff

subplot(2,2,2)
oldfordqqplot(n,Y2partavg)
ylabel('Quantiles of MS2 subject averages')
title('QQ-Plot of MS2 subject averages vs Standard Normal')
% print -r300 -dtiff QQplot2.tiff

%Individuals plots of repeated residual measurements ordered by part-average ('Repeatability Plot')
plotdata1=[Y1-kron(Y1partavg',ones(r,1)),kron(Y1partavg',ones(r,1))]; %second column is used for ordering
[temp,idx1]=sort(plotdata1(:,2));
plotdata1=plotdata1(idx1,:);
subplot(2,2,3)
plot(part,plotdata1(:,1),'.b');
%plot(plotdata1(:,2),plotdata1(:,1),'.b');
xlabel('Subjects ordered by subject-average')
ylabel('Observed Residuals')
title('MS1 Residuals Measurements by Subject (ordered by subject-average)')
% print -r300 -dtiff Repplot1.tiff

plotdata2=[Y2-kron(Y2partavg',ones(r,1)),kron(Y2partavg',ones(r,1))]; %second column is used for ordering
[temp,idx2]=sort(plotdata2(:,2));
plotdata2=plotdata2(idx2,:);
subplot(2,2,4)
plot(part,plotdata2(:,1),'.b');
%plot(plotdata2(:,2),plotdata2(:,1),'.b');
xlabel('Subjects ordered by subject-average')
ylabel('Observed Residuals')
title('MS2 Residuals Measurements by Subject (ordered by subject-average)')
% print -r300 -dtiff Repplot2.tiff

function oldfordqqplot(n,X)

%Ths function creates a QQ-plot with the quantiles of simulated normal datasets (with the mean and variance of the sample data) underlayed
%It serves to better display whether the sample data is normal

%X is a dataset
%n is the length of X

mu=mean(X);
s=std(X);

h=qqplot(X); %initial QQ-plot of data and make black
set(h(1),'marker','.','markersize',15,'markeredgecolor',[0 0 0])
hold on
for i=1:50,
    Y=normrnd(mu,s,n,1); %generate normal samples based on the mean and sd of the observed sample
    temp=qqplot(Y);
    set(temp(1),'markersize',4,'markeredgecolor',[.7 .7 .7]) %change the marker size, shape (+), colour (grey) for underlayed points
    set(temp(2),'color',[1 1 1]); %make QQ-line of generated data invisible
    set(temp(3),'color',[1 1 1]); %make QQ-line of generated data invisible
end
h=qqplot(X); %put original data back on top and make black
set(h(1),'marker','.','markersize',12,'markeredgecolor',[0 0 0])
hold off