Friday, 17 January 2014

Particle approximation to probability density functions: Dirac delta function representation

In the previous post, I have briefly shown the idea of using dirac delta function for discrete data representation. In the second example there, a histogram locations for a given set of points are presented as spike trains, where as heights are somehow given in a second sum. This is hard to follow and visualise, of course if you are not that good in reading formulation with different indexes. Due to pedagocial reasons, an easier representation of arbitrary probability density function (PDF), $p(x)$, one would simply need to couple each discrete points with a corresponding weight.

Hence, a set  $\{x_{i}, \omega^{i}\}_{i=1}^{N}$ would be an estimation of PDF, $\hat{p}(x)$ . At this point we can invoke dirac delta function,

$ \hat{p}(x) = \sum_{i=1}^{N} \omega^{i} \delta(x-x_{i})$

Let's revisit the R code given there, this time let's draw uniform numbers between $[-2, 2]$ to get 100  $x_{i}$ values. Simply these numbers will indicate the locations on the x-axis, a spike train.  For simplicity, let's use Gausian distribution for target PDF, $ \mathcal{N}(0, 1)$.  Than, for weights we need to draw numbers using the spike locations. This approach is easier to understand compare to my previous double index notation.

R Example code

Above explained procedure is trivial to implement in R.

# Generate locations 100 x locations
# out out 1000 points in [-2.0, 2.0]
# Domain where Dirac comb operates
Xj = seq(-2,2,0.002) 
Xi = sample(Xj, 100)
# Now generate weights from N(0,1) at those given locations
Wi = dnorm(Xi)
# Now visualise
plot(Xi, Wi, type="h",xlim=c(-2.0,2.0),ylim=c(0,0.6),lwd=2,col="blue",ylab="p")


Above notation introduces second  abuse of notation while actually there must be a secondary regular grid that pics $x_{i}$ values using dirac delta in practice. Because, the argument of $\hat{p}(x)$, is in the discrete domain. So a little better notation that reflects the above code would be

$ \hat{p}(x_j) = \sum_{i=1}^{N} \omega^{i} \delta(x_{j}-x_{i})$

The set $x_j$ is simply defined in a certain domain, for example regularly. Hence I only recommend not to introduce dirac delta for explaining a particle approximation to PDFs for novice students in the class. It will only confuse them even more.
Figure: Spike trains with weights $\hat{p}(x) = \sum_{i=1}^{N} \omega^{i} \delta(x-x_{i})$

(c) Copyright 2008-2015 Mehmet Suzen (suzen at acm dot org)

Creative Commons Licence
This work is licensed under a Creative Commons Attribution 3.0 Unported License.