Variance of 2D binned uniform density I will start off by saying that I am aware of this post, but applying Sheppard's correction in my case would seem to lead to a negative variance; so either I don't understand that post, or I have a different problem.
Consider $n$ points $P = \{ p_1, ..., p_n \}$ sampled from a 2D uniform distribution in the rectangle $(a_x,b_x)\times (a_y,b_y)$, and a regular grid $G$ spanning that rectangle with $m = m_x \times m_y$ bins all with the same area:
$$ v = \frac{b_x-a_x}{m_x}\cdot \frac{b_y-a_y}{m_y} $$
Then, define the following function:
$$ \forall g\in G,\quad \rho(g) = \frac{\mathrm{Count(P\cap g)}}{n} $$
where we count the number of points in $P$ that fall in any single grid cell $g$. What are the mean and variance of $\rho(g)$?

The formula for $\rho$ looks like a mean estimator to me.
The mean and variance of a continuous uniform distribution over a fixed-size interval are known. With this approach, the mean and variance of $\rho$ across $G$ given $P$ are therefore:
$$
\mu_\rho = \frac{v}{mv} = \frac{1}{m} \qquad
\sigma_\rho^2 = \frac{v^2}{12 n}
$$
(Note that the variance term is divided by $n$ because we are talking about the variance of a mean estimator.)
From simulations though, the variance formula above seems to be wrong. If I understand correctly, applying Sheppard's correction in this case would subtract $v^2/12$ from the variance, which would give a negative result.

The correct approach is to interpret $\rho$ as a binomial r.v., in that case:
$$
    \forall g\in G,\quad \mathrm{Count}(P\cap g) 
    \sim B\left( n, \frac{1}{m} \right)
$$
In that case, the density variable $\rho(g)$ has the following mean and variance:
$$
\begin{align}
    \mu_\rho 
    &= \frac{1}{n}\ n . \frac{1}{m} = \frac{1}{m} \\[3mm]
    \sigma^2_\rho 
    &= \frac{1}{n^2}\ n . \frac{1}{m} . \left(1 - \frac{1}{m}\right) 
    = \frac{m - 1}{m^2 n}
\end{align}
$$
Here is a Matlab program to verify that this is correct via simulations:
    % grid definition
    ax = 0;
    bx = 1;
    ay = 0;
    by = 1;
    
    mx = 11;
    my = 19;
    
    m = mx * my;
    rx = linspace( ax, bx, mx+1 );
    ry = linspace( ay, by, my+1 );
    
    % try with increasing number of points
    npts = ceil(exp(linspace( log(1e2), log(1e4), 20 )));
    
    var_emp = zeros(size(npts)); 
    var_ref = zeros(size(npts));
    for k = 1:length(npts)
        n = npts(k);
        ux = ax + (bx-ax) * rand(n,1);
        uy = ay + (by-ay) * rand(n,1);
        
        H = histcounts2( ux, uy, rx, ry );
        
        var_emp(k) = var(H(:)/n);
        var_ref(k) = (m-1) / (m^2 * n);
    end
    
    % plot results
    plot( npts, var_emp, npts, var_ref, 'LineWidth', 2 ); 
    grid on; xlim(npts([1,end]));
    xlabel('#sample');
    ylabel('Estimator variance');
    legend('Empirical','Theory');


