Let $ X_1, X_2 ..... X_n $ be a random sample from bernoulli distribution with parameter $ \theta $ , Obtain UMVUE of $ \theta^2 (1- \theta ) $


I calculated that T = $ \sum X_i $ is complete , sufficient statistic for $ \theta $ , so I thought of calculating $ E( n \bar X^2 (1 - n \bar X) ) $

Which is equal to $ E( n \bar X^2 ) - E(n \bar X^3) $ The first term can be obtained using the variance formula , but I dont know how to proceed with the second term .

I used the help of one of the question explained here ,using blackwellisation and I got

$ \frac{t(t-1)(n-t)}{n(n-1)(n-2)}$ as the overall UMVUE , where $ t = \sum X_i $

  • 1
    $\begingroup$ It is the difference of the MVUEs of $\theta^2$ and $\theta^3$, and these have been discussed here before. Please search the site. $\endgroup$ Jun 12, 2022 at 10:54
  • $\begingroup$ Does this answer your question? UMVUE for Bernoulli $\endgroup$
    – Xi'an
    Jun 12, 2022 at 16:47
  • $\begingroup$ @Xi'an yes , I took help of the same question to reach my answer , but I am not sure if my answer is correct , I just needed to confirm that once . $\endgroup$
    – simran
    Jun 13, 2022 at 7:24
  • 1
    $\begingroup$ Since this is a self-study question, you need to provide more details on how you got this estimator by Rao-Blackwellisation since RB is enough of an argument for UMVUEs. $\endgroup$
    – Xi'an
    Jun 13, 2022 at 7:44

1 Answer 1


Your question is how to find the expectation of $n\bar X^{k}$ (for $k=3$) when $\bar X$ is the mean of $n$ independent identically distributed Bernoulli variables.

Unless you are familiar with properties of binomial coefficients and adept at algebra, perhaps the easiest way to obtain this goes through the probability generating function. A convenient starting point for this analysis notes that $Y = n\bar X = X_1 + X_2 + \cdots + X_n$ has a Binomial distribution with parameters $(n,\theta).$ That suggests expressing $n\bar X^{k}$ in terms of $Y.$ Simple algebra tells us

$$n\bar X^k = n^{1-k} (n\bar X)^k = n^{1-k}Y^k.$$

Thus, recognizing that the constant $n^{1-k}$ will factor out of the expectation, we have reduced the problem to finding the $k^\text{th}$ moment of $Y.$

The probability generating function of $Y$ is given, by definition, as

$$\begin{aligned} p_Y(t) = \sum_{y=0}^n \Pr(Y=y) t^y = \sum_{y=0}^n \binom{n}{y} (1-\theta)^y \theta^y\, t^y = \left(1-\theta + \theta t\right)^n. \end{aligned}\tag{*}$$

The last equality is an application of the Binomial Theorem.

Let $D$ be the differential operator for $t$ and observe that

$$(tD) p_Y(t) = \sum_{y=0}^n \Pr(Y=y)\, tD(t^y) = \sum_{y=0}^n \Pr(Y=y)\, yt^y.$$

Iterating this $k$ times tells us

$$(tD)^kp_Y(t) = \sum_{y=0}^n \Pr(Y=y)\, y^k t^y.$$

If, after performing this computation, we set $t=1,$ we obtain

$$(tD)^k p_Y(t)\bigg|_{t=0} = \sum_{y=0}^n \Pr(Y=y)\, y^k = E\left[Y^k\right],$$

according to the definition of expectation. Computing these derivatives using the expression in $(*)$ tells us

$$E\left[Y^k\right] = (tD)^k (1-\theta+\theta t)^n\,\bigg|_{t=1}.$$

Exploit this by computing successive derivatives using the sum, product, and chain rules. As a shorthand, let $f(t) = (1-\theta + \theta t):$

$$\begin{aligned} (tD)^1 (1-\theta+\theta t)^n &= t D(1-\theta+\theta t)^n = nt\theta(1-\theta+\theta t)^{n-1} = n\theta tf(t)^{n-1}; \\ (tD)^2 (1-\theta+\theta t)^n &= (tD) \left(n \theta t f(t)^{n-1}\right) = n t \theta f(t)^{n-1} + n(n-1)t^2\theta^2 f(t)^{n-2}; \\ (tD)^3 (1-\theta+\theta t)^n &= \cdots = n t \theta f(t)^{n-1} + 3n(n-1)t^2\theta^2 f(t)^{n-2} + n(n-1)(n-2)t^3\theta^3 f(t)^{n-3}. \end{aligned}$$

Plugging in $t=1$ gives

$$\begin{aligned} E[Y] &= n\theta;\\E[Y^2] &= n\theta + n(n-1)\theta^2\\ E[Y^3] &= n\theta + 3n(n-1)\theta^2 + n(n-1)(n-2)\theta^3.\end{aligned}$$

Don't forget to multiply these by $n^{1-k}$ when computing the expectation of $n\bar X^k.$ It might be worth checking your work by computing the expectations for $k=1$ and $k=2,$ which you already know.

To appreciate how general and practical this approach is, here is an R function to compute these moments.

# Compute the Binomial(n, theta) moments for k = 0, 1, ..., k.max.
# Requires that `n` and `k.max` be natural numbers and it makes sense
# only for 0 <= theta <= 1 (but will work for any `theta`).
moments <- function(n, theta, k.max) {
  # Apply tD to a polynomial a[1] + a[2]*t + a[3]*t^2 + ... + a[n]*t^(n-1).
  # The result is 0*a[1] + 1*a[2]*t + 2*a[3]*t^2 + ... + (n-1)*a[n]*t^{n-1},
  # whose coefficients are the component wise product of `a` and 0:(n-1).
  tD <- function(a) a * (seq_along(a) - 1)
  # Evaluate the probability generating function and its derivatives
  # iteratively.
  i <- seq_len(n+1) - 1                           # 0,1, ... n
  a <- choose(n, i) * theta^i * (1 - theta)^(n-i) # (1-theta+theta*t)^n
  c(1, sapply(seq_len(k.max), function(i) sum(a <<- tD(a))))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.