Cálculo eficiente de $E\left[\left(1+X_1+\cdots+X_n\right)^{-1}\right]$ con $(X_i)$ Bernoulli independiente con parámetro variable

Question

Supongamos que tenemos las variables aleatorias $X_1, \ldots, X_n$ que tienen distribuciones Bernoulli con las probabilidades (posiblemente diferentes) $p_1, \ldots, p_n$ . Por ejemplo, $X_1$ = 1 con probabilidad $p_1$ y 0 con probabilidad $1-p_1$ . ¿Existe una forma eficiente de calcular $$E\left[\frac{1}{1+\sum_iX_i}\right]$$ en tiempo polinómico en $n$ ? Si no es así, ¿existe una solución aproximada?

Answer 1

Una aproximación es a través de las funciones generadoras. Para cada variable aleatoria no negativa $S$ , $$ E\left(\frac1{1+S}\right)=\int_0^1E(t^S)\mathrm{d}t. $$ Si $S=X_1+\cdots+X_n$ y las variables aleatorias $X_i$ son independientes, $E(t^S)$ es el producto del $E(t^{X_i})$ . Si además $X_i$ es Bernoulli $p_i$ , $$ E\left(\frac1{1+S}\right)=\int_0^1\prod_{i=1}^n(1-p_i+p_it)\mathrm{d}t. $$ Esta es una fórmula exacta. No sé cuál es la mejor manera de utilizarla para calcular el LHS de forma eficiente. Por supuesto, se puede desarrollar el integrando en el lado derecho, obteniendo una suma de $2^n$ términos indexados por los subconjuntos $I$ de $\{1,2,\ldots,n\}$ como $$ E\left(\frac1{1+S}\right)=\sum_I\frac1{|I|+1}\prod_{i\in I}p_i\cdot\prod_{j\notin I}(1-p_j). $$ Pero podría ser más útil notar que $$ \prod_{i=1}^n(1-p_i+p_it)=\sum_{k=0}^n(-1)^k\sigma_k(\mathbf{p})(1-t)^k, $$ donde $\sigma_0(\mathbf{p})=1$ y $(\sigma_k(\mathbf{p}))_{1\le k\le n}$ son los polinomios simétricos de la familia $\mathbf{p}=\{p_i\}_i$ . Integrar con respecto a $t$ se obtiene $$ E\left(\frac1{1+S}\right)=\sum_{k=0}^n(-1)^k\frac{\sigma_k(\mathbf{p})}{k+1}. $$ La carga computacional se reduce a la determinación de la secuencia $(\sigma_k(\mathbf{p}))_{1\le k\le n}$ .

Nota 1 La última fórmula es una versión integrada de la identidad algebraica que afirma que, para toda familia $\mathbf{x}=\{x_i\}_i$ de ceros y unos, $$ \frac1{1+\sigma_1(\mathbf{x})}=\sum_{k\ge0}(-1)^k\frac{\sigma_k(\mathbf{x})}{k+1}, $$ truncado en $k=n$ ya que, como máximo $n$ valores de $x_i$ son distintos de cero, $\sigma_k(\mathbf{x})=0$ por cada $k\ge n+1$ . Para demostrar la identidad algebraica, observe que, para cada $k\ge0$ , $$ \sigma_1(\mathbf{x})\sigma_k(\mathbf{x})=k\sigma_k(\mathbf{x})+(k+1)\sigma_{k+1}(\mathbf{x}), $$ y calcular el producto de $1+\sigma_1(\mathbf{x})$ por la serie en el RHS. Para aplicar esta identidad a nuestro entorno, introduzca $\mathbf{X}=\{X_i\}_i$ y observe que, para cada $k\ge0$ , $$ E(\sigma_k(\mathbf{X}))=\sigma_k(\mathbf{p}). $$ Nota 2 Más generalmente, para cada número complejo adecuado $z$ , $$ \frac1{z+\sigma_1(\mathbf{x})}=\sum_{k\ge0}(-1)^k\frac{\Gamma(k+1)\Gamma(z)}{\Gamma(k+1+z)}\sigma_k(\mathbf{x}). $$

Answer 2

En cuanto a las aproximaciones para N grandes (digamos, $N>10$ ) , aunque quizás no sea esto lo que buscas... Un enfoque general para aproximar la expectativa de una función de una variable aleatoria $y=g(x)$ es hacer una expansión de Taylor alrededor de la media (de $x$ ) y la toma de expectativas.

$ \displaystyle y \approx g(\mu_x) + g'(\mu_x) (x-\mu_x) + \frac{1}{2}g''(\mu_x) (x-\mu_x)^2 +\cdots \Rightarrow$

$ \displaystyle \mu_y \approx g(\mu_x) + \frac{1}{2!}g''(\mu_x) \; m_{2,x} + \frac{1}{3!} g'''(\mu_x) \; m_{3,x} \cdots $

donde $m_{k,x}$ es el momento k-ésimo centrado de $x$

En nuestro caso, la aproximación de segundo orden da

$ \displaystyle E(y)= E \left[ \frac{1}{1+\sum x} \right] \approx y_0 + y_0^3 \sum \sigma^2_i$

donde $ \displaystyle y_0 = \frac{1}{1 + \sum E(x_i)} = \frac{1}{1 + \sum p_i}$

(aproximación de primer orden) y

$\sigma^2_i = p_i (1-p_i)$

Experimentalmente, obtengo un error relativo que raramente excede de 0,01, con un error aleatorio (uniforme) $p_i$ y $N=15$ . Con $N=30$ es alrededor de 0,001

Answer 3

Si tienes valores numéricos para las p's y quieres un valor numérico para la expectativa, una aproximación simple O(n*n) es calcular la pdf de la suma. Si S_k es la suma de los primeros k de las X, entonces S_k toma valores en {0..k}; si su pdf se mantiene en el array c, entonces el pdf de S_k+1 puede calcularse en c así (donde p es el parámetro para X_k+1):

c[k+1] = p*c[k]
for j=k .. 1
     c[j] = p*c[j] + (1-p)*c[j-1]
c[0] = (1-p)*c[0]

Un pequeño programa en C basado en esto tarda unos 3,25 segundos (en un PC normal) en calcular la expectativa para n = 30000