Mean conservation for density estimation via diffusion using the finite element method
Keith Y. Patarroyo


Seeing What Is Not There:
 Learning Context to Determine Where Objects Are Missing
    Jin Sun     David W. Jacobs

Department of Computer Science

University of Maryland

{jinsun,djacobs}@cs.umd.edu

    December 30, 2023
=====================================================================================================================



We propose  boundary conditions for the diffusion equation that maintain the initial mean and the total mass of a discrete data sample in the density estimation process. A complete study of this framework with numerical experiments using the finite element method is presented for the one dimensional diffusion equation, some possible applications of this results  are presented as well. We also comment on a similar methodology for the two-dimensional diffusion equation for future applications in two-dimensional domains.






§ INTRODUCTION


Estimating a density function using a set of initial data points in order to find probability information is a very significant tool in statistics<cit.>. The method of Kernel Density Estimation (KDE)<cit.> is now standard in many analysis and applications. Furthermore, this idea  has been applied in multiple fields (Archaeology <cit.>, Economy <cit.>, etc). The author of this article is particularly interested in constructing Perception of Security (PoS) hotspots using (KDE) methods to analyze real data registered by security experts in Bogotá <cit.>. 

Nowadays a wide variety of methods are available to find density functions (KDE) <cit.>,<cit.>. The  method of KDE via difussion is of particular interest for this document; a recent article  <cit.> develops a systematic method for (KDE) using the diffusion equation, also they propose a more general equation to solve some biases for data estimation. However in their analysis, it is only considered the normalization (conservation of mass) of the density function via Neumann boundary conditions, the mean of the sample data is not considered, thus inducing a change of an important initial parameter  from the discrete data sample.
 
In this article, we propose a new set of boundary conditions for the diffusion equation that maintain the initial mean and mass of the the discrete data sample in the density estimation process. A complete study of this framework is performed  using the finite element method (FEM) to solve the one-dimensional diffusion equation for different boundary conditions. We show the induced error on the final density when the mean is not conserved. We also show how this one-dimensional model can be used to simulate a (PoS) in a busy avenue of a city. Lastly the new boundary conditions are presented for the two-dimensional diffusion equation for future applications in two dimensional domains. 



§ DIFFUSION EQUATION WITH DIFFERENT BOUNDARY CONDITIONS


As it was first noted in <cit.> and expanded in <cit.>, solving the diffusion equation with a discrete data sample {b_n}_n=1^N as initial condition (<ref>)  give an estimate of a continuous probability density function. Then by solving the diffusion equation <cit.>,


  ∂u(x,t)/∂t- ∂^2 u(x,t)/∂x^2=0     a<x<b , t>0, 
   

     u(x,0)=1/N∑_i=1^Nδ(x-b_i),   x,b_i  ∈[a,b] , 


with appropriate boundary conditions and then finding the best t (bandwidth) for the initial data sample one obtains a continuous estimation of the experimental density. In this article we do not consider algorithms for bandwidth selection, we consider only the conservation of the mean. For more information on the bandwidth selection see <cit.>.

This one-dimensional toy problem is nevertheless of interest in applications for constructing (PoS). For instance we can model an avenue as a one dimensional domain where predictions of the most dangerous places in a selected  zone can be  accomplished.

In the following sections we present the  non-conservation of the mean for the Neumann boundary conditions for  Problem (<ref>). We also propose new boundary conditions. For the derivations we assume that the functions are sufficiently smooth in order for the theorems of vector analysis to hold. Moreover the following derivations can be done for a more general diffusion equation with a variable diffusion coefficient k(x).



 §.§ Neumann boundary conditions


If we consider the Neumann or natural boundary conditions on the Problem (<ref>), we have


    ∂ u(x,t)/∂ x|_a=0   ,   ∂ u(x,t)/∂ x|_b=0 .


As is widely known, the total mass is conserved over time, see  Section <ref>, however the mean of the initial condition  is, in general, not conserved. Indeed, we have


    d/dt( ∫_a^b x u(x,t) dx )     =  ∫_a^b x  ∂^2 u(x,t)/∂ x^2 dx
        = [ x ∂ u(x,t)/∂ x]_a^b- [u(x,t)]_a^b
        = u(a,t) - u(b,t).


Where we used (<ref>), (<ref>) and integration by parts. Hence the mean is generally not conserved, it depends on the values of u(x,t) at the boundary in a time t.



 §.§ Boundary conditions that conserve the mean


We propose the following boundary conditions for (<ref>),


    ∂ u(x,t)/∂ x|_a=∂ u(x,t)/∂ x|_b  ,   u(b)-u(a)/b-a = ∂ u(x,t)/∂ x|_b .


Note that this boundary conditions are non-local, we need to evaluate in  both boundary points at the same time. Now we show that both the mean and the mass are conserved over time using this boundary conditions. Consider first the conservation of the total mass. We have,


    d/dt( ∫_a^b u(x,t) dx )  =  ∫_a^b∂^2 u(x,t)/∂ x^2 dx
     = [  ∂ u(x,t)/∂ x]_a^b = ∂ u(x,t)/∂ x|_a-∂ u(x,t)/∂ x|_b=0.


