!     Last change:  DOS   3 Aug 2000    6:21 pm
!                                       *** copyright 2000 ***
! *** filename qimhof.f95               *** John F. Monahan **
!                                       **********************
      program pqimhof
!  compute the distribution of a weighted sum of chi-square variables
!   using Imhoff's method of inversion of the characteristic function
!
      implicit none
      real, dimension(10) :: wgt,nc
      real x1,x2,x3,q
      real, parameter :: delta = .01    ! use this delta for all cases
      integer, dimension(10) :: df
      integer icase,kk,m,j
      real qimhof
!
  20  format(i3,i4)
  21  format(f5.2,i3,f5.2)
  22  format(3f6.1)
  23  format(2x,'Case',i2,4x,'m=',i2,3x,' x=',f6.2,5x,'Pr(Q<x)=',f9.6)
!
!                   read cases from Davies (1980) - with number of terms
      open( unit=5, file='qimhof.dat' )
!                   put output here
      open( unit=6, file='qimhof.out' )
!
!                   several test cases 
      do icase = 1,7
!                m = number of RV's in sum
!                kk = number of points in integration
      read(5,20) m,kk
      write(6,20) m,kk
      do j = 1,m
!                wgt = weights
!                df  = degrees of freedom (integer)
!                nc  = noncentrality parameters
      read(5,21) wgt(j),df(j),nc(j)
      write(6,21) wgt(j),df(j),nc(j)
      end do  ! loop on j
!                evaluate df at x1, x2, x3
      read(5,22) x1,x2,x3
      write(6,22) x1,x2,x3
!
      q = qimhof(x1,delta,kk,m,wgt,df,nc)
      write(6,23) icase,m,x1,q
!
      q = qimhof(x2,delta,kk,m,wgt,df,nc)
      write(6,23) icase,m,x2,q
!
      q = qimhof(x3,delta,kk,m,wgt,df,nc)
      write(6,23) icase,m,x3,q
!
      end do  ! loop on i
!
      stop
      end program pqimhof
      real function qimhof(x,delta,kk,m,wgt,df,nc)
! Imhof's method for computing distribution function
!  for weighted sums of independent chi-square random variables
!
!  x     real           evaluate the distribution function here
!  delta real           spacing for integration
!  kk    integer        number of points in integration
!  m     integer        number of RV's in sum
!  wgt   real vector    weights
!  df    integer vector degrees of freedom
!  nc    real vector    noncentrality parameters
!
! See
! J P Imhof (1961) Biometrika, V 48, pp 419-426
! R B Davies (1975) Biometrika, V 60, pp 415-417
! R B Davies (1980) Applied Statistics, V 29, pp 323-333
!
      implicit none
      integer, intent(in) :: kk,m
      real, dimension(m), intent(in) :: wgt,nc
      integer, dimension(m), intent(in) :: df
      real, intent(in) :: x,delta
      real, parameter :: pio2 = 1.570796          ! pi/2
      real e,f,a,ww,u,uu,g,term
      integer k,j
!
!                   initialize sum
      qimhof = 1./2.
!                   main loop on k
!
      do k = 0,kk
!                   Davies (1980) notation
      u = real(2*k+1)*delta/2.
      uu = u * u
!                   inner loops over j      
      e = 0
      f = 0
      a = 0
      do j = 1,m
      ww = wgt(j)*wgt(j)
!                   take care of product by adding logs
      g = 1. + 4.*uu*ww
      e = e - (df(j)/4.)*log( g )
!                   other piece in exponent
      f = f  + ww * nc(j) / g
!                    contribution to angle
      a = a + df(j)*atan(2.*u*wgt(j))/2. + nc(j)*u*wgt(j)/g
!
      end do  ! loop on j
!                   angle
      a = a - u*x
!                   exponent
      term = 0.
      e = e - 2.*uu*f
      if( e .gt. -30. ) term = exp(e)
!                   finish term
      term = term * sin(a) / ( pio2*real(2*k+1) )
!                   add to sum
      qimhof = qimhof - term
!                   done
      end do  ! loop on k
      return
      end function qimhof
! *** end of filename qimhof.f95        *********************