Hellinger Distance between 2 vectors of data points using cumsum in R I have numerous vectors of data points and I want to compute the hellinger distance between the probability distributions of every 2 vectors. I am using this version of the Hellinger distance equation: $\frac{1}{\sqrt{2}} \sqrt{\int\ (\sqrt{f(x)} - \sqrt{g(x)})^2}$ if f and g are density functions. 
This is the code I used to calculate the hellinger distance using the density() function in R and cumsum() to find the integral. Essentially, instead of finding a function to approximate the density, I used the actual y values given by the density() function and made sure the densities are calculated over the same x values. However, I am not sure if this is a correct approximation of the hellinger distance. How can I calculate hellinger distance without approximating functions and just using the density values?
There are R packages such as textmineR calculating the hellinger distance but it requires that the vector of data points being compared have the same length which isn't the case for me. Another package statip, uses the density() function in R and then approximates a linear function to the density when computing the integral but I think using a linear function of a density is problematic.
For example, I compared the hellinger distance from my function and the hellinger distance from the CalcHellingerDist function in textmineR for x,y vectors of equal length and get very different answers (.795 vs .167). Now, I do not know if they are using the same hellinger equation as it's unclear from their github/cran (https://cran.r-project.org/web/packages/textmineR/textmineR.pdf) but regardless I wanted to verify if my function makes sense given how different the values are.
This is my code. 
hellingerD <- function(x,y){ #x and y are numeric vectors
  d1 = density(x)
  d2=density(y)
  #calculating densities for the same x values
  den1 = density(x,to=max(d1$x,d2$x),from=min(d1$x,d2$x),n=16384) 
  den2 = density(y,to=max(d1$x,d2$x),from=min(d1$x,d2$x),n=16384)
  #using 16384 points to have smaller spaced points when calculating integral
  z=((den1$y^0.5) - (den2$y^0.5))^2
  helldist=sqrt(((cumsum(z[1:length(z)-1]) * (den1$x[2]-den1$x[1]))[16383])/2)
  return(helldist)
} 
library(textmineR)
x <- rchisq(n = 100, df = 8)
y <- x^2
hellingerD(x,y)
[1] 0.7954878
CalcHellingerDist(x = x, y = y)
[1] 0.1669515

 A: There are many issues here: some statistical, some numerical.
Statistical issues
The chief statistical issue is that the Hellinger distance between two samples of random distributions is not defined.  We have to decide whether the purpose is (a) to estimate the Hellinger distance of the underlying distributions or (b) to produce some Hellinger-distance like measure of discrepancy between the two empirical distributions.  The problem with (a) is that just about anything we might try is going to be biased (high) due to the random deviations in the samples.  Dealing with that will take us onto a different track, so instead I will address (b).
Numerical issues
The code in the problem replaces the empirical data with kernel density estimates (KDE) and performs a numerical integration of the KDEs which, because they represent continuous distributions, do have a well-defined Hellinger distance.  The numerical issues this code encounters are

*

*The KDE replaces each data value with a multiple of a continuous distribution centered there (the "kernel").  Thus, the KDE extends beyond the original range of the data.


*This linear combination of densities is evaluated on a discrete set of bins, thereby creating some discretization error.


*The KDE depends strongly on the choice of kernel bandwidth.  Default choices in software are based on different objectives and might not be suitable for the present purpose.


*In order to compare two KDEs quantitatively, they must (a) cover a common range of values and (b) discretize that range identically.  In effect, each KDE is represented as a vector of estimated density values in a set of bins.
Analysis
In mathematical notation, given two data sets $x=(x_1,\ldots, x_m)$ and $y=(y_1, \ldots, y_n),$ we need to find (i) a suitable kernel bandwidth $h$ and (ii) an equally-spaced set of points $z_1, z_2=z_1+h, \ldots, z_N = z_1 + (N-1)h$ that covers not only the full range of values in $x$ and $y$ but also extends beyond that range by several multiples of $h.$  Each KDE is a parallel set of density estimates $f_x(i)$ and $f_y(i),$ $i=1,2,\ldots, N.$  Because these are density estimates, we therefore expect their integral to equal unity.  The Riemann sum approximation to that integral is appropriate, given the discrete nature of these data structures, whence we expect that
$$h \sum_{i=1}^N f_x(i) \approx 1.\tag{*}$$
For whatever reason, the density function in R creates densities that sum to $1+1/(2N),$ give or take some floating point error in the calculation.  Although this is a small error, we can compensate by dividing the values returned by density by the value in $(*)$ before using them.
Assuming this normalization has been done (and still calling the normalized densities $f_x$ and $f_y$), the squared Hellinger Distance can be estimated with the Riemann sum
$$\operatorname{HD}(x,y)^2 = 1 - h \sum_{i=1}^n \sqrt{f_x(i)f_y(i)}.$$
Assuming the range of the bins $z_i$ has been extended sufficiently far, so that $f_x(1)\approx f_y(1)\approx f_x(N)\approx f_y(N)\approx 0,$ this is essentially the Trapezoidal Rule calculation (which is equivalent to integrating a linear spline of the data--thus, there should be no concern whatsoever about using such linear approximations).  It therefore should be accurate to $O(h^2).$
Testing issues
The datasets created in the problem are challenging because there is little overlap in their ranges and the extent of $y=x^2$ is hugely greater than the extent of $x.$  This makes it difficult to represent both datasets accurately with the bins $z_i.$

