Una introducción a Pymc y el lenguaje para describir modelos estadísticos

Introducción a Pymc y lenguaje para modelos estadísticos

Foto de Paul Nicklen en Treehugger

En nuestro artículo anterior sobre por qué la mayoría de los ejemplos de inferencia bayesiana tergiversan lo que es, aclaramos un malentendido común entre los principiantes de estadística bayesiana. Es decir, el campo de la estadística bayesiana NO SE define por su uso del teorema de Bayes, sino por su uso de distribuciones de probabilidad para caracterizar la incertidumbre y considerar el rango completo de posibles resultados. Entonces, por ejemplo, en lugar de que nos digan que un dispositivo médico dado es 95% efectivo para detectar una enfermedad dado el hecho de que sabemos con certeza absoluta que ya estamos infectados (lo que podemos llamar nuestra Tasa de Verdaderos Positivos (TPR)), podemos considerar escenarios donde ese dispositivo solo es 24%, 69% o 91% efectivo para detectar la enfermedad y cómo esos números cambian la probabilidad de contraerla dado que has dado positivo en la prueba del dispositivo. Al rechazar las estimaciones puntuales (como la TPR estática) en nuestros modelos y en cambio favorecer un enfoque probabilístico al considerar un rango completo de resultados, nos mantenemos mucho más fieles a la visión bayesiana de la probabilidad.

Dicho esto, la siguiente fase de nuestro aprendizaje hacia el logro de un estado mental bayesiano al pensar en la teoría de la probabilidad es construir modelos estadísticos bayesianos que integren su enfoque probabilístico para caracterizar la incertidumbre. Una de las mejores herramientas para lograr esto es el paquete de Python PyMC, que se construyó específicamente para realizar inferencia bayesiana y construir modelos probabilísticos de aprendizaje automático.

Sin embargo, antes de comenzar a construir modelos bayesianos, necesitamos hacer un pequeño desvío e introducir el lenguaje (es decir, las notaciones matemáticas) utilizado por los profesionales para formular estos modelos. Comprender el lenguaje de construcción de modelos es práctico por dos razones:

Razón 1. La primera razón es que aprender cómo los profesionales, como científicos e investigadores, formulan modelos utilizando la notación de modelo de estadísticas aplicadas se relaciona directamente con cómo escribimos código ejecutable en PyMC. El uso de notaciones permite al modelador estadístico comunicar sus suposiciones subyacentes sobre un modelo utilizando unas pocas líneas cortas de notación. Sostenemos que las notaciones son una forma de comunicación mucho más simple que la alternativa, que pone la responsabilidad en la audiencia de recordar palabras complicadas como homocedasticidad para comprender las características del modelo con el que están trabajando.

Razón 2. La otra razón es que este lenguaje abarca todos los ámbitos, desde la historia del arte hasta la astronomía y la biología de la conservación (McElreath, 2020). Y, al estar familiarizados con el lenguaje de describir modelos, ampliaremos inherentemente los límites en los que podemos comprender y comunicar conocimientos entre la comunidad científica (McElreath, 2020).

Como introducción general a la notación de modelo de estadísticas aplicadas, el Prof. McElreath (2020) describe algunos principios como base para describir y comprender modelos en el lenguaje coloquial científico:

  1. Los profesionales reconocen un conjunto de variables con las que trabajan, que a veces son observables. Estas se llaman datos.a. Los fenómenos no observables que resultan de los datos, como los promedios, se definen como parámetros.
  2. Los profesionales definen cada variable ya sea en términos de otras variables O en términos de una distribución de probabilidad.
  3. La combinación de variables y sus distribuciones de probabilidad define un modelo generativo conjunto, que podemos usar para simular una observación hipotética y analizar las reales.

Ahora que tenemos la justificación detrás del mundo de las notaciones de modelos, ilustremos cómo las usamos con un ejemplo sencillo. Digamos que eres un conservacionista interesado en cuantificar la probabilidad de encontrarte con un lobo costero de la isla de Vancouver (Canis lupus crassodon) macho en estado salvaje en lugar de una hembra.

