Su pregunta es cómo encontrar la esperanza de $n\bar X^{k}$ (para $k=3$) cuando $\bar X$ es la media de $n$ variables de Bernoulli independientes e idénticamente distribuidas.
A menos que esté familiarizado con las propiedades de los coeficientes binomiales y sea hábil en álgebra, quizás la forma más sencilla de obtener esto sea a través de la función generadora de probabilidad. Un punto de partida conveniente para este análisis es observar que $Y = n\bar X = X_1 + X_2 + \cdots + X_n$ tiene una distribución binomial con parámetros $(n,\theta).$ Eso sugiere expresar $n\bar X^{k}$ en términos de $Y.$ Un álgebra simple nos dice
$$n\bar X^k = n^{1-k} (n\bar X)^k = n^{1-k}Y^k.$$
Por lo tanto, reconociendo que la constante $n^{1-k}$ se factorizará fuera de la esperanza, hemos reducido el problema a encontrar el $k^\text{th}$ momento de $Y.$
La función generadora de probabilidad de $Y$ se da, por definición, como
$$\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{*}$$
La última igualdad es una aplicación del Teorema Binomial.
Sea $D$ el operador diferencial para $t$ y observe que
$$(tD) p_Y(t) = \sum_{y=0}^n \Pr(Y=y)\, tD(t^y) = \sum_{y=0}^n \Pr(Y=y)\, yt^y.$$
Iterando esto $k$ veces nos dice
$$(tD)^kp_Y(t) = \sum_{y=0}^n \Pr(Y=y)\, y^k t^y.$$
Si, después de realizar este cálculo, establecemos $t=1,$ obtenemos
$$(tD)^k p_Y(t)\bigg|_{t=0} = \sum_{y=0}^n \Pr(Y=y)\, y^k = E\left[Y^k\right],$$
según la definición de esperanza. Calcular estas derivadas usando la expresión en $(*)$ nos dice
$$E\left[Y^k\right] = (tD)^k (1-\theta+\theta t)^n\,\bigg|_{t=1}.$$
Aproveche esto calculando derivadas sucesivas usando las reglas de la suma, el producto y la cadena. Como abreviatura, sea $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}$$
Sustituyendo $t=1$ da
$$\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}$$
No olvide multiplicar estos por $n^{1-k}$ al calcular la esperanza de $n\bar X^k.$ Puede valer la pena verificar su trabajo calculando las esperanzas para $k=1$ y $k=2,$ que ya conoce.
Para apreciar lo general y práctico que es este enfoque, aquí tiene una función R para calcular estos momentos.
#
# Calcular los momentos Binomiales(n, theta) para k = 0, 1, ..., k.max.
# Requiere que `n` y `k.max` sean números naturales y tenga sentido
# solo para 0 <= theta <= 1 (pero funcionará para cualquier `theta`).
#
momentos <- función(n, theta, k.max) {
#
# Aplicar tD a un polinomio a[1] + a[2]*t + a[3]*t^2 + ... + a[n]*t^(n-1).
# El resultado es 0*a[1] + 1*a[2]*t + 2*a[3]*t^2 + ... + (n-1)*a[n]*t^{n-1},
# cuyos coeficientes son el producto componente por componente de `a` y 0:(n-1).
#
tD <- función(a) a * (seq_along(a) - 1)
#
# Evaluar la función generadora de probabilidad y sus derivadas
# de forma iterativa.
#
i <- seq_len(n+1) - 1 # 0,1, ... n
a <- elija(n, i) * theta^i * (1 - theta)^(n-i) # (1-theta+theta*t)^n
c(1, sapply(seq_len(k.max), función(i) sum(a <<- tD(a))))
}
