!Subroutine hCalc_MCMC_dll
!
!DLL to calculate multivariate bandwidths (h-values) for input into a product 
!kernel estimator.  Calculations use the Bayesian Markov Chain Monte Carlo 
!(MCMC) procedure of Zhang et al. (2006; X.Zhang, M. L. King, and R. J. 
!Hyndman.  A Bayesian approach to bandwidth selection for multivariate kernel 
!density estimation.  Computational Statistics & Data Analysis 50:3009-3031).  
!Arguments and other variables are as follows.
!
!Intent(In):
!  d	number of dimensions in the kernel estimation problem
!  n	number of observations (sample size)
!  u(1:d,1:n)	matrix of observed covariate values for each dimension and
!	sample
!  h0(1:d)	vector of bandwidth values used to initialize the MCMC process
!  sd(1:d)	standard deviation of Gaussian distributions used in selecting
!	new proposal values for h
!  nburnin	number of burn-in iterations to use in the MCMC process
!  niter	number of times to iterate the MCMC process after burn-in
!  kerntype(1:d)  kernel type for each dimension:
!		1 = biweight kernel (for continuous linear covariates)
!              2 = wrapped Cauchy kernel (for continuous circular covariates)
!
!Intent(Out):
!  hchain(1:d,1:niter)	record of the niter values in the MCMC chain for each 
!			dimension
!  h(1:d)   	    	vector of bandwidth values for each dimension
!  acceptance(1:d)	acceptance rates of proposal values for each 
!			dimension 
!
!Other variables:
!  rnd		uniform random number
!  numer		log of the numerator of the Metropolis-Hastings density 
!		ratio
!  denom		log of the denominator of the Metropolis-Hastings density 
!		ratio
!  ratio		Metropolis-Hastings density ratio
!  accepted(1:d)	number of proposal values accepted for each dimension
!  htry(1:d)	proposal bandwidth values
!
!Required subroutines and functions:
!  Gethtry        draws new proposal bandwidth value
!  PostCalc       calculates the value of the posterior distribution
!  Lv1OutKernCalc leave-one-out product kernel calculation
!  d_Median       returns the median value of a double precision vector
!  rnnof()        IMSL function returning a normally distributed random 
!		number (used in Gethtry)
!

subroutine hCalc_MCMC_dll(d,n,u,h0,sd,nburnin,niter,kerntype,hchain,h,acceptance)
!DEC$ attributes dllexport :: hCalc_MCMC_dll
implicit none
	
! Type declarations
integer, intent(in) :: d, n, nburnin, niter, kerntype(1:d)
real (kind=8), intent(in) :: u(1:d,1:n), h0(1:d), sd(1:d)
real (kind=8), intent(out) :: hchain(1:d,1:niter), h(1:d), acceptance(1:d)
integer :: i, j
real (kind=8) :: rnd, numer, denom, ratio
real (kind=8) :: accepted(1:d), htry(1:d)

!Initializations
h = h0
call random_seed ()	!Initialize random number generator using CPU time
call PostCalc(d, h0, n, u, kerntype, denom) !...calculate log of the 
	!posterior estimate for initial h vector; this step initializes 
	!the denominator of the density ratio 

!Metropolis-Hastings algorithm, burn-in iterations (Uses algorithm from A. 
!Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin.  2004.  Bayesian data 
!analysis.  Second edition.  Chapman & Hall/CRC Press, Boca Raton, Louisiana, 
!USA; see pp. 289 et seq.  Burn-in code is separate from chain-generation 
!code below to eliminate overhead of determining acceptance rates during the 
!burn-in period.)
do i = 1, nburnin
  htry = h	!...reintialize htry to current values of h
  do j = 1, d	!...then, taking each of d dimensions, one at a time
    call Gethtry(h(j), sd(j), kerntype(j), htry(j))	!...sample proposal 
				!distribution
    call PostCalc(d, htry, n, u, kerntype, numer)	!...calculate log of 
		!numerator of density ratio
    ratio = exp(numer - denom)	!...calculate density ratio
    if (ratio < 1.d0) then		!...apply acceptance criterion
      call random_number(rnd)
        if (rnd .le. ratio) then
          h(j) = htry(j)
          denom = numer
        end if
    elseif (ratio .ge. 1.d0) then
      h(j) = htry(j)
      denom = numer
    end if	
  end do
end do
	
