3
$\begingroup$

This post gives this formula to compute the empirical distribution:

\begin{equation} \hat{p}(x) = \frac{1}{m} \sum_{i=1}^m \delta(x - x^{(i)}) \tag{3.25} \end{equation}

I am trying to compute this with a concrete example:

set.seed(5)
x <- rnorm(12)
x

x stores a list of 12 variables that generated from a normal distribution.

 [1] -0.84085548  1.38435934 -1.25549186  0.07014277  1.71144087 -0.60290798
 [7] -0.47216639 -0.63537131 -0.28577363  0.13810822  1.22763034 -0.80177945

Calling ecdf

ecdf(x)

gives this output:

Empirical CDF 
Call: ecdf(x)
 x[1:12] = -1.2555, -0.84086, -0.80178,  ..., 1.3844, 1.7114

can some one gives a step by step computation for this empirical distribution with this concrete example

$\endgroup$
2
$\begingroup$

The formula you quoted is empirical generalized PDF, and the R function ecdf gives you the generalized integral of it, which is the empirical CDF (i.e. cumulative distribution function). Specifically, ecdf returns a function object that you can call, e.g. a usage may be as follows:

f = ecdf(x)
f(-1)

As the explanation notes, empirical CDF is just a step function simply calculates the following probability for any given $x$: $$\hat{F_X}(x)=\frac{\text{Number of data points that is } \leq x}{m}$$

$\endgroup$
2
$\begingroup$

Although that post didn't say so, it was using "$\delta$" in the sense of an abstract mathematical object that one can integrate. Its rules are especially simple: for any function $f,$ the integral of "$\delta(x)f(x)\mathrm{d}x$" over any set $\mathcal R$ equals $f(0)$ if $0\in \mathcal R$ and otherwise the integral is zero.

Background

This obviously is not a Riemann integral, nor is it a Lebesgue integral: it is a "generalized function" that extends the operation of integration. The concept it is generalizing is that of the integral-as-functional. That is, the rule that associates with any function $f$ its integral over a fixed region $\mathcal R$ is a "function of a function"--a functional--that enjoys some nice properties. In particular, integration is linear: the integral of a linear combination of functions is the same linear combination of their integrals. $\delta$ is one of the simplest examples of a linear functional that cannot be expressed as a (Riemann or Lebesgue) integral.

For any number $a,$ we can also "shift" $\delta$ by defining $$\int_\mathcal{R}\delta_a(x)f(x)\mathrm{d}x = f(a).$$

This permits us to abuse notation by defining

$$\delta(x-a) = \delta_a(x)$$

because it enables us (somewhat mindlessly) to apply the substitution $y=x-a$ in the integral, giving

$$\int\delta(x-a)f(x)\mathrm{d}x = \int \delta(y) f(y+a)\mathrm{d}y = f(0+a)=f(a)$$

as intended. In particular, letting $f(x)=1$ be the constant function, we have

$$\int_\mathcal{R}\delta_a(x)\mathrm{d}x = \left\{\matrix{1&\text{if }a\in \mathcal{R} \\ 0&\text{otherwise}}\right.\tag{*}$$

Application

By definition, the value of the CDF $F(t)$ of any distribution $\hat p$ is obtained by integrating that distribution up to the number $t,$ which we can perform by using the linearity property:

$$F(t) = \int_{-\infty}^t \hat{p}(x)\mathrm{d}x = \int_{-\infty}^t \sum_{i=1}^m\frac{1}{m}\delta(x-x^{(i)})\mathrm{d}x = \sum_{i=1}^m \frac{1}{m} \int_{-\infty}^t \delta_{x^{(i)}}\mathrm{d}x.$$

The integrals in the sum are immediately obtained from the basic relation $(*)$ above: they are zero except for the $i$ with $x^{(i)}\in (-\infty, t],$ where they equal $1.$ The result can be written

$$F(t) = \sum_{i\mid x^{(i)} \le t} \frac{1}{m}.$$

There is a term $1/m$ for each data value less than or equal to $t.$ Thus,

$F(t) = k/t$ where $k$ counts the number of data values less than or equal to $t.$

The graph of $F$ therefore consists of horizontal segments that jump by some multiple of $1/m$ at each data value $x^{(i)}.$ The multiple is given by the number of times $x^{(i)}$ appears in the dataset. See the figure below for an example.

Computation

An optimal algorithm for computing $F$ begins by efficiently sorting the $n$ data values. For instance, if the $n=5$ data are $1, 3, 2, 5, 3$ then in sorted order they are $1,2,3,3,5.$ The next step is to combine repeated values, summarizing them by their counts. Augment this list with starting and ending values $(-\infty,0)$ (for index $i=0$) and $(\infty,0).$ Finally, compute the cumulative sums $n_i$ of the counts and divide them by $n.$ Generally, write $(z_i, n_i/n)$ for the value $z_i$ and cumulative count $n_i$ at position $i$ in the list. For the example data, this list is

$$(-\infty, 0), (1,1/5), (2,2/5), (3,4/5), (5,5/5), (\infty, 5).\tag{**}$$

Creating this data structure takes at most $O(n \log n)$ time and $O(n)$ storage.

Given $t,$ use a binary search to find the unique $i$ for which $z_i \ge t$ and $z_{i+1}\lt t.$ Return the value of $n_i/n.$ This requires $O(\log n)$ time.

Given an ascending sequence of arguments $t,$ start with a binary search and then search forward through the list with each successive value. For very large sequences of length $N \gg n,$ this takes $O(\log(n)+N)$ time.

To plot the graph of $F,$ draw a horizontal half-open line segment from $z_i$ to $z_{i+1}$ at height $n_i/n$ for all but the last index $i.$ Each segment includes its left endpoint but not its right endpoint. This takes at most $O(n)$ time (even accounting for the clipping to the plotting window).

Example

In R, the output of ecdf is an object consisting of this list (or its logical equivalent). What is reported in the question is merely R's way of printing this object. It's not terribly useful. The help for ecdf says you can find the data for this object in its environment. For instance,

obj <- ecdf(c(1,3,2,5,3))

stores this object in obj which you can inspect via

as.list(environment(obj))

$`nobs`
[1] 5
$x
[1] 1 2 3 5
$y
[1] 0.2 0.4 0.8 1.0
$method
[1] 2
$yleft
[1] 0
$yright
[1] 1
$f
[1] 0

By comparing this to the abstract data structure $(**)$ it should now be obvious, even to those who do not use R, what most of these values represent. (f and method are used by the parent class for approximating functions in more flexible ways.)

The clearest way to inspect an ecdf object is to plot it:

plot(obj)

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.