Combinación lineal de variables discretas $T_i$ con $P(T_i=1)=P(T_i=-1)=1/2$

Question

Dejemos que $T_1,...,T_n$ sea iid con una distribución Rademacher; es decir , $P(T_i=1)=P(T_i=-1)=1/2$ ; y dejar que $w = (w_1,...,w_n) \in \mathbb{R}^n$ sin más limitaciones en $w$ .

¿Existe una forma de calcular $P(T_1\sum_{i=1}^n T_i w_i \geq 0)$ sin pasar por todos los $2^n$ posibles resultados de la $T$ 's?

Answer 1

Sí.

Empecemos por simplificar la cuestión. El evento

$$T_1\sum_{i=1}^n T_i w_i \geq 0$$

es la unión de los eventos disjuntos (determinados por $T_1=\pm 1$ )

$$w_1 + \sum_{i=2}^n T_i w_i \geq 0,\ -w_1 + \sum_{i=2}^n T_i w_i \geq 0.$$

Tras una sencilla manipulación algebraica, y utilizando el hecho de que cada $T_i$ tiene la misma distribución que $-T_i$ (es "simétrico"), estos eventos son ambos equivalentes a los de la forma

$$\sum_{i=2}^n T_i w_i \leq w$$

para $w=\pm w_1$ . Aprovechando de nuevo las simetrías del $T_i$ (para que todos los $w_i$ no negativo) y reordenándolos, se puede escribir

$$\sum_{i=2}^n T_i |w_i| \leq w$$

donde $|w_2| \ge |w_3| \ge \cdots \ge |w_n|$ .

Necesitamos calcular las probabilidades de que se produzcan dos eventos de este tipo y deseamos hacerlo de forma eficiente. Asumiré que la probabilidad debe calcularse con gran precisión, de modo que las aproximaciones serían inaceptables. (Para muchas configuraciones del $|w_i|$ De lo contrario, las aproximaciones normales a las distribuciones podrían funcionar bien. Se puede realizar un análisis más profundo utilizando la función característica $\phi$ de $\sum_{i=2}^n T_i |w_i|$ que por definición es $$\phi(t) = \prod_{i=2}^n\left(\frac{\exp(-j |w_i|t)}{2} + \frac{\exp(j |w_i|t)}{2}\right)=\prod_{i=2}^n\cos(w_it);$$ $j = \sqrt{-1}$ .)

---

A partir de ahora, reduzca $n$ por $1$ , inicie la indexación con $i=1$ y asumir todos los $w_i$ son positivos. (Obviamente, cualquier valor cero puede ser ignorado). $1 \le k \lt n$ . Considere la distribución de $X = \sum_{i=1}^n T_i w_i$ condicionada a la primera $k$ valores de $T_i$ :

$$\Pr(X \le w\,|\, T_1, \ldots, T_k) = {\Pr}_{(T_{k+1},\ldots,T_n)}\left(\sum_{i=k+1}^n T_i w_i \le w - \sum_{i=1}^k T_i w_i = w_0\right).$$

Desde $|T_i| \le 1$ la suma de la izquierda está acotada:

$$-u_{k+1} = -(w_{k+1} + \cdots + w_n) \le \sum_{i=k+1}^n T_i w_i \le w_{k+1} + \cdots + w_n = u_{k+1}.$$

Dejemos que $u_{k+1}$ sea el lado derecho: es una suma acumulativa de los $w_i$ acumulado a partir de los valores más bajos de la derecha. Obviamente, ahora si $u_{k+1} \le w_0$ entonces $X \le w$ es seguro; y si $-u_{k+1} \gt w_0$ entonces $X \le w$ tiene una probabilidad nula. Esto lleva a una simple rama y límite algoritmo, porque no necesitamos buscar más para evaluar la distribución. Mientras que una enumeración exhaustiva habría tenido que examinar $2^{n-k}$ casos posibles, hemos hecho una determinación de su contribución a la distribución en $O(1)$ tiempo.

No hay muchas esperanzas de que esta mejora conduzca a una el peor de los casos algoritmo que es mejor que $O(2^n)$ en el rendimiento (aunque creo que se puede reducir a $O(2^{n/2})$ que -aunque mucho mejor- sigue sin ser polinómico). Sin embargo, la mejora es lo suficientemente buena como para que valga la pena considerarla, especialmente para valores moderados de $n$ donde la enumeración exhaustiva empieza a ser impracticable (en algún punto por encima de $20$ y ciertamente por encima de $40$ ). Por lo tanto, pasemos a evaluar el algoritmo. ¿Cuál es su eficacia?

El peor caso es cuando la mayoría de los $w_i$ tienen tamaños comparables, ya que entonces la heurística branch-and-bound rara vez logra algo. Afortunadamente, esta es exactamente la situación en la que el Teorema Central del Límite puede proporcionar excelentes aproximaciones. Por lo tanto, quizás sea aún más interesante explorar situaciones en las que los tamaños de las $w_i$ están muy repartidos.

