#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(X<x), which gives greater relative accuracy of p. This is found
#by the same method but using a contour that cuts the real axis at the saddle 
#point to the right of the origin.
#The method achieves high accuracy with a modest number N of quadrature points,
#which has been set to 100. From experimentation this appears to give close to
#machine accuracy for all x in the most challenging case of a single chi-squared
#component, but little precision is lost on setting N=50. The number L of 
#bisections used to determine the saddle points is set to 10.
if (x>1e-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
}