#! /bin/sh
# This is a shell archive.  Remove anything before this line, then unpack
# it by saving it into a file and typing "sh file".  To overwrite existing
# files, type "sh file -c".  You can also feed this as standard input via
# unshar, or by typing "sh <file", e.g..  If this archive is complete, you
# will see the following message at the end:
#		"End of shell archive."
# Contents:  README mcv.q mcv.f
# Wrapped by scottdw@stat.rice.edu on Fri Sep 13 16:45:24 1996
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f 'README' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'README'\"
else
echo shar: Extracting \"'README'\" \(417 characters\)
sed "s/^X//" >'README' <<'END_OF_FILE'
X# 9-13-96
X
XThese routines compute 3 (multivariate) cross-validation
Xfunctions for a data vector (matrix) using the Gaussian
Xkernel.  In the multivariate case, these routines assume
Xthe same smoothing parameter along each axis.
X
XReference:  Sain, Baggerly, and Scott (1994).
X	    "Cross-Validation of Multivariate Densities."
X	    JASA, 89, 807-817. 
X
XTo install and test:
X
X% Splus
X> source("mcv.q")
X> setup()
X> q()
X%
END_OF_FILE
if test 417 -ne `wc -c <'README'`; then
    echo shar: \"'README'\" unpacked with wrong size!
fi
# end of 'README'
fi
if test -f 'mcv.q' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'mcv.q'\"
else
echo shar: Extracting \"'mcv.q'\" \(957 characters\)
sed "s/^X//" >'mcv.q' <<'END_OF_FILE'
X# 9-13-96
X# function to compute cross-validation values
X
Xsetup <- function() {
X	unix("f77 -c -O mcv.f")
X	dyn.load("mcv.o")
X	source("mcv.q")
X	print("testing")
X	x <- rnorm(100)
X	h <- seq(.1,1,by=.1)
X	motif()
X	par(mfrow=c(2,2))
X	ans <- mcv(h,x)
Xbrowser()
X	plot(ans$h,ans$ucv,type="l",main="ucv")
X	plot(ans$h,ans$bcv,type="l",main="bcv")
X	plot(ans$h,ans$boot,type="l",main="boot")
X}
X
Xmcv <- function(h,data){
X    data <- as.matrix(data)
X    d <- ncol(data)
X    n <- nrow(data)
X    nh <- length(h)
X    bcv <- double(nh)
X    ucv <- double(nh)
X    boot<- double(nh)
X    for( i in seq(along=h) ) {
X	ucv[i] <- .Fortran("ucv",as.double(data),as.integer(d),
X	    as.integer(n),as.double(h[i]),a=double(1))$a
X	bcv[i] <- .Fortran("bcv",as.double(data),as.integer(d),
X	    as.integer(n),as.double(h[i]),a=double(1))$a
X	boot[i]<- .Fortran("boot",as.double(data),as.integer(d),
X	    as.integer(n),as.double(h[i]),a=double(1))$a
X	}
X    list(h=h,ucv=ucv,bcv=bcv,boot=boot)
X}
END_OF_FILE
if test 957 -ne `wc -c <'mcv.q'`; then
    echo shar: \"'mcv.q'\" unpacked with wrong size!
fi
# end of 'mcv.q'
fi
if test -f 'mcv.f' -a "${1}" != "-c" ; then 
  echo shar: Will not clobber existing file \"'mcv.f'\"