Como nota al margen, una de las razones por las que nos interesa utilizar el lobo costero de la isla de Vancouver en nuestro ejemplo es porque se descubrieron recientemente y se clasificaron como una subespecie separada. Por lo tanto, la cantidad de datos e información sobre ellos es limitada en comparación con otras subespecies de lobos, por lo que sería interesante generar estimaciones sobre las características de esta subpoblación aquí y en futuros artículos. No solo los datos son limitados, sino que, en nuestra opinión, son animales muy interesantes y únicos en comparación con otras subespecies de lobos, ya que pueden nadar largas distancias e incorporar alimentos marinos en su dieta natural (Muñoz-Fuentes et al., 2009). Desde una perspectiva más amplia, la investigación sobre ecología animal y estadísticas de población puede tener implicaciones en el estado de conservación de un animal, la gestión de poblaciones y la comprensión de su papel dentro del ecosistema más amplio de una área.

Volviendo a nuestro ejemplo, basándonos solo en esta simple pregunta de investigación, ahora podemos generar suposiciones y construir un modelo que cuantifique la proporción de género entre los lobos costeros de la Isla de Vancouver. Para empezar, podemos suponer que la distribución de los resultados en este experimento mental es binomial, lo que significa que solo puede haber dos resultados que pueden ocurrir como resultado de nuestra pregunta de investigación. O bien nos encontramos con un lobo costero macho durante nuestros viajes por la Isla de Vancouver o una hembra. De hecho, podemos simular cómo esperamos que se vea una distribución binomial cuando hay un 50% de probabilidad de que ocurra cualquiera de los dos resultados binomiales utilizando el siguiente código:

# Importando nuestras bibliotecas requeridasimportar arviz como azimportar matplotlib.pyplot como pltimportar numpy como npimportar pandas como pdimportar pymc como pmimportar scipy.stats como statsaz.style.use("arviz-darkgrid")

x_binom = np.random.binomial(n=10, p=0.5, size=100)sns.histplot(x_binom, kde=True, discrete=True)plt.xlabel(“Número de veces que la moneda cayó en cara (N=10)”)plt.ylabel(“Número de muestras”)plt.title('Visualizando una Distribución Binomial')

En el gráfico, encontraremos una aproximación cercana a una distribución normal a partir del resultado equivalente a 10 lanzamientos de monedas, con el eje x describiendo el número de veces que la moneda cayó en cara. Y si volvemos a ejecutar este experimento 100 veces como hemos hecho en el código, el eje y detalla con qué frecuencia obtenemos el resultado en el eje x. En otras palabras, de las 100 veces que volvimos a ejecutar este experimento, aproximadamente 30 de esos experimentos resultaron en que la moneda cayera en cara alrededor de 5 de cada 10 lanzamientos. Ten en cuenta que los resultados del gráfico variarán ligeramente cada vez que ejecutes el código anterior.

También vamos a incluir otra suposición en nuestro modelo donde asumiremos que cualquiera de los dos resultados tiene una oportunidad igual de ocurrir. Entre las distribuciones binomiales, esto no siempre es el caso, por lo que no es una suposición que debamos tomar siempre por defecto. Las distribuciones donde hay una posibilidad igual de que ocurra cualquier resultado se llaman distribución uniforme. A menudo se utiliza en las notaciones de modelos como una forma de comunicar la falta de conocimiento previo sobre los posibles resultados en nuestro modelo. Sin embargo, con suficiente práctica y experiencia a lo largo del tiempo, descubriremos que es raro que un practicante bayesiano que construye el modelo no tenga conocimiento previo cero del fenómeno que pretende modelar. También podemos simular esta distribución si tienes curiosidad por saber cómo se ve en un gráfico:

x_uni = np.linspace(-.25, 1.25, 100)plt.plot(x_uni, stats.uniform.pdf(x_uni, 0, 1))plt.xlabel("Previa (Proporción de Género)")plt.ylabel("Plausibilidad de la proporción")plt.title("Visualizando una distribución uniforme")

Una vez que hayamos identificado nuestras suposiciones sobre el modelo de proporción de género en los lobos costeros, vamos a formalizarlo en la notación de un modelo estadístico aplicado:

  • M ~ Binomial(N, p): Podemos interpretar nuestra primera línea como “el número de machos (M) que encontraremos está distribuido binomialmente con un tamaño de muestra de N y una probabilidad previa de p.” Esta línea representará nuestra función de verosimilitud, también conocida como la probabilidad de la evidencia P(E | H) en el Teorema de Bayes, que es la probabilidad de observar un resultado dado. En este caso, el “resultado dado” será el resultado de nuestros datos de muestra sobre los lobos costeros, sobre lo cual entraremos en más detalle en breve.
  • p ~ Uniform(0, 1): Por último, podemos interpretar nuestra segunda línea como “la probabilidad previa P(H) está distribuida uniformemente entre 0 y 1″.

