
\documentclass[12pt]{article}
\begin{document}

\subsection*{Using Maple programs for ``A simple general formula...''
by Jianrong Wu and Nancy Reid}

This note describes a set of Maple programs for calculating third order
$p$-values for a continuous statistical model $f(y;\theta)$ and is
based on ``A simple general formula for tail probabilities for 
frequentist and Bayesian inference", by Fraser, Reid and Wu  to appear in 
{\it Biometrika} in June 1999. 

You are welcome to freely distribute and modify this code, but please
refer to the source on each occasion, and please send a note to 
{\tt jian@utstat.toronto.edu} if you use or change the programs.

The programs consist of seven files: mle\_template, mle\_file, pvalue\_template,
bayes\_template,
cmle\_file, pvalue\_file, bayes\_file.  The template files need to be modified
to accomodate each new problem, as explained below.   The template
files are used in the Examples section below.  The other files are given
in a separate section called Code, which is the last section of this file.

After starting Maple you may need {\tt with(stats)}.  This may be followed
by {\tt read mle\_template}, {\tt read pvalue\_template}, or 
{\tt read bayes\_template}, assuming you have already saved
the files {\tt mle\_file}, {\tt cmle\_file}, {\tt pvalue\_file} and {\tt bayes\_file}
in the same subdirectory.
The full maximum likelihood estimate can be calculated using {\tt mle\_file},
or using another program.

\subsection*{Modifying the template files}

\subsubsection*{\tt mle\_template}

The input needed in {\tt mle\_template} consists of:

\begin{itemize}
\item the number of data points, $n$
\item the number of parameter components, $p$
\item the observed data $y^0$
\item the statistical model $f=\prod f[i], i= 1..n$ (assuming independent compon
ents)
\item a starting vector of values for the parameter $\theta$, called {\tt x0}
\end{itemize}

\subsubsection*{\tt pvalue\_template}

The input needed in {\tt pvalue\_template} consists of:

\begin{itemize}
\item the number of data points, $n$ 
\item the number of parameter components, $p$
\item the observed data $y^0$
\item the statistical model $f=\prod f[i], i= 1..n$ (assuming independent components)
\item the parameter of interest $\psi$, here called {\tt g}
\item the full maximum likelihood estimate {\tt xhat}, which can be
computed using {\tt mle\_file}
\item a starting vector of values for the extended parameter $(\theta,\alpha)$,
called here {\tt x0}.
\item a pivotal quantity $z=(z_1(y_1,\theta), \cdots, z_n(y_n,\theta))$.
\item the size of the grid ($N$) and 
a grid of values for $\psi$ over which to compute the approximation,
in a vector {\tt g0[i], i=1..N}.  This grid is used twice, once to compute the
restricted mle, and in a separate loop to compute the $p$-values.  If 
another program is used to compute the restricted mle, 
the first loop can be omitted and the vector of values $\hat\theta_\psi$ 
provided in a vector called {\tt b}. The grid
defined in the statement {\tt g0:= ...} must be compatible with {\tt b}.
\end{itemize}

\subsubsection*{\tt bayes\_template}

The input needed in {\tt bayes\_template} consists of:

\begin{itemize}
\item the number of data points, $n$ 
\item the number of parameter components, $p$
\item the observed data $y^0$
\item the statistical model $f=\prod f[i], i= 1..n$ (assuming independent components)
\item the parameter of interest $\psi$, here called {\tt g}
\item the full maximum likelihood estimate {\tt xhat}, which can be
computed using {\tt mle\_file}
\item a starting vector of values for the extended parameter $(\theta,\alpha)$,
called here {\tt x0}.
\item a prior density $\pi(\theta)$, here called {\tt a}
\item the size of the grid ($N$) and 
a grid of values for $\psi$ over which to compute the approximation,
in a vector {\tt g0[i], i=1..N}.  This grid is used twice, once to compute the
restricted mle, and in a separate loop to compute the $p$-values.  If 
another program is used to compute the restricted mle, 
the first loop can be omitted and the vector of values $\hat\theta_\psi$ 
provided in a vector called {\tt b}. The grid
defined in the statement {\tt g0:= ...} must be compatible with {\tt b}.
\end{itemize}

