Smoothing histograms with kernel methods I have a problem where I can receive as output, multidimensional counts in "histogram" form. I can also adjust the size of the bins I receive (i.e., many or few bins). I want to smooth the data and thought kernel methods would be a good solution. When I read about kernel methods in the literature, they spoke of histograms as the simplest kernel method (which I understand conceptually), and that is it. But what I want is this problem formalized in a way that looks to smooth histograms specifically. Does anyone know where I can look to find this discussed more formally? 
Also, to put my question differently. I am looking for methods to find the first derivative of a multidimensional empirical cumulative distribution. 
 A: I do not know if you are still interested in solving this issue. If this is the case, I would suggest to use a different approach if it is ok with you...anyway this answer might be useful to somebody else in the future (I hope at least).
An option to smooth multivariate histograms, is to use P-splines and fit the array of counts as suggested in the comment at your question (see this reference for example). P-splines combine B-spline bases and finite difference penalties (see Eilers and Marx, 1991).
To keep the notation simple, suppose that we want to smooth a 2d hist (the method generalizes to more dim as you will read in the reference I mentioned above). Let define $Y$ the 2-dim array of counts collected on a grid of $n_{1} \times n_{2}$ grid of bins. We can model the expected counts as Poisson using tensor product B-splines as:
$$
\mu = E[\mbox{vec}(Y)] = \exp[(B_{2} \otimes B_{1}) \eta]
$$
where $B_{1}$ and $B_{2}$ are $(n_{1} \times m_{1})$ and $(n_{2} \times m_{2})$ B-spline matrices and $\eta$ is a $m_{1}m_{2}$ vector of unknown spline coefficients to be estimated from the data. The optimal spline coefficients can be estimated by maximizing the following penalized log-likelihood (foe example by using IWLS):
$$
\ell = \sum\log(\mu^{\mbox{vec(Y)}} \exp(-\mu)) - \lambda (\eta^{\top}P\eta)
$$
Again, to keep notation simple, I will consider the isotropic 2d penalty
$$
P = \lambda(I_{n_{2}} \otimes \Delta_{1}^{\top} \Delta_{1} + \Delta_{2}^{\top} \Delta_{2} \otimes I_{n_{1}})
$$
where $\Delta$ is the finite difference matrix operator (I will use here 2nd order differences) and the $\lambda$ parameter regulates the amount of smoothness required for the final estimates (please note that here I consider a single $\lambda$...it is possible to have 2 of them in the definition of $P$ above).
Below you will find a small code reproducing an example (I left some comment in the R-code to help readability, I hope). In this case, I will use a "2D Whittaker smoother" which is a special case of P-spline (you can read about the relationship between these two methods in Eilers (2003) for example).
rm(list = ls())
library(colorout)
library(gplots)
library(fields)

# Simulate some data
set.seed(1234)
N = 500
x = 5 + rnorm(N) * 0.5
y = 2 + rnorm(N) * 0.5

# Bi-variate histogram
h = hist2d(data.frame(x, y), nbins = 25, show = FALSE)

# Data handling
x_mid   = h$x
y_mid   = h$y
cnt     = h$counts
cnt_vec = as.vector(cnt)
n       = dim(cnt)[1]
m       = dim(cnt)[2]

# Bases stuffs (Whittaker) & penalties (isotropic)
dd = 2
lm = 1e2 # Fixed lambda - can be selected
E1 = diag(n) 
E2 = diag(m)
Eo = E1 %x% E2
P1 = crossprod(diff(E1, diff = dd))
P2 = crossprod(diff(E2, diff = dd))
Po = lm * diag(m) %x% P1 + lm * P2 %x% diag(n) # isotropic smoothing, can be made anisotropic

# Smoothing
tol = 1e-8
eta = log(cnt_vec + 1)
for (it in 1:50) 
{
    mu = exp(eta)
    z = cnt_vec - mu + mu * eta
    a = solve(t(Eo) %*% (c(mu) * Eo) + Po, t(Eo) %*% z)
    etnew = Eo %*% a
    de = max(abs(etnew - eta))
    cat('Crit', it, de, '\n')
    if(de < tol) break
    eta = etnew
}

cnt_fit = matrix(exp(eta), nc = m, nr = n)

# Plot results
image.plot(x_mid, y_mid, cnt)
contour(x_mid, y_mid, cnt_fit, add = TRUE, lwd = 2, nlevels = 25)


Finally, I hope I got the notation right and that my answer replies (at least partially) the original question.
A: Actually, I would personally  rather use a kernel smoothing approach than searching to improve the histogram. This is because the AMISE error of kernel smoothing is $O\left(n^{-\frac{4}{5}}\right)$ while histogram is $O\left(n^{-\frac{2}{3}}\right)$ "only". More details on kernel smoothing can be found in M.P. Wand and M.C. Jones. "Kernel Smoothing" (1994).
However, the chapter 4 of the book by David W. Scott "Multivariate density estimation" (1992) is on "frequency polygons", a generalization of the histogram. Here is how he motivates it:
"The discontinuities in the histogram limit its usefulness as a graphical tool for multivariate data. The frequency polygon is a continuous density estimator based on the histogram, with some form of linear interpolation."