Ten en cuenta que el símbolo ~ significa “se distribuye como” y comunica que el parámetro del lado izquierdo del modelo tiene una relación estocástica (lo que significa “determinado al azar” o tiene una propiedad de “distribución de probabilidad aleatoria” que se puede predecir estadísticamente PERO NO con precisión) con su distribución. En otras palabras, siempre que veamos el símbolo ~ que representa una relación estocástica, debemos esperar que el resultado de la declaración matemática a la derecha sea una distribución capturada por la variable a la izquierda (en nuestro caso, ya sea p o M). Expresar relaciones estocásticas contrasta fuertemente con las declaraciones matemáticas que estamos acostumbrados a ver con un solo número (en lugar de una distribución) que es el resultado de una ecuación que consiste en un símbolo =, como uno común que a menudo vemos para Regresiones Lineales, que es y = mx + b.

Lo último que necesitamos antes de poder traducir nuestro modelo a código es generar una muestra de observaciones para nuestra función de verosimilitud ( P(E | H) ), la cual, nuevamente, representa la probabilidad de observar el resultado dado de los eventos. Por ejemplo, si utilizamos los datos de muestra adquiridos por Muñoz-Fuentes y sus colegas (2009) en su estudio sobre los lobos costeros de la Isla de Vancouver, encontraremos que de los 28 lobos con género etiquetado, 15 eran machos. En este caso, podemos describir nuestra función de verosimilitud como la probabilidad de encontrar 15 machos de los 28 lobos que encontramos en nuestro conjunto de datos de muestra.

Con nuestras variables establecidas, ahora podemos usar PyMC para incorporar las notaciones del modelo que acabamos de formalizar directamente en el código:

# Leyendo nuestro conjunto de datos de muestra del estudio de Muñoz-Fuentes et al. (2009).wolf_samples = "https://raw.githubusercontent.com/vanislekahuna/Statistical-Rethinking-PyMC/main/Data/Vancouver-Island-wolf-samples.csv"islandwolves_df = pd.read_csv(wolf_samples)islandwolves_df# Determinando la proporción macho-hembra en el conjunto de datosprint(f"El conteo de valores para el sexo de los lobos costeros de la Isla de Vancouver es: ")print(islandwolves_df["Sexb"].value_counts()

# Generando las observaciones de muestra para nuestro modelo.machos = islandwolves_df["Sexb"].value_counts()["M"]hembras = islandwolves_df["Sexb"].value_counts()["F"]muestras_totales = machos + hembras

### Nuestro modelo PyMC ###with pm.Model() as wolf_ratio:  # Una distribución uniforme de posibles resultados para los valores previos.  priors = pm.Uniform("uniform_prior", lower=0, upper=1)   # Genera una distribución binomial a partir de nuestra distribución uniforme de valores previos y nuestros datos de muestra, donde 15 de los 28 lobos resultaron ser machos.  n_machos = pm.Binomial("sex_ratio", p=priors, observed=machos, n=muestras_totales)  # Muestras extraídas del modelo para generar nuestra distribución posterior.   trace_wolves = pm.sample(1000, tune=1000)######################

Ahora que nuestro modelo ha sido actualizado, no nos perdamos en los detalles y olvidemos el único propósito de construir el modelo en primer lugar. La razón por la que estamos construyendo el modelo es cuantificar la probabilidad de encontrarnos con un lobo macho a través de la distribución de posibilidades y, en esencia, clasificarlos por su verosimilitud. En otras palabras, nos interesa la probabilidad de observar 15 machos de los 28 lobos totales que hemos encontrado (nuestra función de verosimilitud) dada cada proporción plausible macho-hembra, como se representa en nuestra distribución previa (es decir, el 1% de todos los lobos en la naturaleza siendo machos, 2%, 3%, hasta el 100%). El resultado final debe ser una distribución de probabilidad posterior que muestre las proporciones macho-hembra más plausibles basadas en los datos de muestra que hemos observado. A partir del código anterior, la distribución posterior se captura en el objeto `trace_wolves`, que genera las probabilidades posteriores resultantes mediante el muestreo del modelo que construimos. En caso de que te lo estés preguntando, los métodos de muestreo utilizados para generar la distribución posterior es un tema que podemos cubrir en un artículo posterior, ya que PyMC proporciona varias opciones, cada una siendo un agujero de conejo profundo y expansivo por sí mismo.