Para ello, he creado cuatro conjuntos de datos con $n=15$ de cuatro distribuciones con formas radicalmente diferentes: Beta $(1/5,1/5)$ (fuertemente bimodal), Exponencial (modo alto cerca de cero), Uniforme y Gamma $(20)$ (casi Normal, con valores muy cercanos a $20$ ). Cada conjunto de datos se normalizó a un máximo de $1$ . Utilizando el método de rama y límite, $\Pr(X \le w)$ se calculó para $19$ valores de $w$ que van desde $0$ hasta $2\sqrt{n}$ . (Los valores negativos de $w$ no es necesario mostrarlo porque las distribuciones de $X$ son casi simétricos). La figura muestra los resultados, graficando las probabilidades (la función de distribución acumulativa empírica de $X$ ) en la fila superior y el eficiencia (en escalas log-lineales) en la fila inferior. La "eficiencia" es la relación entre el número de configuraciones necesarias para la enumeración exhaustiva ( $2^n$ ) al número de configuraciones consideradas por el algoritmo branch-and-bound.

Para la distribución Beta, cuya fuerte bimodalidad limita esencialmente los cálculos a la mitad superior de los datos, las eficiencias son mayores. Como era de esperar, son menores para la distribución Gamma, cuyos valores son todos comparables. Sin embargo, incluso en este caso difícil, todas las eficiencias superan $3$ .

Las eficiencias más pequeñas suelen ser para $w=0$ . Con el tiempo, la eficiencia aumentará a medida que $|w|$ aumenta. En conjunto, estos resultados son una fuerte evidencia de que el enfoque descrito aquí no sólo es computacionalmente más eficiente que la enumeración exhaustiva, sino que tiende a ser mucho más eficiente.

---

La aplicación del algoritmo es directa y sencilla. R El código es el siguiente. Incluye una sección que prueba el algoritmo (comparando su salida con una enumeración exhaustiva) y otra sección para reproducir la figura.

#
# The algorithm.
#
f <- function(w0, w) {
  w <- sort(abs(w), decreasing=TRUE)
  w.sum <- c(rev(cumsum(rev(w)))[-1], 0)
  count <- 0                              # Counts calls to f()
  f <- function(u, u.sum, w0) {
    count <<- count + 1
    y <- sapply(w0 + c(-1,1)*u[1], function(w1) {
      if (w1 < -u.sum[1]) return(0)
      if (w1 >= u.sum[1]) return(1)
      return(f(u[-1], u.sum[-1], w1))
    })
    return(mean(y)) # (This could easily be changed to accommodate 
                    # other probabilities for the T[i].)
  }
  list(Value=f(w, w.sum, w0), Count=count)
}
#
# Test with complete enumeration.  The plot pairs should exactly coincide.
#
binary <- function(a, b, zero=-1, one=1) rep(c(rep(zero, 2^a), rep(one, 2^a)), 2^b)
n <- 9
b <- sapply(1:n, function(a) binary(n-a, a-1))
par(mfrow=c(2,2))
set.seed(17)
for (i in 1:4) {
  w <- rexp(n)
  x <- b %*% w
  y <- sapply(x, function(w0) f(w0, w)$Value) #$
  plot(ecdf(x))
  points(x, y, pch=16, cex=1/2, col="Red")
}
#
# Explore the efficiencies actually achieved.
#
n <- 15
qb <- function(q) qbeta(q, 1/5, 1/5); sigma <- sqrt(1/(1+2*1/5))/2
qg <- function(q) qgamma(q, 20)
dist <- list(Gamma=qg, Normal=qnorm, Uniform=qunif, Exponential=qexp, Beta=qb)
par(mfcol=c(2,4))
for (s in c("Beta", "Exponential", "Uniform", "Gamma")) {
  d <- dist[[s]]
  w <- d(((1:n)-1/2)/n)
  w <- w / max(abs(w))
  print(c(system.time({
    y <- sapply(x <- seq(0, 2*sqrt(n), length.out=19), 
                function(w0) {x <- f(w0, w); c(x$Value, x$Count)})
  })["elapsed"], count=sum(y[2, ])))

  plot(x, y[1, ], ylim=c(1/2, 1), type="b", ylab="Probability", main=s)
  plot(x, 2^n/y[2, ], ylim=c(1, 2^n), type="b", log="y", ylab="Efficiency")
}

Answer 2

$$T_1\sum_{i=1}^n T_i w_i = w_1 + T_1\sum_{i=2}^n T_i w_i$$

desde $T_1^2 = 1$ .

Escribe la suma como $S_{2w}$ . Así que estamos viendo la probabilidad