Where we used (<ref>), (<ref>) and integration by parts. This shows that the total mass is conserved. Consider now the conservation of the mean. We have,


    d/dt( ∫_a^b x u(x,t) dx )     =  ∫_a^b x  ∂^2 u(x,t)/∂ x^2 dx
        = [ x ∂ u(x,t)/∂ x]_a^b- [u(x,t)]_a^b
        = (b-a)∂ u(x,t)/∂ x|_b -u(b,t) + u(a,t)
        = 0.


Again (<ref>), (<ref>) and integration by parts were used to obtain the desired result.

This shows that the boundary conditions (<ref>) for  problem (<ref>) conserve both mean and mass. Now we proceed to make some numerical simulations using FEM to show the consequences of the application of this boundary conditions in the process of estimation a probability density for a data sample (<ref>). 



§ NUMERICAL STUDY OF MEAN CONSERVATION 


Now the problem (<ref>),(<ref>) is written in a weak formulation <cit.>  in order to apply the finite element method to the  problem. Now for all v(x) ∈ C^∞(a,b) we have,


    ∫_a^b∂ u(x,t)/∂ t v(x)dx + ∫_a^b∂ u(x,t)/∂ xd v(x)/d xdx =(v(b)-v(a))∂ u(x,t)/∂ x|_b.


We solve this weak formulation using FEM with low order elements in theinterval [a,b]=[0,10], where the number of elements is M. Then  Problem (<ref>),(<ref>),(<ref>) yields the problem in the discretised space V^h. Find u(x,t) ∈ V^h, such thatfor all v(x) ∈ V^h: 




   ∫_a^b ∂u(x,t)/∂t v(x)dx + ∫_a^b ∂u(x,t)/∂x d v(x)/d xdx =(v(b)-v(a))∂u(x,t)/∂x|_b, 
   

    u(x,0)=M/(b-a)N∑_i=1^Nδ(x-b_i),   x,b_i  ∈[a,b], 
      

    ∂u(x,t)/∂x|_a=∂u(x,t)/∂x|_b   ,   u(b)-u(a)/b-a = ∂u(x,t)/∂x|_b.  


Where we represent delta measures by the closest base element of the finite element approximation. Note that (<ref>) contains a normalization factor, since now the elements integral are not one (since they are not delta measures).

Now we use the Galerkin method of mean weighted residuals for the spatial part of the problem choosing low order elements ϕ_i. This formulation can be found in <cit.>. For our numerical studies we solve the temporal part of the problem (element coefficients) using the implicit-Euler Galerkin Discretization <cit.>, thus the problem is reduced to solve a linear system iteratively for every timestep Δ t. 

In order to implement the previous formulation numerically, we use  to do all the calculations for the simulation. The code is available publicly in <cit.>. There we start by generating a list of {b_n}_n=1^N=500   uniformly distributed points in theinterval  [0,10]. This points are located in the closest interval of the spatial FEM partition { (0+ (n-1)/500,n/500) }_n=1^5000. The histogram of this points, Figure <ref>  can be seen for instance as the number of times a certain criminal act was informed in a zone from the avenue. See Figure <ref>.

 



If we represent this data as an initial condition (<ref>) we obtain the Figure <ref>. Where we plotted alternatively each consecutive FEM basis function red and black.

Now  we solve numerically the problem using the implicit-Euler Galerkin discretization for the problem (<ref>),(<ref>),(<ref>) and we evolve the solution until time t=0.1 using either Neumann boundary conditions, see Figure <ref> and  mean conserving boundary conditions, see Figure <ref>. The solution for the mean conserving boundary condition is positive for this numerical experiment, see Figure <ref>, this fact  is currently being explored for future analytical studies.

As the Figure <ref> shows, the solutions are similar and therefore we can see that for this example the new boundary condition does not generate a noticeable change on the generation of the continuous density distribution. Nevertheless we present the plots of  change of mass Δ m(t) = m(0) - m(t), Figures <ref>, <ref> and change of mean Δμ(t) = μ(0) - μ(t), Figures <ref>, <ref> for both Neumann and mean conserving boundary conditions.








Figures <ref> and <ref> present the real difference in the evolution of the density. We effectively see that the mean conserving boundary conditions conserve the mean in the density estimation process. On the other hand if we where to have an initial condition that is biased to one of the boundaries, the differences of the estimated densities by both boundary conditions would differ significantly. However there is no evidence to think that this phenomena occurs in real avenues.

For the numerical experiment presented here we can see that the mean for the Neumann boundary conditions has changed about 0.4% in t = 0.1. This change is small, in fact, for an avenue of 10 km, the change in mean would be about 40 m.  We conclude that for this numerical experiment  for the process of density estimation (when the data has not change to much due to the smoothing process) the Neumann boundary condition provide a very fast (since they are easy to implement) and  accurate way to estimate a continuous probability density.  Nevertheless the mean of the sample is not preserved exactly, on the other hand, the mean conserving boundary condition, apart from being also easily implementable, is accurate and do preserve the mean of the sample.



