How can I derive the function curve from a histogram of observed data

Question

I'm analysing some datasets that produce heavy tailed data when plotted as a histogram. My initial goal was to attempt to fit a known distribution to my dataset. Thereafter I use to the properties of a known distribution to determine proportions of values between x and say x+n.

After little success with data transformation, hurdle techniques and phase distributions I abandoned the idea of distribution fitting.

From working with KDE techniques in the past I am aware i can determine a bandwidth parameter h, using existing techniques, but I am not aware of how i could use KDE to define a curve function that would fit my data nor how i could determine the area under the curve between two distinct points on the x-axis.

I'm currently considering if it would be possible to use integral calculus to first define a curve function based on my histogram then use that result to determine the area under the curve.

My knowledge of integration is at undergrad level and sadly i never covered this this of application. From doing some initial reading it appears that Lebesgue integration may be a useful approach.

So the question that arises is, what is the most useful approach to a) determine the curve function of a dataset based on a plotted histogram and b) plot the area under the curve function.

Links to textbooks or papers would be greatly appreciated.

Answer 1

To understand what you are doing it might be useful to imagine three scenarios, in which you want to know the proportion of observations that lie in the interval $(10, 20):$

(a) You know that your data are from an exponential population with rate $\lambda = 0.1,$ thus $\mu = 10.$ Here you could use the CDF of $\mathsf{Exp}(\lambda = 0.1).$ The answer from R statistical software (or from a simple formula) is 23.25%. No data are required.

diff(pexp(c(10,20), .1))
[1] 0.2325442

(b) You know your data are from some exponential population, but do not know the rate. Then you could take $n = 100$ observations and use their sum to estimate the rate $\lambda$ as $\hat \lambda = (n - 2)/[\sum_{i=1}^{n} X_i]$ (An unbiased estimate, see the Wikipedia article on 'exponential distribution' under estimation.) and then use the estimated rate to find the desired proportion.

(c) You know nothing about the population distribution. Then you could take $n = 100$ observations and simply count the ones between 10 and 20.

In summary: (a) Gives an exact proportion directly. (b) Uses data and a knowledge of the distribution family to estimate the proportion. (c) Looks at the data and estimates the proportion by a simple count. It turns out that (c) works reasonably well, but not quite as well as (b), if feasible, because (b) is based on more information.

Here is a simulation in R with a million iterations of the 100-observation experiments mentioned above. All methods give estimates of about 23.3%:

set.seed(1115)
m = 10^6;  lam = .1; est.frac = obs.frac = lam.est=numeric(m);  n = 100
for(i in 1:m) {
  x = rexp(n, lam)
  lam.est[i] = 98/sum(x)
  est.frac[i] = diff(pexp(c(10,20),lam.est[i]))
  obs.frac[i] = mean(x > 10 & x < 20) }
diff(pexp(c(10,20), lam))
[1] 0.2325442             # exact value directly from exponential CDF
mean(est.frac);  sd(est.frac)
[1] 0.2326386             # value based on estimate of exponential rate
[1] 0.00960211            # SD of estimate
mean(obs.frac);  sd(obs.frac)
[1] 0.2325385             # proportion based on counting obs in (10,20)           
[1] 0.04217533            # SD of estimate

In any one 100-observaton experiment the 95% margin of error for the estimated proportion would be about $0.02$ for method (b) and about $0.08$ for method (c).

Note: I agree with the comments by @CloseToC and @Glen_b that it is better to use the empirical CDF than a KDE. In my R code mean(x > 10 & x < 20) is a quick version of that. However, there is relevant information in a KDE that might be useful in some applications. I illustrate below using the default KDE in R with the default bandwidth. (There are probably better choices of KDE, but this one will do for a demo.)

First, a simulated dataset, its histogram, and its default KDE:

x = rexp(100, .1);  hist(x, prob=T, col="skyblue2");
lines(density(x), lwd=2, col="red")

Here is relevant information that can be retrieved from the KDE function density. Clearly, it could be used for numerical integration, which might be worthwhile, if the fit were better.

kde.info = density(x);  kde.info

Call:
        density.default(x = x)

Data: x (100 obs.);     Bandwidth 'bw' = 2.618

        x                y            
 Min.   :-7.850   Min.   :1.716e-05  
 1st Qu.: 7.646   1st Qu.:3.146e-03  
 Median :23.143   Median :7.267e-03  
 Mean   :23.143   Mean   :1.611e-02  
 3rd Qu.:38.640   3rd Qu.:2.427e-02  
 Max.   :54.137   Max.   :6.159e-02  

Here are specific x and y-values of the KDE relevant to the interval $(10, 20).$

cond=(kde.info$x > 10 & kde.info$x < 20)
round(kde.info$x[cond],3)
 [1] 10.103 10.224 10.346 10.467 10.588 10.709 10.831 10.952 11.073 11.195
[11] 11.316 11.437 11.559 11.680 11.801 11.923 12.044 12.165 12.286 12.408
[21] 12.529 12.650 12.772 12.893 13.014 13.136 13.257 13.378 13.500 13.621
[31] 13.742 13.863 13.985 14.106 14.227 14.349 14.470 14.591 14.713 14.834
[41] 14.955 15.077 15.198 15.319 15.440 15.562 15.683 15.804 15.926 16.047
[51] 16.168 16.290 16.411 16.532 16.654 16.775 16.896 17.017 17.139 17.260
[61] 17.381 17.503 17.624 17.745 17.867 17.988 18.109 18.230 18.352 18.473
[71] 18.594 18.716 18.837 18.958 19.080 19.201 19.322 19.444 19.565 19.686
[81] 19.807 19.929
round(kde.info$y[cond],3)
 [1] 0.041 0.041 0.040 0.040 0.039 0.039 0.039 0.038 0.038 0.037 0.037 0.036
[13] 0.036 0.035 0.035 0.034 0.033 0.033 0.032 0.032 0.031 0.031 0.030 0.029
[25] 0.029 0.028 0.028 0.027 0.026 0.026 0.025 0.025 0.024 0.024 0.023 0.023
[37] 0.022 0.022 0.021 0.021 0.020 0.020 0.019 0.019 0.018 0.018 0.018 0.017
[49] 0.017 0.017 0.016 0.016 0.016 0.015 0.015 0.015 0.015 0.014 0.014 0.014
[61] 0.014 0.014 0.013 0.013 0.013 0.013 0.013 0.013 0.013 0.012 0.012 0.012
[73] 0.012 0.012 0.012 0.012 0.012 0.012 0.012 0.012 0.011 0.011

There are 82 numbers in each list, which we can interpret to get widths and heights (of sides) of 81 trapezoids. One method of numerical integration is to find the sum 0.2217 of the areas of these trapezoids, which is a relatively poor estimate of the proportion of values in $(10,20).$