In this histogram of the combined data, the $y$ values are distinguished by the red ticks in the rug plot at the bottom; the $x$ values are shown with black ticks (all at the far left).
Statistical issues redux
Now that we have reviewed some of the numerical issues, let's return to the basic statistical one: what to do about the dependence of the KDE--and thence the Hellinger Distance--on the bandwidth?  I have no answer, but do have a suggestion: study how the distance depends on the bandwidth so you can determine how sensitive the result is.
Such a study is relatively simple to carry out: control the bandwidth in the calculation through an argument to the function.  Systematically vary that argument around the default bandwidth.  Plot the resulting distances against the amount of variation.  If the plot is reasonably unvarying (horizontal), then bandwidth doesn't matter.  But if the plot varies enough to matter, you will have to study this dependency more closely--and you won't want to rely on a single default result.
Results
The function f in the R code below handles all the issues raised here.  As a test, I generated a dataset of 150 $x$ values and 50 $y$ values from the Standard Normal distribution and created four datasets of values $y+\mu,$ $\mu=0,1,2,3,$ to compare to those $x$ values by means of the Hellinger distance plot described in the preceding section.  As a reference, I also plotted the Hellinger distances between the underlying Normal$(\mu,1)$ distributions (which appear as dotted horizontal lines).

"Multiplier" is the multiple of a default bandwidth selected by density.  We should expect distances to increase with smaller multiples, because they distinguish individual data points, thereby enlarging the distances.  (In the limit of a zero multiple the distances will be $1$ unless some of the $x$ and $y$ values are identical.)  For multiples larger than $1,$ the KDE is making both datasets look more and more Normal, so in the limit the distance will plateau at a value of $\sqrt{1-\exp(-d^2/8)}$ where $d-\bar y - \bar x$ is the difference in sample means.
Evidently the empirical distance is fairly constant and close to the underlying Hellinger Distance for larger $\mu.$  For smaller $\mu,$ where $x$ and $y$ are not easy to distinguish, the empirical calculation depends noticeably on the bandwidth and appreciably overestimates the underlying Hellinger Distance.  I expect these qualitative observations to hold generally.
Answer to the question
Let's study the data in the question.  Here is the code:
# Generate the datasets
set.seed(17)
x <- rchisq(n = 100, df = 8)
y <- x^2

# Compute points in the HD plot
extent <- exp(seq(log(1/10), log(10), length.out=51))
hd <- f(x, y, extent=extent, cut=20, n=2^14)

# Graph those points
plot(extent, hd, log="x", type="l", lwd=2,
     xlab="Multiplier", ylab="Hellinger Distance")
abline(v=1, lwd=2, col="Gray", lty=3)
abline(h=0.7955, col="Red", lwd=2)


The plot is in black.  The red line shows the value of $0.7955$ reported by HellingerD.  It's near the high end of this plot, corresponding to a bandwidth only $0.2$ times as great as the default bandwidth.  Over this range of multiples, the computed Hellinger distance ranges from over $0.8$ down to $0.5.$  If this amount of variation is important in your application, then you better watch out: the answer is sensitive to the choice of bandwidth.  The value computed for the default R bandwidth (shown at a multiple of $1$) is
f(x,y)


[1] 0.7505968



#
# Compare two empirical distributions.
#
f <- function(x, y, n=2^9, extent=1, tol=1e-2, cut=3, ...) {
  # Estimate bandwidths for KDEs
  z <- c(x,y)
  d <- density(z, n=n, cut=cut, ...) # Compute a default bandwidth d$bw
  r <- range(z) + cut*c(-1,1)        # Range of future calculations
  
  # Compute Hellinger distances for various bandwidths.
  HD <- function(x, y, dx) {
    # Normalize the densities
    c.x <- sum(x)*dx * 2*length(x) / (2*length(x) + 1)
    c.y <- sum(y)*dx * 2*length(y) / (2*length(y) + 1)
    if (abs(c.x - 1) > tol || abs(c.y - 1) > tol) 
      warning("Normalization factor errors are ", 1/(c.x-1), " and ", 1/(c.y-1))
    x <- x / c.x
    y <- y / c.y
    
    # Return the HD between the KDEs
    sqrt(max(0.0, 1 - sum(sqrt(x*y)) * dx))
  }
  
  # Apply `HD` to `extent` times the estimated bandwidth in `d$bw`.
  sapply(extent*d$bw, function(bw) {
    d.x <- density(x, from=r[1], to=r[2], n=n, width=bw, ...)
    d.y <- density(y, from=r[1], to=r[2], n=n, width=bw, ...)
    dx <- diff(d.x$x[1:2])
    HD(d.x$y, d.y$y, dx)
  })
}
#
# Test the comparison.
#
set.seed(17)
n <- 50
x <- rnorm(3*n)
y <- rnorm(n)

# Create a plotting region, axis labels, etc.
extent <- seq(log(1/3), log(3), length.out=11)
plot(range(exp(extent)), c(0,1), type="n", log="x", lwd=2, 
     xlab="Multiplier",
     ylab="Hellinger Distance")
abline(v=1, lty=3, lwd=2, col="Gray")

# Study how HD varies with bandwidth multiplier for four shifted versions of `y`.
invisible(lapply(0:3, function(mu) {
  n0 <- 3*(length(x)+length(y))
  h <- f(x, y+mu, extent=exp(extent), cut=3, n=n0, tol=0.5/n0)
  lines(exp(extent), h, col=hsv(mu/5,.8,.8), lwd=2)
  abline(h = sqrt(1 - exp(-mu^2/8)), lty=3, col=hsv(mu/5,.8,.6), lwd=2)
}))