else
echo shar: Extracting \"'mcv.f'\" \(7334 characters\)
sed "s/^X//" >'mcv.f' <<'END_OF_FILE'
XC
XC  Fortran subroutines for calculating the 
XC  three multivariate cross-validation criterion
XC  for multivariate product kernel density
XC  estimation given in Sain, Baggerly, and
XC  Scott (1994), "Cross-Validation of 
XC  Multivariate Densities," JASA, 89, 807-817. 
XC
XC  Coded and copywrite (c) by:
XC
XC  Stephan R. Sain
XC  Applied Physics Center   
XC  Pacific Northwest Laboratory
XC  P.O. Box 999 / MS K5-12
XC  Richland, WA 99352
XC  (509) 375-2383
XC  sr_sain@pnl.gov
XC
XC  The Pacific Northwest Laboratory is operated
XC  by the Battelle Memorial Institute for the
XC  Department of Energy.
XC
XC  This code was first written while the author
XC  was a graduate student in the Department of
XC  Statistics at Rice University, Houston, Tx.
XC
XC  Permission is hearby granted to freely
XC  use, copy, modify and distribute this 
XC  software for non-commercial purposes.
XC
XC  Last modified: November 2, 1994.
XC
XC  Please send bug reports or any other problems 
XC  with the this code to the above address.
XC
XC---------------------------------------------------------------------------
XC---------------------------------------------------------------------------
XC
X      subroutine ucv (data, d, n, h, b)
XC
XC---------------------------------------------------------------------------
XC  Least-Squares Cross-Validation estimation of smoothing parameters for 
XC  multivariate density estimation using a Gaussian product kernel.
XC
XC  Reference: Sain, S.R., Baggerly, K.A., and Scott, D.W. (1994),
XC  "Cross-Validation of Multivariate Densities", JASA, 89, 807-817.
XC---------------------------------------------------------------------------
XC
XC  Variables:
XC     data - a n x d matrix containing n observations of the 
XC            d-dimensional data.
XC     h    - a d-vector of smoothing parameters.
XC     b    - the value of the ucv criterion (output).
XC
XC---------------------------------------------------------------------------
XC
X      implicit double precision (a-h, o-z)
XC
X      dimension data(n,d), h(d)
X      integer d
XC
X      parameter (pi = 3.1415926535897932385)
XC
XC  Calculate constants.
X      dn = n
X      dd = d
X      pi2 = dsqrt(pi)
X      d2 = dsqrt(2.0d0)
X      hprod = 1.0d0
X      d2d = 1.0d0
X      pi2d = 1.0d0
X      do 300 k = 1, d
X         hprod = hprod*h(k)
X         h(k) = 1.0d0/h(k)
X         d2d = d2d*d2
X         pi2d = pi2d*2.0d0*pi2
X 300  continue
XC
XC  Calculate sums.
X      sum1 = 0.0d0
X      sum2 = 0.0d0
X      do 200 j = 2, n
X         do 100 i = 1, j-1
X            dij2 = 0.0d0
X            do 50 k = 1, d
X               dij = h(k)*(data(i,k) - data(j,k))
X               dij2 = dij2 + dij*dij
X 50         continue
X            if ( dij2 .gt. 240.0d0 ) goto 100
X            sum1 = sum1 + dexp(-0.25d0*dij2)
X            if ( dij2 .lt. 120.0d0 ) sum2 = sum2 + dexp(-0.5d0*dij2)
X 100     continue
X 200  continue
XC
XC  Calculate ucv estimate.
X      b = 1.0d0/(dn*hprod*pi2d)*(1.0d0+(2.0d0*sum1-4.0d0*d2d*sum2)/dn)
XC
X         return
X         end
XC
XC--------------------------------------------------------------------------
XC--------------------------------------------------------------------------
XC
X      subroutine bcv (data, d, n, h, b)
XC
XC---------------------------------------------------------------------------
XC  Biased Cross-Validation estimation of the smoothing parameters for 
XC  multivariate density estimation using a Gaussian product kernel.
XC
XC  Reference: Sain, S.R., Baggerly, K.A., and Scott, D.W. (1994),
XC  "Cross-Validation of Multivariate Densities", JASA, 89, 807-817.
XC---------------------------------------------------------------------------
XC
XC  Variables:
XC     data - a n x d matrix containing n observations of the 
XC            d-dimensional data.
XC     h    - a d-vector of smoothing parameters.
XC     b    - the value of the bcv criterion (output).
XC
XC---------------------------------------------------------------------------
XC
X      implicit double precision (a-h, o-z)
XC
X      integer d
X      dimension data(n,d), h(d)
XC
X      parameter (pi = 3.1415926535897932385)
XC
XC  Calculate constants.
X      dn = n
X      dd = d
X      pi2 = dsqrt(pi)
X      d2pi2 = dsqrt(2.0d0)*pi2
X      hprod = 1.0d0
X      dp = 1.0d0
X      pi2d = 1.0d0
X      do 100 k = 1, d
X         hprod = hprod*h(k)
X         h(k) = 1.0d0/h(k)
X         dp = dp*d2pi2
X         pi2d = pi2d*2.0d0*pi2
X 100  continue
XC 
XC  Calculate sums.
X      sum1 = 0.0d0
X      sum2 = 0.0d0
X      sum3 = 0.0d0
X      do 300 j = 2, n
X         do 200 i = 1, j-1
X            dijk = 0.0d0
X            do 150 k = 1, d
X               dij = h(k)*(data(i,k) - data(j,k))
X               dijk = dijk + dij*dij
X 150        continue
X            prod = dexp(-0.5d0*dijk) 
X            sum1 = sum1 + prod*dijk*dijk
X            sum2 = sum2 + prod*dijk
X            sum3 = sum3 + prod
X 200     continue
X 300  continue
XC 
XC  Calculate bcv estimate.
X      b = 1.0d0/(dn*hprod*pi2d)
X     &        + 1.0d0/(2.0d0*dn*dn*hprod*dp)
X     &        * (sum1 - (2.0d0*dd + 4.0d0)*sum2 
X     &        + dd*(dd + 2.0d0)*sum3)
XC
X         return
X         end
XC
XC--------------------------------------------------------------------------
XC--------------------------------------------------------------------------
XC
X      subroutine boot (data, d, n, h, b)
XC
XC---------------------------------------------------------------------------
XC  Bootstrap estimation of the smoothing parameters for multivariate
XC  density estimation using a Gaussian product kernel.
XC
XC  Reference: Sain, S.R., Baggerly, K.A., and Scott, D.W. (1994),
XC  "Cross-Validation of Multivariate Densities", JASA, 89, 807-817.
XC---------------------------------------------------------------------------
XC
XC  Variables:
XC     data - a n x d matrix containing n observations of the 
XC            d-dimensional data.
XC     h    - a d-vector of smoothing parameters.
XC     b    - the value of the boot criterion (output).
XC
XC---------------------------------------------------------------------------
X
XC
X      implicit double precision (a-h, o-z)
XC
X      integer d
X      dimension data(n,d), h(d)
XC
X      parameter (pi = 3.1415926535897932385) 
XC
XC  Calculate constants.
X      dn = n
X      dd = d
X      pi2 = dsqrt(pi)
X      d2 = dsqrt(2.0d0)
X      d3 = dsqrt(3.0d0)
X      hprod = 1.0d0
X      d2d = 1.0d0
X      d3d = 1.0d0
X      pi2d = 1.0d0
X      do 10 k = 1, d
X         hprod = hprod*h(k)
X         h(k) = 1.0d0/h(k)
X         d2d = d2d*d2
X         d3d = d3d*d3
X         pi2d = pi2d*2.0d0*pi2
X 10   continue
XC
XC  Calculate sums.
X      sum1 = 0.0d0
X      sum2 = 0.0d0
X      sum3 = 0.0d0
X      do 200 j = 2, n
X         do 100 i = 1, j-1
X            dij = 0.0d0
X            do 50 k = 1, d
X               dijk = h(k)*(data(i,k)-data(j,k))
X               dij = dij + dijk*dijk
X 50         continue
X            if ( dij .gt. 480.0d0 ) goto 100
X            sum1 = sum1 + dexp(-0.125d0*dij)
X            if ( dij .gt. 360.0d0 ) goto 100
X            sum3 = sum3 + dexp(-dij/6.0d0)
X            if ( dij .lt. 240.0d0 ) sum2 = sum2 + dexp(-0.25d0*dij)
X 100     continue
X 200  continue
XC 
XC  Calculate bootstrap estimate.
X      b = 1.0d0/(dn*hprod*pi2d)
X     &        * (1.0d0+2.0d0*(sum1/d2d+sum2-2.0d0*d2d*sum3/d3d)/dn)
XC
X         return
X         end
XC
XC--------------------------------------------------------------------------
XC--------------------------------------------------------------------------
X
X      
X
X      
X
X
X
END_OF_FILE
if test 7334 -ne `wc -c <'mcv.f'`; then
    echo shar: \"'mcv.f'\" unpacked with wrong size!
fi
# end of 'mcv.f'
fi
echo shar: End of shell archive.
exit 0