§ TWO-DIMENSIONAL DENSITIES


We now present the problem for the diffusion equation <cit.> in two dimensions,


    ∂ u( x,t)/∂ t- ∇^2 u( x,t) =0,      x=(x_1,x_2) ∈Ω⊂ℝ^2 , t>0.


Again we want the conservation of mass and mean in the time evolution of the density. Consider first the conservation of the total mass. We have,


    d/dt( ∫_Ω u( x,t) d x)  =  ∫_Ω∇^2 u( x,t)  d x
     = ∫_∂Ω∂ u( x,t)/∂ν  dσ,


where ∇ u ·ν= ∂ u( x,t)/∂ν, and ν denotes  the outward normal unit vector to ∂Ω.   To deduce this relation we used (<ref>), and the first Green identity <cit.>. Consider now the conservation of the mean. We have,


    d/dt( ∫_Ω x_i u( x,t) d x)     =   ∫_Ω  x_i ∇^2 u( x,t)  d x
        =∫_∂Ω     x_i ∂ u( x,t)/∂ν  dσ - ∫_Ω∇_i u( x,t) d x,


where ∇_i u( x,t) =   e_i  ·∇ u( x,t), assuming Cartesian unit vectors. Again (<ref>) and the first Green's identity   were used to obtain the desired result.

Then the conditions that we have to impose on u( x,t) in order to conserve mean and mass are:


    ∫_∂Ω     x_i ∂ u( x,t)/∂ν  dσ = ∫_Ω∇_i u( x,t) d x  ,    i = 1,2,     and    ∫_∂Ω∂ u( x,t)/∂ν  dσ =0.


The advantage of two dimensional domains is that we are not restricted to impose only two conditions for the boundary(mean and mass conservation). For these domains we can in principle conserve additional higher moments of the density distribution that are meaningful for the particular problem. Applications on two dimensional domains are of special interest for the author since a two dimensional map of the city can generate really robust results in the field of  Perception of security(PoS). 



§ CONCLUSIONS


The proposed mean conserving boundary conditions were shown to effectively maintain the mean of the initial data sample over the continuous density estimation process. This was also confirmed by the numerical simulation of the estimation process where we used a list of uniformly distributed points in the interval [0,10] as an initial condition.

The numerical experiments presented here show that even though Neumann boundary conditions do not conserve the mean over time, they are accurate enough to maintain the mean in a very restricted interval before the over-smoothing of the density estimation process. 

We showed the application and some of the consequences of both the idea of (KDE) and the new boundary conditions to avenues in a city. The consequences of implementing the diffusion equation with the proposed boundary conditions in companion of more special initial conditions and in 2D domains remains to be analyzed. 



§ ACKNOWLEDGMENTS


I would like to express my gratitude to  Juan Galvis, whose guidance was essential for the completion of this manuscript. I also want to thank Francisco A. Gómez and Zdravko Botev, whose comments were really appreciated for the analysis of the results.

9
  
 
  bib1 Silverman, B. W. (1986). Density Estimation for Statistics and Data Analysis. Chapman and Hall, London.
  bib2 Wand, M. P. and Jones, M. C. (1995). Kernel Smoothing. Chapman and Hall, London.
  
  bib3 Baxter, M. J., Beardah, C. C. and Westwood, S.  (2000). Sample size and related issues in the analysis of lead isotopedata. J. Archaeological Science 27 973-980.
  
  bib4   DiNardo, J., Fortin, N. M., and Lemieux, T. (1996). Labor market institutions and the distribution of wages, 1973 1992: A semiparametric approach. Econometrica 64 1001- 1044.
  
  bib0 Gómez, F., Torres, A. , Galvis, J., Camargo, J. and Martínez, O.  (2016). Hotspot mapping for Perception of Security. Smart Cities Conference (ISC2), 2016 IEEE International
  
  bib5 Botev, Z. I. 
, Grotowski, J. F.  and Kroese, D. P.  (2010). Kernel density estimation via Diffusion. Ann. Statist, Vol. 38, No. 5, 2916–2957.
  
  
  
  bib6 Chaudhuri, P. and Marron, J. S. (2000). Scale space view of of curve estimation. Ann. Statist. 28 408–428. MR1790003.
  
  bib7 Olver, P.  (2014). Introduction to Partial Differential Equations. Springer-Verlag, 1st Ed.
  
  bib8  Thomee, V. (1997). Galerkin Finite Element Methods for Parabolic Problems. Springer-Verlag, 2nd Ed.
  
  bib9 1D Inhomogenious Diffusion Equation solution with FEM. Patarroyo K. available online. <http://nbviewer.jupyter.org/github/MrKeithPatarroyo/Example/blob/master/FEM-Public/1D%20Heat%20Equation%20FEM-Public.ipynb> 
   
   bib10 Salsa, S. (2008). Partial Differential Equations in Action From Modelling to Theory. Springer-Verlag Italia, Milano.
   


