UMVUE of $ \theta^2 (1- \theta ) $ X is random sample from bernoulli distribution Let $ X_1, X_2  ..... X_n $ be a random sample from  bernoulli distribution with parameter $ \theta $ , Obtain   UMVUE of $ \theta^2 (1- \theta ) $
MY APPROACH
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 $
 A: 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))))
}

