Goal

We sample $n$ individuals with replacement and measure their trait values $x_i$ for all individual $i$. The variable $x$ is bounded [0,1].

Let $\mu$ be the average trait value in the population and $\sigma^2$ be the variance in the population.

What is the unbiased estimator for $\frac{\sigma^2}{\mu (1 - \mu)}$?

Issue

Such non-linear statistics can be really hard to solve. I managed to find an unbiased estimators for both the nominator and the denominator but I can't combine them easily as their covariance term is not null.

If you are curious you can find some calculations here where I also considered a case where sampling is without replacement (which is actually my ultimate goal).

Question

What Mathematica tools can I use to calculate / approximate an good estimator for $\frac{\sigma^2}{\mu (1 - \mu)}$?

• Could you expand on your need for an unbiased estimator? Some biased estimators have a bad rap. I'd certainly take a biased estimator over an unbiased estimator if the mean square error was much smaller. – JimB Sep 14 '17 at 23:39
• @JimBaldwin Good point! It is more by ignorance than anything else that I put more interest on accuracy than precision. I edited to write "good" instead of "unbiased" – Remi.b Sep 15 '17 at 0:24

This is an extended comment and certainly not an answer.

You might want to just stick with plugging in the sample estimates for $\sigma^2$ and $\mu$. But to feel more or less comfortable doing so, you should perform some simulations with known/generated populations. Below is some Mathematica code to obtain "nearly exact" estimates of bias and precision based on the sample version of your population statistic. (By "nearly exact" I mean within internal roundoff error.) And the code below only works for small population sizes and small sample sizes.

(* Generate a finite population from a beta distribution *)
nPop = 20; (* Population size *)
population = RandomVariate[BetaDistribution[5, 3], nPop];

(* Population statistic of interest *)
θ = Variance[population]/(Mean[population] (1 - Mean[population]));

(* Generate all-possible samples *)
n = 5;  (* Sample size *)
samples = Subsets[population, {n}];
sampleSpace = (Variance[#]/(Mean[#] (1 - Mean[#])) & /@ samples);

(* Calculate bias *)
bias = Mean[sampleSpace] - θ;

(* Show results *)
Print["Population statistic: " <> ToString[θ]]
Print["Bias: " <> ToString[bias]]
Print["Bias as a percent of the population statistic: " <>
ToString[100 bias/θ]]
Print["Bias as a percent of the standard deviation of the estimates: \
" <> ToString[100 bias/StandardDeviation[sampleSpace]]]

(* Population statistic: 0.0963137 *)
(* Bias: -0.000364909 *)
(* Bias as a percent of the population statistic: -0.378876 *)
(* Bias as a percent of the standard deviation of the estimates: -0.654203 *)

(* Display histogram of sample space *)
Histogram[sampleSpace, Automatic, "PDF"] For this small and made-up example, bias is the least of your worries.