\subsection*{Output}

The constrained maximum likelihood estimate over the grid for $\psi$
is stored in {\tt b}, a vector of length $N$ (the number of grid points).
The approximations to the tail probabilities provided are:
\begin{itemize}
\item the (first order) normal approximation to the likelihood root $r$,
returned in {\tt plr}
\item the (third order) normal approximation to using the Barndorff-Nielsen 
$r^*$ combination of $r$ and $Q$,
returned in {\tt pBN}
\item the (third order) approximation using the Lugannani and Rice combination
of $r$ and $Q$, returned in {\tt pLR}
\end{itemize}

\subsection*{Examples:}

\subsubsection*{Example 1: Maximum likelihood estimate for a log normal model}

{\tt
\begin{verbatim}
 #Example using mle_file to find maximum likelihood estimates of log
 #normal model

 with(linalg):
 #load linear algebra maple procedure to Maple session
 n:=23: # sample size
 p:=2:  # parameter dimension
 y0:=[17.88,28.92,33.00,41.52,42.12,45.60,48.48,51.84,51.96,54.12,55.58,67.80,
 68.64,68.64,68.68,84.12,93.12,98.64,105.12,105.84,127.92,128.04,173.40]: # data
 f[i]:=(1/((2*Pi)^0.5*x[2]))*exp(-(y[i]-x[1])^2/(2*x[2]^2)); # model
 f:=product(f[i],i=1..n);
 x0:=[70, 38]: # starting value
 read mle_file;
 xhat(20); # m=20 iterations of the Newton-Raphson procedure
\end{verbatim}
}

\subsubsection*{Example 2: Normal coefficient of variation (example 1 in the paper)}