!Metropolis-Hastings algorithm, final iterations (iterate niter times and 
!record values accepted at each iteration)
accepted=0  !...initialize counter of number of accepted proposal values
do i = 1, niter
  do j = 1, d	!...for each of d dimensions
    call Gethtry(h(j), sd(j), kerntype(j), htry(j))!...sample proposal 
		   !distribution
    call PostCalc(d, htry, n, u, kerntype, numer)  !...calculate log of 
		!   numerator of density ratio
    ratio = exp(numer - denom)	!...calculate density ratio
    if (ratio < 1.d0) then			!...apply acceptance criterion
      call random_number(rnd)
      if (rnd .le. ratio) then
        h(j) = htry(j)
        denom = numer
        accepted(j) = accepted(j) + 1	!...tally acceptances of htry(j)
      end if
    elseif (ratio .ge. 1.d0) then
      h(j) = htry(j)
      denom = numer
      accepted(j) = accepted(j) + 1	!...tally acceptances of htry(j)
      end if	
      hchain(j, i) = h(j)  	!...record ith value of h(j)
  end do
end do
	
!Calculate acceptance rates for proposal values
acceptance = dble(accepted)/dble(i)
	
!Set bandwidths equal to the median of values accepted over all niter 
!iterations
do i = 1, d
  call d_Median(niter, hchain(i,1:niter), h(i))
end do

return
end subroutine hCalc_MCMC_dll


!PostCalc
!
!Calculates a value proportional to the posterior log(probability) of the 
!bandwidth vector, h, given the matrix of observed sampling data, u.  The 
!calculation follows Eq. 6 of Zhang et al. (2006).  Arguments and other 
!variables are as follows.
!
!Intent(In):
!  d		number of dimensions in the kernel estimation problem
!  h(1:d)   	vector of bandwidth values for each dimension
!  n              number of observations (sample size)
!  u(1:d,1:n)	matrix of observed covariate values for each dimension and 
!		sample
!  kerntype(1:d)	kernel type for each dimension:
!	  1 = biweight kernel (for continuous linear covariates)
!             2 = wrapped Cauchy kernel (for continuous circular covariates)
!
!Intent(Out):
!  logp	log(probability) of the bandwidth vector h, conditioned on u
!
!Other variables:
!  i	counter variable
!  lambda	hyperparameter controlling shape of the prior density (Zhang et 
!	al. 2006:3014, Eq. (5))
!  logkde	product of the leave-one-out kernel density estimates, calculated 
!	as the sum of logs (Zhang et al. 2006:3013-3014, Eq. (6) and last
!	equation on p. 3013)
!  sumx	product of the prior densities of each component of h, calculated 
!	as sum of logs (Zhang et al. 2006:3014, Eq. (6))
!
subroutine PostCalc(d, h, n, u, kerntype, logp)
implicit none

! Type declarations
integer :: d, n, kerntype(1:d)
real (kind=8) :: logp
real (kind=8) :: h(1:d), u(1:d,1:n)
integer :: i
real (kind=8), parameter :: lambda = 1.d0
real (kind=8) :: logkde, sumx

!Calculate product of prior densities (as sum of logs)
sumx = 0.d0
do i = 1, d
  sumx = sumx + log(1.d0/(1.d0 + lambda * h(i)**2))
end do

!Calculate product of sampling distribution (as the sum of logs [= logkde])
call Lv1OutKernCalc(d, n, u, kerntype, h, logkde)

!Return logp
logp = sumx + logkde

return
end subroutine PostCalc


!Subroutine Gethtry
!
!Draws a proposal bandwidth value, htry, from a normal proposal distribution
!with location h and standard deviation sd, while also ensuring that htry > 0 
!when using a biweight kernel, and 0 < htry <= 1 when using a wrapped Cauchy 
!kernel.  Arguments and other variables are as follows.
!
!Intent(In):
!  h	current bandwidth value
!  sd	standard deviation of the normal proposal distribution 
!  kerntype	kernel type:
!	  1 = biweight kernel
!             2 = wrapped Cauchy kernel
!
!Intent(Out):
!  htry        	new proposal value for h
!
!Required subroutines and functions:
!  rnnof()   IMSL function returning a normally distributed random number
!
subroutine Gethtry(h, sd, kerntype, htry)
use rnnof_int
implicit none

!Type declarations
integer, intent(in) :: kerntype
real (kind=8), intent(in) :: h, sd
real (kind=8), intent(out) :: htry
	
!Draw new proposal value with appropriate constraints
select case (kerntype)
  case (1)	!...biweight kernel
100	htry = d_rnnof() * sd + h
     if (htry .le. 0.d0) goto 100  				!...ensure htry>0
  case (2)	!...wrapped Cauchy kernel
200  htry = d_rnnof() * sd + h
     if (htry .le. 0.d0 .or. htry .gt. 1.d0) goto 200	!...ensure 0<htry<=1
end select
	
return
end subroutine Gethtry


