A partir de una distribución aleatoria en SQL | de Sami Abboud | febrero de 2024

De una función de densidad de probabilidad a muestras aleatorias

tAquí hay diferentes métodos para actualizar la política de un agente de aprendizaje por refuerzo en cada iteración. Hace unas semanas comenzamos a experimentar reemplazando nuestro método actual con un paso de inferencia bayesiano. Algunas de las cargas de trabajo de datos dentro de nuestro agente están escritas en SQL que se ejecuta en el motor BigQuery de GCP. Usamos esta pila porque proporciona capacidades computacionales escalables, paquetes de aprendizaje automático y una interfaz SQL sencilla.

El paso de inferencia bayesiana que buscamos implementar utiliza uno de los parametrizaciones alternativas de una distribución beta. Lo que significa que necesitaremos poder utilizar una distribución beta en SQL. Mientras trabajaba en esto, descubrí que extraer de la distribución aleatoria en SQL es un tema con muy pocos ejemplos bien documentados. Así que estoy escribiendo sobre ello aquí…

BigQuery no tiene una distribución beta. No tiene la capacidad de extraer cualquier distribución aleatoria. Entonces mi primera intuición fue tomar la definición de la distribución beta, escribirla en SQL, establecer los parámetros usando un CTA, dibujar un número aleatorio entre 0 y 1 y calcular el valor de esa función.

Pero estamos en 2024, así que le pregunté a ChatGPT cómo lo haría:

A mí: “¿​Cómo se crean sorteos aleatorios a partir de una distribución beta en BigQuery?

ChatGPT:

WITH beta_draws AS (
SELECT
(alpha — 1) / (alpha + beta — 2) AS a,
(beta — 1) / (alpha + beta — 2) AS b
FROM
(SELECT 2 AS alpha, 5 AS beta) — Set your desired alpha and beta values
)
SELECT
(a * POW(RAND(), alpha - 1) * POW(1 - RAND(), beta - 1) + b) AS beta_draw
FROM
beta_draws,
UNNEST(GENERATE_ARRAY(1, 1000)) AS rand_num

A mí pensando para mi mismo: Correcto, entonces claramente eso no funcionará.

¿Ves el problema en el código? ChatGPT sortea dos diferente Valores x para el PDF de distribución beta presunta. Arreglé esto, limpié un poco la consulta y probé 1000 valores. Y aquí está el código SQL para hacer eso:

WITH raw_data AS (
SELECT
‘user_a’ AS contact_id
, 0.674 AS probability
, 72 AS lifetime_messages_received
), parameters AS (
SELECT
contact_id
, probability*lifetime_messages_received AS alpha
, (1.0 — probability)*lifetime_messages_received AS beta
, RAND() AS x
FROM raw_data
CROSS JOIN UNNEST(GENERATE_ARRAY(1, 1000)) AS draw_id
)
SELECT
contact_id
, ARRAY_AGG(POW(x, alpha — 1.0) * POW(1.0 — x, beta — 1)) AS beta_x
FROM parameters
GROUP BY contact_id

Gracias a todos, eso es todo 🎁 ¡Nos vemos en la próxima publicación!

¡EQUIVOCADO! 🔴

Tomemos una implementación confiable de extracción de una distribución beta usando los mismos parámetros y comparemos los resultados. He usado SciPy beta.rvs() en Python y aquí hay dos histogramas de 100 bins que permitirán comparar las dos distribuciones dibujadas.

from scipy.stats import beta
alpha_param = 0.674 * 72
beta_param = (1–0.674) * 72
scipy_beta_draw = beta.rvs(alpha_param, beta_param, size=1000)

Bueno, no hace falta una lupa para darse cuenta de que las distribuciones son diferentes. Regresé a la definición de la distribución beta y me di cuenta de que podría deberse a que la distribución beta también tiene una constante de escala que depende de la función gamma que no incluí en el cálculo 🤦.

Problema: La función gamma no tiene expresión en forma cerraday BigQuery no proporciona una implementación que se aproxime a ella. Entonces, en este punto decidí cambiar a Python, un lenguaje con el que estoy más familiarizado y que hará que mi experimentación sea más eficiente. La idea era que si lo concretaba en Python, podré traducirlo a SQL. Todavía necesitaría alguna forma de aproximarme a una función gamma, pero paso a paso.

Implementemos un sorteo manual de una distribución beta en Python, pero ahora con la constante correcta usando la función gamma de SciPy:

import numpy as np
from scipy.special import gamma
from scipy.stats import uniform
alpha_param = 0.674 * 72
beta_param = (1–0.674) * 72
constant = gamma(alpha_param + beta_param) / (gamma(alpha_param) * gamma(beta_param))
scipy_manual_beta_draw = np.array([
constant*pow(x, alpha_param-1)*pow(1-x, beta_param-1)
for x in uniform.rvs(size=1000)
])

Examinemos nuevamente la distribución usando un histograma de 100 contenedores:

Lo primero que notamos es que la escala ahora es diferente, pero la distribución aún se parece a la dibujada en BigQuery.

… algo está mal… es momento de dar un pequeño paseo para pensar 🚶

…

Después de una corta caminata:

