As much as I know, in the usual machine learning setting, given n data points $(x_1, y_1),...,(x_n, y_n)$ one tries to estimate the function $f$ for which $E(y|x_i) = f(x_i)$, or $E(y) = f(x)$ (pardon me for my lack of mathematical precision here, unless it makes the question more difficult to answer)
There are known ways to come up with a an estimator $\hat{f}$ for $f$. This is simply us trying to estimate the expected value given the features. However, I was recently asked how to come up with an unbiased estimator of the variance, that is: $Var(y|X)=E(y-E(y|X)|X)$. In other words, we need a function $\hat{h}$ that estimates $h(X)=E(Var(y|X))=E(E(y-E(y|X)|X))=E(y^2|X)-E^2(y|X)$
What I tried: the first term, $E(y^2|X)$ is straightforward to predict, as it is just a different target variable with the same features, but I have failed to derive an unbiased estimator of the second term (coming up with an unbiased estimator of y given X and squaring it is just not gonna cut it, due to bias of that estimator if I am not mistaken, and I have no idea how to estimate that bias itself), and I am not even sure this is the right way to go about it.