# Splines free of Gibbs phenomenon for smoothing

I'm trying to do 3-dimensional smoothing of some observations we made. To simplify things, here's a visualisation of one slice of the data we have:

Note that the observations are gridded and contain a lot of zeros. When I run a cubic smoothing spline (I use the package csaps in python) we get something like the plot below.

(My apologies for the changing colour scales).

As can be seen on the second plot, the splines oscillate and become even negative, which in our application does not make sense and is an issue. I researched a bit and found that this phenomenon is often referred to as Gibbs phenomenon that can affect splines (for example in this or that application).

I have a vague understanding of the math behind it and am more interested in the practical application. Does anybody know of an implementation of splines that are free of the Gibbs phenomenon in Python (or R or something else really) that we could use for our application?

Some more background on the application:

We are estimating the probability that two particles are in the same volume. The observations are occupation rates per volume (in some arbitrary unit, such as hours/year). We would like to fit the smooth to the data in an analogous way as you can fit a function to a 1D histogram (in our case it’s just 3D). Since the observed data is noisy (and we don't have access to the underlying generating process), we would like to use the fitted function to draw or to estimate some probabilities instead of using the observations directly. The probabilities (or occupancy rate, which is somewhat similar in our case) are then used to estimate the probability that two particles are in the same volume (the other particle is generated by a different process which we know how to model).

• What is the purpose of the smoothing? Ultimately, to recommend any solution, we really would need to know how to compare one solution to another to assess which might be better.
– whuber
Commented Jan 17, 2022 at 20:02
• The value in each cell is the count of observations we have for the volume. Since the data is somewhat noisy, we would like to smooth it. Commented Jan 17, 2022 at 20:12
• Of course--but that tells us nothing useful for recommending smooths! We need to know what you will be doing with those smooths. In many applications, this Gibbs phenomenon will be harmless. If you wish to reduce it, you have to give up some desirable property of the cubic smoother. You need to give us some information to guess how you want to make such trade-offs.
– whuber
Commented Jan 17, 2022 at 20:14
• why not try KDE, e.g. with a Gaussian kernel? Commented Jan 17, 2022 at 21:25
• Spitballing, (as it looks like the underlying data aren't smooth) but could you model the data as Poisson with a log link & a 2d thin plate spline of the x,y coordinates? Essentially you want the resulting smooth to approach zero (but be slightly none zero) in the large region with all the 0s. Isn't this an issue of the smoother having some amount of overall wiggliness that is needed to get close to the data, while over most of the domain you really want wiggliness to be 0 (or very very close to it) so the fit is ~constant there? This implies an adaptive smooth might work also Commented Jan 19, 2022 at 10:15

Thanks to the hints in the comments, I managed to find a suitable solution to my problem. More precisely, I found two that overcome the issue with the splines. But neither of my solutions uses splines anymore.

I think it's worth highlighting that the suggestion of @Gavin is probably also possible, but I didn't try this since I'm quite happy with what I have so far.

Gaussian Filter

I applied a 3D gaussian filter to the data and the results look reasonable, as shown in the figure below.

From my point of view, the drawback of this method is that the choice of sigma (standard deviation of the Gaussian kernel) is somewhat arbitrary. Instead, I settled on the next method.

Kernel Density Estimation

Thanks to the hint from @Aksakal I looked into kernel density estimation (KDE) and it seems to me suitable for the problem at hand. The result is shown in the figure below.

What I like about this method is that the bandwidth of the kernel can be estimated from the data with cross validation (in my opinion this is an advantage over just simply going for the value that seems to be looking fine).

If anybody is interested in KDE in Python, an excellent comparison of the different implementations to be found here.