¿Qué significa realmente extraer de una distribución aleatoria? Lo que he implementado hasta ahora es un muestreo aleatorio de la función de densidad de probabilidad beta (PDF) y no funcionó.

Así que tuve que buscar algunas clases de estadística.

Aquí hay un par de buenos repasos sobre:

En resumen, la conclusión es que extraer de una variable aleatoria en realidad significa muestreo de la función de distribución acumulativa inversa (CDF)no de la función de densidad de probabilidad (PDF) como estaba haciendo hasta ahora.

Por supuesto 🤦. Mi profesor de probabilidadde quien me acabo de enterar que falleció a causa de una enfermedad en 2020, me habría animado a “repasar los conceptos básicos” en este punto.

De acuerdo. Revisemos el código Python y ahora extraigamos muestras del CDF inverso (que también se llama la función cuantil) de nuestra distribución beta, y compararla con la distribución dibujada usando beta.rvs() de SciPy:

import numpy as np
from scipy.special import gamma
from scipy.stats import uniform, beta
alpha_param = 0.674 * 72
beta_param = (1–0.674) * 72
n_draws = 1000
# Use SciPy RVS for comparison
scipy_beta_draw = beta.rvs(alpha_param, beta_param, size=n_draws)
# Manual beta draw with the help of the SciPy Gamma function
# We start with a discrete analogue of the Beta PDF we wish to draw from.
# This is just sampling from the PDF at fixed intervals but do check out
# this review for a more in-depth treatment of the subject:
# https://jsdajournal.springeropen.com/articles/10.1186/s40488-015-0028-6
# Set the resolution for generating the discrete PDF
n_samples = 1000
# The beta distribution is supported on the range [0, 1], so we set the
# pdf min and max parameters accordingly
pdf_min = 0.0
pdf_max = 1.0
x_span = np.linspace(pdf_min, pdf_max, n_samples)
constant = gamma(alpha_param + beta_param) / (gamma(alpha_param) * gamma(beta_param))
beta_pdf = np.array([
constant * pow(x, alpha_param — 1) * pow(1 — x, beta_param — 1)
for x in x_span
])
# Using the discrete Beta PDF, we now compute a discrete Beta CDF.
# To do that, we integrate the PDF. For each point x, we sum the PDF until
# that point and multiple with the width of each sample.
freq = 1.0 / n_samples
beta_cdf = beta_pdf.cumsum() * freq
def inv(cdf, q):
“””Return inverse CDF for value q using the quantile function”””
return x_span[np.argmin(cdf < q)]
# Finally, we can now draw n_draws from the discrete inverse of CDF, aka
# generate random samples from a beta distribution
manual_beta_draw = np.array([
inv(beta_cdf, x)
for x in uniform.rvs(size=n_draws)
])

*uf* esto se ve mucho mejor:

Ahora que ya hemos extraído muestras de una variable aleatoria, es hora de volver a SQL. En aras de la simplicidad, y debido a que BigQuery no viene fácilmente con una implementación de una función Gamma¹, me basaré en el distribucion logistica (con parámetros a=0 y b=1).

 — The following 3 parameters need to be adjusted based on the support of the
— PDF of the distribution you wish to draw from. This values are set for a logistic
— distribution with a=0 and b=1
DECLARE pdf_min INT64 DEFAULT -10;
DECLARE pdf_max INT64 DEFAULT 10;
DECLARE n_samples INT64 DEFAULT 5000;
DECLARE sampling_step FLOAT64 DEFAULT (pdf_max — pdf_min) / n_samples;
— The number of random draws you wish to perform
DECLARE n_draws INT64 DEFAULT 1000;
WITH pdf AS (
— The discrete sampling of the logistic distribution PDF
SELECT
x
, exp(-x) / pow(1 + exp(-x), 2) AS y — a=0, b=1
FROM UNNEST(GENERATE_ARRAY(pdf_min, pdf_max, sampling_step)) AS x
), cdf AS (
— The discrete CDF
SELECT
x
, SUM(y)
OVER (
ORDER BY x
) * (1.0 / n_samples) AS y
FROM pdf
), random_draws AS (
— Random draws in the range of [0, max(cdf)]
SELECT
RAND() * (SELECT MAX(y) FROM cdf) as q
, draw_id
FROM UNNEST(GENERATE_ARRAY(1, n_draws)) AS draw_id
)
— Calculate the inverse CDF per draw using the quantile function by generating
— and array of the discrete support of the distribution and returning the value
— of the index just before the randomly generated number is larger than the CDF
SELECT
ARRAY_AGG(x ORDER BY x)[OFFSET(SUM(CAST(y < q AS INT64)))] AS x
FROM random_draws
JOIN cdf
ON TRUE
GROUP BY draw_id;

Comparemos ahora las distribuciones de los tres métodos de muestreo:

1. SciPy’s logistic.rvs()
2. Muestrear manualmente el PDF de distribución logística en Python y extraer una muestra aleatoria según el Paso 2 anterior
3. Haciendo lo mismo en SQL

¡Esto me parece un éxito! 💪

Este código SQL anterior muestra la distribución logística, pero debería funcionar en cualquier distribución en la que pueda obtener una representación discreta del PDF al muestrearlo a intervalos consistentes.