$$P\left(T_1S_{2w} \geq -w_1\right) = \\P\left(T_1S_{2w} \geq -w_1 \mid T_1 =1\right)\cdot P(T_1 =1) + P\left(T_1S_{2w} \geq -w_1 \mid T_1 =-1\right)\cdot P(T_1 =-1)$$

$$\Rightarrow P\left(T_1S_{2w} \geq -w_1\right)=\frac 12 P\left(S_{2w} \geq -w_1 \right)+ \frac 12P\left(-S_{2w} \geq -w_1 \right)$$

$$=\frac 12 \Big[1-P\left(S_{2w} \leq -w_1 \right)+P\left(S_{2w} \leq w_1 \right)\Big] \tag{1}$$

Ahora defina $Z_i \equiv (T_i + 1)/2$ . Entonces el $Z_i$ son Bernoulli $(1/2)$ variables aleatorias.

Así que

$$S_{2w}=\sum_{i=2}^n (2Z_i-1) w_i = 2\sum_{i=2}^n Z_iw_i - \sum_{i=2}^n w_i \tag{2}$$

Inserción en $(1)$ y reordenando obtenemos

$$P\left(T_1S_{2w} \geq -w_1\right) = \frac 12 \left[1-P\left(\sum_{i=2}^n Z_iw_i \leq (1/2)\sum_{i=2}^n w_i-(w_1/2) \right)+P\left(\sum_{i=2}^n Z_iw_i \leq (1/2)\sum_{i=2}^n w_i+(w_1/2) \right)\right]$$

Así que necesitamos probabilidades relacionadas con $S_{Zw}\equiv \sum_{i=2}^n Z_iw_i$ . El Función generadora de probabilidad , PGF de un Bernoulli $(1/2)$ es

$$G_Z(z) = \frac 12(1+z)$$

y de la combinación lineal de Bernoullis independientes

$$G_{S_{Zw}}(z) = \prod_{i=2}^{n}\left(\frac 12(1+z^{w_i})\right) = \frac 1{2^{n-1}}\prod_{i=2}^{n}(1+z^{w_i})$$

Desde el PGF podemos recuperar las probabilidades de la función de masa de probabilidad y, por tanto, de la función de distribución acumulativa que es nuestro objetivo.

UN EJEMPLO

A efectos de exposición, supongamos que $n=3$ y $(w_1 =2, w_2 =2, w_3 = 4)$ .

Entonces $$(1/2)\sum_{i=2}^n w_i-(w_1/2) = 2, \;\;\; (1/2)\sum_{i=2}^n w_i+(w_1/2) = 4$$

Así que queremos evaluar

$$P\left(T_1S_{2w} \geq -1\right) = \frac 12 \left[1-P\left(S_{Zw} \leq 2 \right)+P\left(S_{Zw} \leq 4 \right)\right]$$

$$=\frac 12 \left[1+P\left(S_{Zw} =3 \right)+P\left(S_{Zw} = 4 \right)\right]$$

En un ejemplo tan simple, no necesitamos realmente el PGF Dado los valores supuestos de los pesos en el Bernoullis, es fácil ver que $P\left(S_{Zw} =3 \right) = 0$ y que $P\left(S_{Zw} = 4\right) = 1/4$ . Así obtenemos la respuesta, para el ejemplo concreto

$$P\left(T_1S_{2w} \geq -1\right) = 5/8$$

Verificación a través del FGP
Para el ejemplo numérico concreto

$$G_{S_{Zw}}(z) = \frac 1{2^{n-1}}\prod_{i=2}^{n}(1+z^{w_i}) = \frac 1{4}(1+z^2)(1+z^4) = \frac 1{4}(1+z^2+z^4+z^6)$$

La relación que une el PGF con probabilidades es ( $(k)$ denotando el orden de la derivada tomada

$$P(S_{Zw} = k) =\frac{ G_{S_{Zw}}^{(k)}(0)}{k!}$$

Tenemos

$$G_{S_{Zw}}^{(1)}(z) = \frac 1{4}(2z+4z^3+6z^5)$$ $$G_{S_{Zw}}^{(2)}(z) = \frac 1{4}(2+12z^2+30z^4)$$ $$G_{S_{Zw}}^{(3)}(z) = \frac 1{4}(24z+120z^3)$$ $$G_{S_{Zw}}^{(4)}(z) = \frac 1{4}(24+360z^2)$$

Así que

$$P\left(S_{Zw} =3 \right) = \frac{ G_{S_{Zw}}^{(3)}(0)}{3!} = 0$$ $$P\left(S_{Zw} =4 \right) = \frac{ G_{S_{Zw}}^{(4)}(0)}{4!} = \frac 1{4}$$

Por lo tanto,

$$P\left(T_1S_{2w} \geq -1\right) = \frac 12 \left[1+0+\frac 14 \right] = 5/8$$

Ahora imagina lo que supondría si el ejemplo no fuera tan sencillo.