With design weights \(w_i = 1/\pi_i\) the Horvitz–Thompson total is \(\hat t_y = \sum_i w_i y_i\) and the estimated population size is \(\hat N = \sum_i w_i\). Their ratio \(\bar y_w = \sum_i w_i y_i / \sum_i w_i\) is the Hajek mean. It is a ratio of two estimated totals rather than a linear statistic, so its variance comes from linearisation rather than directly.
Details
The variance estimator is the standard with-replacement linearisation
of that ratio,
$$v(\bar y_w) = \frac{n}{(n-1)(\sum_i w_i)^2}
\sum_i w_i^2 (y_i - \bar y_w)^2,$$
the residual technique of Sarndal, Swensson & Wretman (1992,
Section 5.6) applied to the Hajek ratio. This is the same quantity
survey::svymean reports for a design declared with
ids = ~1; the identity was checked numerically against that
package when this function was written.
weights are design weights \(1/\pi_i\), not frequencies. Zero
weights drop a unit; negative weights are rejected.
Mirrors morie.fn.wmeansr on the Python side.