El gráfico a la izquierda a continuación visualiza cómo se ve nuestra distribución posterior resultante, donde el eje x representa nuestra distribución previa de probabilidad de las proporciones de género mientras que el eje y representa la plausibilidad de esa proporción:

az.plot_trace(trace_wolves)

Y si ayuda estructurar el modelo `wolf_ratio` que acabamos de construir como variables en el Teorema de Bayes, podemos hacerlo:

Nuevamente, es importante enfatizar que el objetivo aquí es generar una distribución posterior de plausibilidades para cada valor previo en su distribución. NO ESTAMOS tratando de encontrar la proporción de género con la probabilidad más alta. Dicho esto, si estuvieras interesado en cuál sería ese valor, podríamos usar la función `find_map` de PyMC para darnos el valor previo más plausible (es decir, la proporción más plausible de machos a hembras) en la distribución posterior:

pm.find_MAP(model=wolf_ratio)

O la función `summary` de Arviz, que también nos da el máximo a priori (MAPA) de la distribución posterior, además de los intervalos de mayor densidad (a veces también llamados intervalos de compatibilidad), que son los dos valores anteriores que contienen el 89% de la distribución.

az.summary(trace_wolves, round_to=2, kind="stats", hdi_prob=0.89)

Ahí lo tienes. Acabamos de pasar por una introducción básica al lenguaje de las notaciones de modelos estadísticos utilizados en la mayoría de los estudios académicos. Luego integramos las notaciones de los modelos en nuestro primer modelo PyMC, donde simulamos la plausibilidad de encontrar un lobo costero de la isla de Vancouver en estado salvaje, dado que no sabemos nada sobre su proporción de género masculino a femenino. Familiarizarse con la notación del modelo formará una base sólida para construir modelos bayesianos más complejos utilizando PyMC. Planeamos lanzar más blogs sobre este tema en un futuro cercano, comenzando con un artículo sobre Predicción Lineal Bayesiana, ¡así que mantente atento!

Referencias:

McElreath, R. (2020). Statistical Rethinking: A Bayesian Course with examples in R and Stan. Routledge.

Muñoz‐Fuentes, V., Darimont, C. T., Wayne, R. K., Paquet, P. C., & Leonard, J. A. (2009). Ecological factors drive differentiation in wolves from British Columbia. Journal of Biogeography, 36(8), 1516–1531. https://www.raincoast.org/files/publications/papers/Munoz_et_al_2009_JBiogeog.pdf

Muñoz-Fuentes, V., Darimont, C.T., Paquet, P.C., & Leonard, J.A. (2010). The genetic legacy of extirpation and re-colonization in Vancouver Island wolves. Conservation Genetics, 11, 547–556. https://digital.csic.es/bitstream/10261/58619/1/evol.pdf

We will continue to update Zepes; if you have any questions or suggestions, please contact us!

Share:

Was this article helpful?

93 out of 132 found this helpful

Discover more

Inteligencia Artificial

Mejorando los Modelos de Lenguaje con Indicaciones Analógicas para Mejorar el Razonamiento

En los últimos años, los modelos de lenguaje han demostrado una notable habilidad para entender y generar texto simil...

Inteligencia Artificial

AI Sesgo Desafíos y Soluciones

¿De dónde proviene el sesgo en la inteligencia artificial? Una vez que lo encontramos, ¿cómo podemos reducirlo o elim...

Inteligencia Artificial

Potenciando la IA en Dispositivos Qualcomm y Meta colaboran con la tecnología Llama 2

El lanzamiento de Llama 2, la nueva versión de código abierto de Meta, ha generado discusiones sobre los casos de uso...

Inteligencia Artificial

Gigantesco Telescopio Adopta Robots de Mantenimiento Inteligentes

Cinco sistemas y plataformas de robots inteligentes han sido autorizados por las autoridades para mantener el Telesco...