{\tt
\begin{verbatim}
 with(linalg):
 n:=5: # sample size
 p:=2: # parameter dimension
 y0:=[326,302,307,299,329]: # data
 f[i]:=(1/((2*Pi)^0.5*x[2]))*exp(-(y[i]-x[1])^2/(2*x[2]^2)); # model
 f:=product(f[i],i=1..n);
 g:=x[2]/x[1]: # parameter of interest
 z:=array([seq((y[i]-x[1])/x[2],i=1..n)]); #pivotal quantity
 xhat:=[312.6,12.46755739]: # full maximum likelihood estimate
 read cmle_file;
 N:=70: # grid length 
 x0:=[ 313.0, 4.7, 0]: #starting values for $(\theta,\alpha)$

# Calculate the restricted mle $\hat\theta_\psi$ for $N$ values of $\psi$
# on a grid specified by g0

 for i from 1 to N do;
 g0:=0.0125+0.0025*i;
 b[i]:=xchat(20);
 x0:=b[i]:
 od;

# Calculate approximate $p$-value for $\psi$ for $N$ values of $\psi$
# on the same grid as used in the loop above

 read pvalue_file;
 for i from 1 to N do;
 g0:=0.0125+0.0025*i;
 xchat:=b[i];
 alphahat:=xchat[p+1]:
 pplr[i]:=plr();
 ppLR[i]:=pLR();
 ppBN[i]:=pBN();
 od;


\end{verbatim}


\subsubsection*{ Example 3.  Bayesian example: Normal coefficient of variation}

{\tt
\begin{verbatim}
 with(linalg):
 n:=5: # sample size
 p:=2: # parameter dimension
 y0:=[326,302,307,299,329]: # data
 f[i]:=(1/((2*Pi)^0.5*x[2]))*exp(-(y[i]-x[1])^2/(2*x[2]^2)); # model
 f:=product(f[i],i=1..n):
 a:=1/x[2]: # prior
 g:=x[2]/x[1]: # parameter of interest
 xhat:=[312.6,12.46755739]: # full maximum likelihood estimate
 read cmle_file;
 N:=70: # grid length
 x0:=[ 313.0, 4.7, 0]: #starting values for $(\theta,\alpha)$

# Calculate the restricted mle $\hat\theta_\psi$ for $N$ values of $\psi$
# on a grid specified by g0

for i from 1 to N do;
 g0:=0.0125+0.0025*i;
 b[i]:=xchat(20);
 x0:=b[i]:
 od;

# Calculate approximate posterior probabilities 
# for $\psi$ for $N$ values of $\psi$
# on the same grid as used in the loop above

 read bayes_file;
 for i from 1 to N do;
 g0:=0.0125+0.0025*i;
 xchat:=b[i];
 alphahat:=xchat[p+1]:
 pplr[i]:=plr();
 ppLR[i]:=pLR();
 ppBN[i]:=pBN();
 od;
\end{verbatim}
}

\pagebreak

\subsection*{Code}
{\tt 

\begin{verbatim}
###########################################################################
# This is a program to calculate mle information and should be in a file
# called  mle_file
#######################

# create parameter vector theta
x:=proc(p):
xp:= array([seq(x[i],i = 1 .. p)]):
end:

# create variable y
y:=proc(n):
yn:= array([seq(y[i],i = 1 .. n)]):
end:

# calculate the score function
score:=proc(f);
l:=simplify(log(f));
l0:=subs([seq(y[i]=y0[i],i=1..n)],l);
score:=grad(l0,x(p));
end:

# calculate the observed information matrix
j:=proc(f);
l:=simplify(log(f));
l0:=subs([seq(y[i]=y0[i],i=1..n)],l);
j:=evalm(-1*hessian(l0,x(p)));
end:

# procedure for substituting theta values to a matrix function
subxs:=proc(m,s,t);
h:=seq(col(m,i),i=1..s);
hy0:=subs([seq(x[i]=xs[i],i=1..s)],array([seq(h[i],i=1..s)]));
hy0m:=evalf(blockmatrix(1,s,seq(transpose(matrix(1,t,hy0[i])),i=1..s)));
end:

# procedure for substituting theta values to a vector function
subvxs:=proc(h,k);
subs([seq(x[i]=xs[i],i=1..k)],array([seq(h[i],i=1..k)]));
end:

# procedure for calculating maximum likelihood estimate.
xhat:=proc(m)
global xn,xs,ts,jxs,jxsinv,scorexs;
xn:=x0;
for i from 1 to m do;
xs:=xn;
scorexs:=evalf(subvxs(score(f),p));
jxs:=evalm(subxs(j(f),p,p));
jxsinv:=inverse(jxs);
ts:=multiply(jxsinv,scorexs);
xn:=evalm(xs+ts);
od;
end:
############################################################End of mle_file

    
###########################################################################
# This is a program to calculate constrained mle information and should be 
# in a file  called  cmle_file
#######################

# create parameter vector theta
x:=proc(p):
xp:= array([seq(x[i],i = 1 .. p)]):
end:

# create variable y
y:=proc(n):
yn:= array([seq(y[i],i = 1 .. n)]):
end:

# procedure for substituting theta values to a matrix function
subxs:=proc(m,s,t);
h:=seq(col(m,i),i=1..s);
hy0:=subs([seq(x[i]=xs[i],i=1..s)],array([seq(h[i],i=1..s)]));
hy0m:=evalf(blockmatrix(1,s,seq(transpose(matrix(1,t,hy0[i])),i=1..s)));
end:

# procedure for substituting theta values to a vector function
subvxs:=proc(h,k);
subs([seq(x[i]=xs[i],i=1..k)],array([seq(h[i],i=1..k)]));
end:

# calculate the score function of the tilted log-likelihood function.
scoret:=proc(f)
global x,xp;
l:=simplify(log(f));
l0:=subs([seq(y[i]=y0[i],i=1..n)],l);
l0t:=l0+alpha*(g-g0);
x := proc(p) local xp,i; xp := array([seq(x[i],i = 1 .. p)]) end;
x[p+1]:=alpha;
xp:=array([seq(x[i],i = 1 .. p+1)]);
grad(l0t,xp);
end:

# calculate the observed information matrix of the tilted log-likelihood.
jt:=proc(f)
global x,xp;
l:=simplify(log(f));
l0:=subs([seq(y[i]=y0[i],i=1..n)],l);
l0t:=l0+alpha*(g-g0);
x := proc(p) local xp,i; xp := array([seq(x[i],i = 1 .. p)]) end;
x[p+1]:=alpha;
xp:=array([seq(x[i],i = 1 .. p+1)]);
j:=evalm(-1*hessian(l0t,xp));
end:

#  procedure for calculating constrained maximum likelihood estimate.
xchat:=proc(m)
global xn,xs,scoretxs,jtxs,jtxsinv,ts,rn;
xn:=x0;
for i from 1 to m do;
xs:=xn:
scoretxs:=evalf(subvxs(scoret(f),p+1));
jtxs:=subxs(jt(f),p+1,p+1);
jtxsinv:=inverse(jtxs);
ts:=multiply(jtxsinv,scoretxs);
xn:=evalm(xs+ts);
od;
end:

###########################################################End of cmle_file

###########################################################################
# This is a program to calculate p-values and should be in a file 
# called  pvalue_file.  It computes the frequentist version of Q,
# $q_f$ in equation (3.6) and the first and third order approximations
###########################################################################

# create parameter vector theta
x:=proc(p):
xp:= array([seq(x[i],i = 1 .. p)]):
end:

# create variable y
y:=proc(n):
yn:= array([seq(y[i],i = 1 .. n)]):
end:

# calculate observed log-likelihood function
l0:=proc()
l:=log(f):
l0:=simplify(subs([seq(y[i]=y0[i],i=1..n)],l)):
end:

# calculate the gradient of log-likelihood function
l0y:=proc(y0)
l:=simplify(log(f));
ly:=map(simplify,grad(l,y(n)));
l0y:=subs([seq(y[i]=y0[i],i=1..n)],array([seq(ly[i],i=1..n)]));
end:

# procedure for substituting observed data into a matrix function
suby0:=proc(m,s,t);
h:=seq(col(m,i),i=1..s);
hy0:=subs([seq(y[i]=y0[i],i=1..n)],array([seq(h[i],i=1..s)]));
hy0m:=evalf(blockmatrix(1,s,seq(transpose(matrix(1,t,hy0[i])),i=1..s)));
end:

# procedure for substituting maximum likelihood estimate in a matrix function
subhat:=proc(m,s,t);
h:=seq(col(m,i),i=1..s);
hhat:=subs([seq(x[i]=xhat[i],i=1..p)],array([seq(h[i],i=1..s)]));
hhatm:=evalf(blockmatrix(1,s,seq(transpose(matrix(1,t,hhat[i])),i=1..s)));
end:

# procedure for substituting constrained maximum likelihood estimate 
# into a matrix function
subchat:=proc(m,s,t);
h:=seq(col(m,i),i=1..s);
hhat:=subs([seq(x[i]=xchat[i],i=1..p)],array([seq(h[i],i=1..s)]));
hhatm:=evalf(blockmatrix(1,s,seq(transpose(matrix(1,t,hhat[i])),i=1..s)));
end:

# calculate the signed square root of the log-likelihood (1.4)
R:=proc()
l0hat:=evalf(subs([seq(x[i]=xhat[i],i=1..p)],l0())):
l0chat:=evalf(subs([seq(x[i]=xchat[i],i=1..p)],l0())):
ghat:=subs([seq(x[i]=xhat[i],i=1..p)],g):
R:=signum(ghat-g0)*(2*(l0hat-l0chat))^(1/2):
end:

# calculate ancillary directions from pivotal quantity (2.1)
ancillary:=proc(z)
z0:=subs([seq(y[i]=y0[i],i=1..n)],array([seq(z[i],i=1..n)])):
z0x:=jacobian(z0,x(p)):
zy:=jacobian(z,y(n)):
z0y:=suby0(zy,n,n):
z0xhat:=subhat(z0x,p,n):
z0yhat:=subhat(z0y,n,n):
ancillary:=evalm(-1*multiply(inverse(z0yhat),z0xhat)):
end:

# derivative of the nominal reparametrization $\phi_\theta$
phix:=proc()
phi:=transpose(multiply(transpose(augment(l0y(y0))),ancillary(z))):
phii:=array([seq(phi[i,1],i=1..p)]):
phix:=jacobian(phii,x(p));
end:

# calculate the difference of the nominal reparametrization at the
# maximum likelihood estimate and at the constrained maximum likelihood 
# estimate: $\phi(\hat\theta)-\phi(\hat\theta_\psi)$
gradient:=proc(f)
phi:=transpose(multiply(transpose(augment(l0y(y0))),ancillary(z))):
phii:=array([seq(phi[i,1],i=1..p)]):
phihat:=evalf(subs([seq(x[i]=xhat[i],i=1..p)],array([seq(phii[i],i=1..p)]))):
phichat:=evalf(subs([seq(x[i]=xchat[i],i=1..p)],array([seq(phii[i],i=1..p)]))):
phid:=evalm(phihat-phichat):
end:

# the determinant of the observed information matrix $\hat j_{\theta\theta}$
jhatdet:=proc(xhat)
j:=evalm(-1*hessian(l0(),x(p))):
jhat:=subhat(j,p,p):
jhatdet:=det(jhat):
end:

# the determinant of recalibrated observed full information matrix at
# maximum likelihood estimate $\hat\j_{(\theta\theta)}$
jjhatdet:=proc()
phixhat:=subhat(phix(),p,p); 
phixdet:=det(phixhat);
jjhatdet:=jhatdet()*phixdet^(-2);
end:

# the determinant of recalibrated observed nominal information matrix at
# constrained maximum likelihood estimate $\tilde j_{(\theta\theta)}$
jlowtiltechat:=proc();
l0tilte:=l0()+alphahat*(g-g0);
jtilte:=evalm(-1*hessian(l0tilte,x(p)));
jtiltechat:=subchat(jtilte,p,p);
phixchat:=subchat(phix(),p,p);
phixchatinv:=inverse(phixchat);
jlow1:=multiply(transpose(phixchatinv),jtiltechat);
jlowtiltechat:=multiply(jlow1,phixchatinv);
end:

# calculate $\psi_{\phi}(\hat\theta_{\psi})$ (3.2)
gphichat:=proc()
phixchat:=subchat(phix(),p,p);
phixchatinv:=inverse(phixchat);
gx:=map(simplify,grad(g,x(p)));
gxchat:=subs([seq(x[i]=xchat[i],i=1..p)],array([seq(gx[i],i=1..p)]));
gphichat:=multiply(transpose(phixchatinv),gxchat);
end:

# calculate Q = special standardized maximum likelihood departure
# w1 is $\bar\psi(\hat\theta)-\bar\psi(\hat\theta_\spi)$
# w2 is $\psi_\phi'(\hat\theta_\psi) \tilde j^{(\theta\theta)}
#                                     \psi_\phi(\hat\theta_\psi)$
Q:=proc()
w1:=dotprod(gphichat(),gradient(f)); 
jupertiltechat:=inverse(jlowtiltechat());
w2:=innerprod(gphichat(),jupertiltechat,gphichat());
ghat:=subs([seq(x[i]=xhat[i],i=1..p)],g):
Q:=signum(ghat-g0)*abs(w1)*(jjhatdet()/(w2*det(jlowtiltechat())))^(0.5);
end:

# pdf of the standard normal
pdf:=proc(r):
stats[statevalf,pdf,normald](r):
end:

# cdf of the standard normal
cdf:=proc(r):
stats[statevalf,cdf,normald](r):
end:

# significance function of the signed square root of log-likelihood ratio
plr:=proc():
cdf(R());
end:

# significance function using Barndorff-Nielsen formula
pBN:=proc():
cdf(R()+log(Q()/R())/R());
end:

# significance function using Lugannani-Rice formula
pLR:=proc():
cdf(R())+pdf(R())*(1/R()-1/Q());
end:

###########################################################End of pvalue_file


###########################################################################
# This is a program to calculate marginal posterior
# p-values and should be in a file 
# called  bayes_file.  It computes the Bayesian version of Q,
# $q_B$ in equation (3.6), and the first and third order approximations
###########################################################################

# create parameter vector theta
x:=proc(p):
xp:= array([seq(x[i],i = 1 .. p)]):
end:

# create variable y
y:=proc(n):
yn:=array([seq(y[i],i = 1 .. n)]):
end:

# calculate observed log-likelihood function
l0:=proc():
l:=log(f):
l0:=simplify(subs([seq(y[i]=y0[i],i=1..n)],l)):
end:

# procedure for substituting maximum likelihood estimate into a matrix function
subhat:=proc(m,s,t);
h:=seq(col(m,i),i=1..s);
hhat:=subs([seq(x[i]=xhat[i],i=1..p)],array([seq(h[i],i=1..s)]));
hhatm:=evalf(blockmatrix(1,s,seq(transpose(matrix(1,s,hhat[i])),i=1..s)));
end:

# procedure for substituting constrained maximum likelihood estimate
# into a matrix function
subchat:=proc(m,s,t);
h:=seq(col(m,i),i=1..s);
hhat:=subs([seq(x[i]=xchat[i],i=1..p)],array([seq(h[i],i=1..s)]));
hhatm:=evalf(blockmatrix(1,s,seq(transpose(matrix(1,s,hhat[i])),i=1..s)));
end:

# calculate the signed square root of the log-likelihood (1.4)
R:=proc():
l0hat:=evalf(subs([seq(x[i]=xhat[i],i=1..p)],l0())):
l0chat:=evalf(subs([seq(x[i]=xchat[i],i=1..p)],l0())):
ghat:=subs([seq(x[i]=xhat[i],i=1..p)],g):
R:=signum(ghat-g0)*(2*(l0hat-l0chat))^(1/2):
end:

# procedure for substituting constrained maximum likelihood estimate
# into a vector function
subvchat:=proc(h,k);
subs([seq(x[i]=xchat[i],i=1..k)],array([seq(h[i],i=1..k)]));
end:

# calculate the score vector at the constrained maximum likelihood estimate
scorechat:=proc();
l:=log(f):
l0:=subs([seq(y[i]=y0[i],i=1..n)],l);
score:=grad(l0, x(p));
scorechat:=evalm(subvchat(score, p));
end;

# the determinant of information matrix at maximum likelihood estimate
jhatdet:=proc(xhat)
j:=evalm(-1*hessian(l0(),x(p))):
jhat:=subhat(j,p,p):
jhatdet:=det(jhat):
end:

# the determinant of information matrix at constrained maximum
# likelihood estimate
jtiltechat:=proc();
l0tilte:=l0()+alphahat*(g-g0);
jtilte:=evalm(-1*hessian(l0tilte,x(p)));
jtiltechat:=subchat(jtilte,p,p);
end:

# calculate the $\psi_{\theta}(\hat\theta_{\psi})$
gxchat:=proc();
gx:=map(simplify,grad(g,x(p)));
gxchat:=subs([seq(x[i]=xchat[i],i=1..p)],array([seq(gx[i],i=1..p)]));
end:

# prior evaluated at maximum likelihood estimate
ahat:=proc():
ahat:=subs([seq(x[i]=xhat[i],i=1..p)],a):
end:

# prior evaluated at constrained maximum likelihood estimate
achat:=proc():
achat:=subs([seq(x[i]=xchat[i],i=1..p)],a):
end:

# calculate Q = special standardized maximum likelihood
# departure $q_{B}$ in (3.6)
Q:=proc();
jupertiltechat:=inverse(jtiltechat());
w1:=dotprod(scorechat(),scorechat())/dotprod(gxchat(),gxchat());
w2:=innerprod(gxchat(), jupertiltechat, gxchat());
ghat:=subs([seq(x[i]=xhat[i],i=1..p)],g):
jtiltedet:=det(jtiltechat());
Q:=signum(ghat-g0)*abs((w1*w2*(jtiltedet()/jhatdet()))^(0.5))*(ahat()/achat());
end:

# pdf of the standard normal
pdf:=proc(r):
stats[statevalf,pdf,normald](r):
end:

# cdf of the standard normal
cdf:=proc(r):
stats[statevalf,cdf,normald](r):
end:

# significance function of the signed square root of log-likelihood ratio
 alphahat:=xchat[p+1]:
 pplr[i]:=plr();
 ppLR[i]:=pLR();
 ppBN[i]:=pBN();
 od;


\end{verbatim}


}
\end{document}

