#source this BEFORE sourcing multisteptestfun.r
Pwtdchisq=function(x,lam)
{
#Returns the value of P=P(X>x) where X is a linear combination, with positive
#weights held in the vector lam, of independent chi-squared variates each on
#one degree of freedom.
#Uses the method of Talbot (J. Inst. Maths. Applics.1979,vol23, pp 97-120)
#for inversion of a Laplace Transform. Let the mgf of X, with negative argument
#to be consistent with Laplace Transform methods, be g(s). The required P is
#given by the integral of g(s)exp(sx)/s along a line parallel to the imaginary
#axis and located to the left of the origin but to the right of the highest
#negative singularity. Talbot's method replaces this line with a contour cutting
#the real axis at the same location, symmetric about the real axis, but with
#asymptotes at +/- sigma -i infinity for a scale factor sigma.
#The point of intersection of the contour with the real axis is determined close
#to a saddle point by applying bisection to find the approximate zero of the
#derivative of the integrand on the negative real axis. The second derivative
#at this point is used to determine the scale factor sigma of the contour.
#In fact this contour is only used if X > sum(lam). For X<=sum(lam) we evaluate
#P=1-p where p=P(X1e-30*sum(lam)){
lam2=lam*2
nlam=length(lam)
L=10; # loop number for bisection
twopi=2*pi
case=(x>sum(lam))
if (case)
{spt=-(1/max(lam2))/2 # semi-bound on negative saddle point
step=-spt}
else
{spt=((1+nlam/2)/x)/2 # semi-bound on positive saddle point
step=spt}
for (k in (1:L)){
fval=(sum(0.5*lam2/(1+lam2*spt)) +1/spt - x)
step=step/2
spt=spt-( 2*(fval<0)-1)*step
}
sigma=1/sqrt(0.5*sum((lam2/(1+lam2*spt))^2)+1/spt^2)
#P1 below is a Laplace type saddle point appproximation to P which is not
#returned but may be printed for comparison.
P1=exp(-sum(log(1+lam2*spt))/2) *(1/ abs(spt)) * exp(spt*x) * sigma/sqrt(twopi)
tau=4*sigma # tau and sigma used for transforming contour
shift=spt-tau
N=50 # number of quadrature intervals
grid=seq(0,N)
N1=N+1
grid=grid/N
rho=rep(0,N1)
theta=dg=zgrid=sgridrl=sgridim=dsgridrl=dsgridim=
tgridrl=tgridim=rad=ftrl=ftim=rho
zgrid=twopi*grid #grid values of z
#determine contour grid points s as function of grid points of z
sub = seq(2,N)
sgridrl[sub]=zgrid[sub]/(2*tan(zgrid[sub]/2)) # radius of polar form of s-grid
sgridrl[c(1,N1)]=c(1,1) # set points at zero and infinity to 1
sgridim=zgrid/2
# form derivative of grid map, imaginary then real part
dsgridim[sub]=
-(sin(zgrid[sub]/2)*cos(zgrid[sub]/2)-zgrid[sub]/2)/(2*sin(zgrid[sub]/2)^2)
dsgridim[c(1,N1)]=c(0,1)
dsgridrl=rep(1/2,N1)
# linear transformation for contour
tgridrl=tau*sgridrl+shift
tgridim=tau*sgridim
# form squared reciprocal product of terms in G
rho=tgridrl^2+tgridim^2
theta=2*acos(tgridrl/sqrt(rho))
for (dlam in lam2) {
# accumulate products of function terms
rad=sqrt((1+dlam*tgridrl)^2 + (dlam*tgridim)^2)
arg=acos((1+dlam*tgridrl)/rad)
rho=rho*rad
theta=theta+arg
}
# incorporate other factors
rho=tau*exp(tgridrl*x)/sqrt(rho)
theta=-theta/2+tgridim*x
# incorporate derivative of S
ftrl=rho*cos(theta)
ftim=rho*sin(theta)
gtrl=-(ftrl*dsgridrl-ftim*dsgridim)
#plot(gtrl) #should be approximatedly 'normal' in shape
gtrl[1]=gtrl[1]/2
gtrl[N1]=0
P2=sum(gtrl)*2/N
if (1-case) {
P1=1-P1
P2=P2+1
}
} # end of comptations for x>1.e-30*sum(lam)
else {P2=1}
P2
# end of function Pwtdchisq
}