I am not sure how to render this as a KDE task for practical purposes. Instead, I will propose an alternative regression type approach.
In my example, I will use a bivariate normal point cloud to get at the rotational symmetry of your problem.
$$
(x,y)\sim N(0,Diag(1))
$$
We then bin the data according to Euclidean coordinates and divide the x and y axis into bins with some equal spacing. For each bin, record the median $x,y,r,\theta$ values, and then the number of observations in each bin for a Poisson regression. In the case of the Multivariate normal I have presented above, we do not need to worry about $\theta$ because the polar representation of the the above is:
$$
f(r,\theta)=\frac{1}{2\pi}r exp(-r^2/2)
$$
This suggests in the log-linear regression, we need to enter $r^2$ as the relevant co-variate. Here is the implementation.
require(mvtnorm)
require(dplyr)
require(ggplot2)
set.seed(1234)
m=rep(0,2);s=diag(rep(0.1,2))
data=data.frame(rmvnorm(10000,mean = m,sigma = s))
colnames(data)<-c('x','y')
data$r<-sqrt(data$x^2+data$y^2)
data$theta<-atan(data$y/data$x)
binLim<-round(max(abs(data[,c('x','y')])))
binsize=0.1
data$x_interval<-findInterval(data$x,seq(-binLim,binLim,binsize))
data$y_interval<-findInterval(data$y,seq(-binLim,binLim,binsize))
binnedData=data%>%group_by(x_interval,y_interval)%>%summarise(
n=n(),x=median(x),y=median(y),r=median(r),theta=median(theta)
)
binnedModel<-glm(n~I(r^2),
data=binnedData,
family=poisson)
binnedData$fitted<-binnedModel$fitted.values
plot(binnedData$fitted,binnedData$n)
binnedData$density_est=binnedData$fitted/sum(binnedData$fitted)
For checking against the known density function in this case: take the median values we saved and compute the integral where $b$ is the binsize used for binning the data. They seem to be roughly correct.
$$
\int_{x-b}^{x+b}\int_{y-b}^{y+b} N(0,\Sigma)
$$
binnedDensities<-apply(binnedData%>%data.frame,1,function(bin){
pmvnorm(lower = c(bin[['x']],bin[['y']])-binsize/2,
upper = c(bin[['x']],bin[['y']])+binsize/2 ,mean = m, sigma = s)
})
ggplot(binnedData,aes(x=x,y=y))+geom_point(aes(colour=binnedDensities),size=5)+
scale_color_viridis_b()+theme_bw()+ggtitle('Analytical Density')
ggplot(binnedData,aes(x=x,y=y,colour=density_est))+geom_point(size=5)+
scale_color_viridis_b()+theme_bw()+ggtitle('Estimated Density')
plot(binnedData$density_est,binnedDensities)