!Subroutine Lv1OutKernCalc
!
!Calculates product kernel density estimates, as per Scott (1992:150, Eq. 
!6.37), but implementing a leave-one-out strategy in order to estimate the 
!probability density at each of the observed sample points, conditioned on 
!the bandwidth vector (see Zhang et al. 2006:3013, last equation at bottom of 
!page).  Arguments and other variables are as follows.
!
!Intent(In):
!  d	number of dimensions in the kernel estimation problem
!  n	number of observations (sample size)
!  u(1:d,1:n) matrix of observed covariate values for each dimension and 
!	  sample
!  kerntype(1:d)	kernel type for each dimension:
!	  1 = biweight kernel (for continuous linear covariates)
!             2 = wrapped Cauchy kernel (for continuous circular covariates)
!  h(1:d)	vector of bandwidth for each dimension
!
!Intent(Out):
!  logkde	sum of the natural logarithms of the product kernel density 
!	estimates
!
!Other variables:
!  v(1:n)	 univariate point values, standardized to observed values, given 
!	 h(j) for that dimension
!  x	 equal to 1 – h(j) for circular covariates
!  kern(1:n) vector of length n containing univariate kernel calculations for 
!	 the current point at which density is being estimated
!  prodh	 product of bandwidths over all d dimensions
!  prodkern(1:n)	product kernel calculation at current data point for each 
!                 of the n observations
!  prodsum	 sum of prodkern over all n, minus prodkern(i), where i is the
!	 leave-one-out record
!
subroutine Lv1OutKernCalc(d, n, u, kerntype, h, logkde)
implicit none

!Type declarations
integer, intent(in) :: d, n, kerntype(1:d)
real (kind=8), intent(in) :: u(1:d,1:n), h(1:d)
real (kind=8), intent(out) :: logkde
integer :: i, j
real (kind=8) :: x
real (kind=8), parameter :: pi=3.14159265359d0
real (kind=8) :: v(1:n)
real (kind=8) :: kern(1:n), prodkern(1:n), prodsum
	
!Initializations
logkde = 0.d0

!Calculate kernel density estimates for each of n points in observed sample 
!data, leaving out the ith point
do i = 1, n
  prodkern = 1.d0	!...reinitialize product kernel estimate before estimating 
	!density at next of the n points
  !Calculate product kernel for ith point by...
  do j = 1,d
    !...calculating univariate kernels...
    v = (u(j,i) - u(j,1:n))/h(j)
    select case (kerntype(j))
      case (1) !Biweight kernel (Silverman 1986:43, Table 3.1)
        kern = 1.d0/h(j) * 15.d0/16.d0 * (1.d0 - v**2)**2 * abs(abs(v) < 1.d0)
      case (2) !Wrapped Cauchy kernel, from Batschelet (1982:282, Eq. 15.4.2)
        x = 1.d0 – h(j)
        kern = (1.d0-x**2)/(2.d0*pi*(1.d0 + x**2 - 2.d0 * x * cos(v * h(j))))
    end select
    !...and multiplying them over all d dimensions...
    prodkern = prodkern * kern
  end do
  !...then sum the resulting product kernel, subtracting the value for the 
  !ith (leave-one-out) record, and dividing by n - 1 to obtain final density 
  !estimate for the ith point, while also ensuring against a log(0) error
  prodsum = max(sum(prodkern) - prodkern(i), tiny(prodsum))
  logkde = logkde + log( 1.d0/dble(n - 1) * prodsum )
end do

return
end subroutine Lv1OutKernCalc


!Subroutine d_Median
!
!Calculates the median of a double precision vector x of length n.  Arguments 
!and other variables are as follows.
!
!Intent(In):
!  n	number of observations in x
!  x(1:n)	vector of values
!
!Intent(Out):
!  median	median of x
!
!Other variables:
!  y(1:n)	vector containing sorted values of x
!  i, j	index values
!  swap	temporary residence for values swapped during sort
!
subroutine d_Median(n, x, median)
implicit none
    
!Type declarations
integer, intent(in) :: n
real (kind=8), intent(in) :: x(1:n)
real (kind=8), intent(out) :: median
integer :: i, j
real (kind=8) :: swap
real (kind=8) :: y(1:n)
    
!Sort values
y = x
do i = 1, n-1
  do j = 2, n
    if (y(j) .lt. y(i)) then
      swap = y(i)
      y(i) = y(j)
      y(j) = swap
    end if
  end do
end do
    
!Calculate median
if (mod(n,2) .eq. 0) then   !...if n is an even number, calculate median via 
			    !linear interpolation
  i = n/2
  median = dble(y(i) + y(i+1))/2.d0
else			!...else n is odd number, so assign median directly
  i = int(dble(n)/2.d0 + 0.5d0)
  median = y(i)
end if    

return
end subroutine d_Median