Programación probabilística con NumPyro

Programación probabilística con NumPyro Probabilistic Programming with NumPyro Trabajo de Fin de Grado Curso 2023–2024 Autor Gonzalo García Fernández Director Miguel Palomino Tarjuelo Doble grado en Ingeniería Informática y Matemáticas Facultad de Informática Universidad Complutense de Madrid Programación probabilística con NumPyro Probabilistic Programming with NumPyro Trabajo de Fin de Grado en Ingeniería Informática Autor Gonzalo García Fernández Director Miguel Palomino Tarjuelo Convocatoria: Septiembre 2024 Doble grado en Ingeniería Informática y Matemáticas Facultad de Informática Universidad Complutense de Madrid Agradecimientos En primer lugar me gustaría agradecer a mi tutor y director del trabajo de fin de grado, Miguel Palomino Tarjuelo, por su dedicación y orientación a lo largo de todo el proyecto. Su experiencia, apoyo y paciencia han sido cruciales durante la redacción de todo el proyecto. También me gustaría agradecer a mi familia y amigos, en especial a Andrea Chivite y Sergio Javier Hernando, que me han acompañado y ayudado durante todo el proceso. v Resumen Programación probabilística con NumPyro La programación probabilística es un enfoque innovador para el desarrollo de modelos probabilísticos complejos y la ejecución eficiente de inferencias. Este trabajo presenta una introducción a la programación probabilística con NumPyro, una biblioteca de Python que aprovecha las capacidades de cálculo automático y el paralelismo de JAX para realizar inferencias bayesianas avanzadas. A lo largo de este trabajo, se proporcionará una explicación detallada de las funcionalidades más importantes de NumPyro, acompañada de ejemplos prácticos y explicaciones claras para facilitar la comprensión de la biblioteca. Se irán introduciendo ejemplos progresivamente para mejorar el entendimiento del lector. Por último, se presentará una posible aplicación de la programación probabilística y se expondrán las conclusiones derivadas del trabajo. Palabras clave NumPyro, programación probabilística, manejadores de efectos, JAX, Bayes. vii Abstract Probabilistic Programming with NumPyro Probabilistic programming is a new approach to the development of complex probabilistic models and the efficient execution of inferences. This paper presents an introduction to probabilistic programming with NumPyro, a Python library that uses automatic computation capabilities and parallelism of JAX to perform efficient Bayesian inference. Throughout this paper, there will be a detailed explanation of the most important functionalities of NumPyro. It will also be complemented by practical examples and clear explanations to facilitate the understanding of the language. Finally, an application of probabilistic programming will be presented and the conclusions derived from the work will be presented. Keywords NumPyro, Probabilistic programming, Effect Handlers, JAX, Bayes. ix Cap´ ıtulo 1 Introducción 1.1. Motivación Hoy en día, todo está basado en los datos y en la tecnología. Por lo tanto, entender y predecir problemas complejos es más importante que nunca. Muchos de los problemas a los que nos enfrentamos como predecir el clima, encontrar patrones en datos financieros o detectar enfermedades tienen algo en común: la incertidumbre. Nos encontramos en un entorno donde los datos a menudo son incompletos, ruidosos y cambian constantemente. Entonces, ¿cómo podemos hacer inferencias precisas y tomar decisiones acertadas en medio de tanta incertidumbre? Ahí es donde entra en juego la programación probabilística. La programación probabilística es una potente herramienta que mezcla la teoría de la probabilidad con la programación. Con esta combinación, podemos crear modelos que reflejen la incertidumbre presente en los datos y en las diferentes acciones que se llevan a cabo con esos datos. En lugar de trabajar solo con datos fijos y deterministas, la programación probabilística permite que los modelos consideren variables que pueden tomar diferentes valores, asignando a cada uno una probabilidad. Esta capacidad de manejar la incertidumbre y hacer inferencias a partir de datos incompletos o inciertos hace que la programación probabilística sea una herramienta de gran utilidad. 1.2. Objetivos Este trabajo es una introducción a la programación probabilística y, específicamente, a un lenguaje de programación llamado NumPyro. El objetivo principal de este trabajo es hacer una introducción clara y concisa a la programación probabilística. En concreto, me gustaría mostrar el funcionamiento de la programación probabilística con NumPyro: qué es, sus componentes y cómo funciona internamente de tal forma que si alguien sin ningún conocimiento previo de programación probabilística leyera el trabajo, obtendría un conocimiento claro sobre cómo utilizarlo. Además, me gustaría incluir alguna de las aplicaciones para entender cómo se puede aplicar la programación probabilística en nuestro día a día. 1 2Capítulo 1. Introducción Este trabajo complementa y es complementado por mi trabajo de fin de grado en Matemáticas (Fernández (2024)) que detalla exhaustivamente los métodos de inferencia desde un punto de vista mucho más matemático. 1.3. Plan de trabajo y estructura El plan de trabajo seguido durante el trabajo de fin de grado es el siguiente: Para el segundo capítulo, se ha estudiado la situación de la programación probabilística hoy en día, un marco más general de la programación probabilística de por qué surgió y por qué es necesario. Para el tercer capítulo, que hace referencia al lenguaje probabilístico de NumPyro, se ha estudiado principalmente la documentación de NumPyro. Además, la he complementado con otros artículos hacen referencia al lenguaje o a diferentes aspectos del mismo. En el capítulo cuatro, se realiza el estudio de la red neuronal bayesiana. Este se ha realizado investigando sobre redes neuronales bayesianas para aprender su funcionamiento, investigando sobre diferentes aspectos económicos que pueden ser relevantes en la predicción de acciones y por último investigando cómo se programa una red neuronal bayesiana utilizando NumPyro. Por último, haremos referencia a los capítulos anteriores explicando las conclusiones que hemos sacado del trabajo y analizando los resultados de la red neuronal bayesiana. Cap´ ıtulo 2 Programación probabilística 2.1. Introducción a la programación probabilística La creación de modelos es una práctica que comienza cuando somos niños y se prolonga a lo largo de nuestra vida. Los niños construyen modelos de aviones y luego los queman con petardos solo para ver qué pasa. Los científicos utilizan ratones como un modelo de organismo para simular cómo se comportaría un experimento similar en humanos. Un modelo probabilístico es una representación matemática que utiliza conceptos de probabilidad para describir y predecir eventos inciertos. Estos modelos se basan en la idea de que no podemos predecir con seguridad el resultado de un evento, pero podemos asignar probabilidades a los posibles resultados para calcular la probabilidad de que ocurran. Es decir, un modelo probabilístico nos permite calcular la probabilidad de diferentes resultados en situaciones donde existe incertidumbre. Entonces, en nuestro día a día muchas veces estamos considerando modelos probabilísticos al pensar en una situación y en lo que podría suceder. La programación probabilística se inventó para simplificar y establecer un entorno, para poder realizar modelos probabilísticos y, principalmente, ejecutarlos. Antes de la aparición de la programación probabilística, crear modelos probabilísticos complejos requería de tener un conocimiento amplio en programación y matemáticas, lo que hacía muy difícil realizar esta práctica debido a que cada nuevo modelo necesitaba un código prácticamente nuevo. Además era muy complicado representarlos debido a qué la forma habitual es a través de un grafo. El problema con la representación a través de un grafo, es que cuando aumento la complejidad del modelo, aumentan las dependencias y el número de variables y el grafo se vuelve inmanejable. La programación probabilística utiliza la estadística y la informática para crear programas que incluyen variables aleatorias y muestran relaciones probabilísticas entre estas variables en un entorno que permite representarlo como si fuera un programa informático y ejecutarlo a través de métodos de inferencia. A diferencia de la programación tradicional, donde se espera que un programa produzca una salida determinada a partir de una entrada dada, en la programación probabilística los programas generan distribuciones de probabilidad como salidas. Estas distribucio3 4Capítulo 2. Programación probabilística nes de probabilidades extraídas se llaman distribuciones posteriores y nos permiten además de conocer el valor más probable para la variable, la incertidumbre, es decir, cómo de seguros estamos sobre la salida que hemos obtenido. Estas distribuciones posteriores se obtienen a partir de métodos de inferencia, en concreto, métodos de inferencia bayesiana que nos permiten obtener la distribución de probabilidad de una variable. Además, todos los lenguajes probabilísticos tienen en común que añaden dos nuevas funcionalidades principales a la programación tradicional. La primera sample, que nos permite extraer muestras de una función de probabilidad. Esto nos será de gran utilidad, porque muchas variables pasarán de ser variables deterministas a variables aleatorias que tendrán asignada una distribución previa basada en nuestras creencias. La otra operación será obs, que nos permite condicionar por una variable observada, es decir, nos permite ajustar la salida de un programa a unos datos que han sido observados con anterioridad. Estas dos funcionalidades y los métodos de inferencia que mencionábamos en el párrafo anterior son los pilares de la programación probabilística. En el capítulo 3, donde trataremos en detalle un lenguaje probabilístico, extenderemos todos los conceptos de la programación probabilística, excepto los métodos de inferencia que los veremos brevemente debido a que están explicados en mucho mayor detalle en mi trabajo de fin de grado de Matemáticas. 2.2. Evolución de la programación probabilística El primer lenguaje de programación probabilística fue BUGS(Bayesian inference using Gibbs Sampling) inventado en la Universidad de Cambridge alrededor de 1989, unos años más tarde de la introducción del algoritmo de muestreo de Gibbs inventado por los hermanos Stuart y Donald Geman alrededor de 80 años después de la muerte de Gibbs. BUGS es un programa que lleva a cabo inferencia bayesiana en problemas estadísticos usando el algoritmo de muestreo de Gibbs. BUGS asume un modelo de probabilidad, donde todas las variables involucradas se tratan como variables aleatorias. El objetivo principal es calcular la distribución posterior de los parámetros desconocidos y los datos no observados, dados los datos observados. Esto se logra mediante el método Montecarlo, específicamente a través del muestreo de Gibbs. El programa cuenta con un conjunto básico de comandos que permiten controlar una sesión en la que se analiza un modelo estadístico expresado en el lenguaje de BUGS. Un compilador procesa el modelo y los datos disponibles en una estructura interna adecuada para la computación eficiente, y un muestreador usa esta estructura para generar valores apropiados de las cantidades desconocidas. Adicionalmente, se proporciona un conjunto de funciones en para el análisis y la visualización de los resultados, así como para asegurar la convergencia del muestreo(Spiegelhalter et al. (1996)). BUGS fue un herramienta muy novedosa debido a que fue la primera que permitió establecer un entorno para poder realizar inferencia bayesiana. Sin embargo, al ser el primero de estos lenguajes probabilísticos nos encontramos que no tiene algunas características que hoy en día tienen otros lenguajes probabilísticos. Por ejemplo, el muestreo de Gibbs es un método sencillo y robusto pero puede ser ineficiente si 2.2. Evolución de la programación probabilística 5 hay modelos de altas dimensiones. Otra desventaja es que BUGS proporciona un lenguaje de modelado gráfico declarativo para describir modelos probabilísticos y aunque permite realizar un gran número de modelos, la flexibilidad es limitada. Por último, al ser un lenguaje más antiguo, puede llegar a ser un problema para algunos usuarios. El primero de los lenguajes más novedosos es Stan, que pretende establecer un lenguaje de programación probabilística para el modelado para usando con métodos de inferencia poder realizar un modelo robusto, escalable y eficiente. Stan se diferencia del lenguaje BUGS anteriormente descrito en que está basado en un nuevo lenguaje imperativo de programación probabilística que es más flexible y expresivo que los lenguajes de modelado gráfico declarativo en los que se basa BUGS, en aspectos como la declaración de variables con tipos y el soporte para variables locales y declaraciones condicionales. Además, también se diferencia en que pasamos de usar el método de muestreo de Gibbs a utilizar el Montecarlo Hamiltoniano (HMC), un muestreador más eficiente y robusto que el mencionado anteriormente(Carpenter et al. (2017)). A continuación, podemos observar de un código de ejemplo de Stan. data { int<lower=0> N; // N >= 0 int<lower=0,upper=1> y[N]; // y[n] in { 0, 1 } } parameters { real<lower=0,upper=1> theta; // theta in [0, 1] } model { theta ~ beta(1,1); // prior y ~ bernoulli(theta); // likelihood } Por último, me gustaría recalcar de Stan que todas las expresiones son de tipo estático, incluidas las variables. Esto significa que su tipo se declara en tiempo de compilación como parte del modelo y no cambia durante la ejecución. Este comportamiento es el mismo que se encuentra en lenguajes de programación como C++ o Java pero es diferente al comportamiento de los lenguajes Python o R. Stan ha sido un gran referente en la programación probabilística, sin embargo, actualmente han surgido nuevos lenguajes y bibliotecas de programación probabilística que han comenzado a ganar popularidad como PyMC3, Pyro, o NumPyro. Estas nuevas herramientas aprovechan las capacidades de las nuevas librerías de machine learning y computación en hardware acelerado como JAX. Estas ofrecen enfoques nuevos y flexibles para el modelado bayesiano. En nuestro caso, hemos escogido NumPyro para este trabajo y en el capítulo 3 desglosaremos todo lo relativo a este lenguaje de programación probabilística. 6Capítulo 2. Programación probabilística 2.3. Posibles aplicaciones de la programación probabilística Una vez, hemos visto una pequeña introducción a la programación probabilística y a su desarrollo a lo largo de la historia me gustaría introducir algunas de las aplicaciones más relevantes que puede ofrecer. 1. Simulación restringida: La generación de gráficos procedurales es una aplicación de la programación probabilística que resulta muy visual y fácil de entender. Por ejemplo, imaginése que en una película necesitan crear un bosque con gráficos por ordenador. Para eso, contratan a un programador de gráficos procedurales que hace un modelo generativo, que básicamente es un simulador aleatorio que crea un árbol diferente cada vez que se ejecuta. Ahora, supongamos que en ese bosque cada árbol tiene que cumplir con ciertas reglas o restricciones. Las decisiones aleatorias que toma el modelo generativo, como cuánto se alargan las ramas, cuántas ramas salen del tronco, los ángulos de las ramas nuevas, cuándo termina una rama, etc., determinan cómo se ve cada árbol. Cada árbol que se genera es el resultado de una configuración específica de esas variables aleatorias en el programa. Cuando se aplican las restricciones, la distribución original de los árboles generados se ajusta para que todos los árboles cumplan con esas reglas, es decir, al condicionar por las restricciones se transforma la distribución previa en una distribución posterior en la que todos los árboles posteriores cumplen con las restricciones. 2. Predicción de acciones: La programación probabilística se utiliza para predecir el precio de las acciones porque permite manejar de manera clara la variabilidad y complejidad que existe en los mercados financieros. Los precios de las acciones pueden cambiar mucho y están influenciados por muchos factores impredecibles, como cambios en la economía, eventos políticos o cómo actúan los inversores. Con la programación probabilística se pueden crear modelos que capturen esta incertidumbre y la reflejen en las predicciones. Por ejemplo, la programación probabilística nos da una distribución de posibles resultados en lugar de una sola predicción fija. Esto es muy importante para entender los riesgos que se están tomando al comprar o vender una acción. Además, la programación probabilística permite que estos modelos se actualicen continuamente a medida que se obtienen nuevos datos del mercado, ajustando la distribución previa para mantener las predicciones actualizadas con las condiciones actuales de los mercados y los consumidores. 3. Reconocimiento de patrones: En robótica, las imágenes capturadas por cámaras o señales de sensores suelen ser a menudo datos ruidosos o incompletos que presentan un gran problema para que los robots interpreten lo que ocurre a su alrededor de manera precisa. Los modelos probabilísticos permiten a los robots gestionar esta incertidumbre al considerar múltiples opciones sobre lo que pueden estar detectando y al tomar decisiones más informadas en función de las circunstancias en las que se encuentre. Por ejemplo, en el reconocimiento de objetos, un robot puede identificar correctamente un objeto, incluso en 2.4. Explicación de un programa probabilístico 7 condiciones de poca luz o desde ángulos complicados, gracias a la capacidad de los modelos probabilísticos para medir la incertidumbre de las posibles interpretaciones. Un ejemplo de esto pueden ser los sistemas autónomos (Shamsi et al. (2020)). 4. Variabilidad de Diagnósticos: En el campo de la medicina, la diferencia entre pacientes y la incertidumbre en los diagnósticos son comunes. La programación probabilística permite modelar esta incertidumbre al construir modelos que permitan saber la probabilidad de diferentes diagnósticos y obtener resultados basados en los síntomas y datos del paciente. Esto permitiría a los médicos evaluar diferentes opciones teniendo una forma de saber cómo de seguros están de su decisión. En capítulos posteriores, se mostrará un ejemplo utilizando una red neuronal bayesiana de cómo se utiliza la programación probabilística para la predicción de acciones y cómo su interpretabilidad ayuda en la toma de decisiones. 2.4. Explicación de un programa probabilístico A continuación, vamos a explicar un programa probabilístico muy básico extraído del artículo van de Meent et al. (2018), que se ha utilizado también a lo largo de esta introducción. El modelo es es el siguiente: (let [prior (beta a b) x (sample prior) likelihood (bernoulli x) y (sample likelihood)] (observe likelihood y) x) Este modelo podría modelar un experimento para comprobar si una moneda está trucada o no, es decir, si hay mayor probabilidad de que salga cara o cruz o es neutral como debería ser. Primero se establece la variable prior que representa una distribución de probabilidad, en concreto, beta(a,b)donde aybson parámetros. A continuación, xse convierte en una muestra de la variable prior. La extracción de la muestra la hacemos a través de una de las nuevas primitivas de las que hemos hablado a lo largo de esta introducción, la función sample. Además la variable likelihood representa una distribución de probabilidad, en concreto, bernoulli(x)eyrepresenta una muestra de la variable likelihood. Por último, lo que hace este programa a través de la primitiva observe es condicionar y actualizar la distribución previa de la variable likelihood, lo cual actualizará la distribución de la variable prior. Es decir, si y, en lugar de ser una única variable observada fuera un conjunto de variables observadas de tamaño nsuficientemente grande y repitiéramos este proceso un número de veces n, la información que se obtiene es una distribución posterior para la variable x, que nos indica cuál es su valor más probable. Si es 0.5 la moneda no está trucada y si es diferente de 0.5 la moneda está trucada. Además, nos permite saber con qué certeza lo sabemos. Cap´ ıtulo 3 NumPyro En este capítulo explicaremos en detalle el lenguaje de programación probabilística escogido. Veremos qué es NumPyro, sus características principales y fundamentos y cómo utilizarlo con algunos modelos de programación probabilística. Para ello, utilizaremos distintas fuentes a lo largo de los dos siguientes capítulos, principalmente NumPyro (2017), Pyro (2017). 3.1. ¿Qué es NumPyro? NumPyro es un entorno de programación probabilística que permite a los usuarios construir modelos estadísticos complejos y realizar inferencias sobre ellos de manera eficiente. NumPyro es una biblioteca de Python que surge como una versión más rápida de Pyro, un entorno de modelado probabilístico desarrollado por Uber AI Labs. NumPyro se basa en Pyro en cuanto a que utilizan la misma interfaz, las mismas primitivas y las mismas abstracciones con manejadores de efectos. La gran diferencia que encontramos es que mientras Pyro esta contruido sobre Pytorch, NumPyro está construido sobre JAX, una biblioteca desarrollada por Google que permite la diferenciación automática y la compilación just-in-time (JIT) lo que resulta en mejoras significativas en el rendimiento. La utilización de JAX por NumPyro será detallada durante este capítulo, pero la relación entre NumPyro y JAX es crucial para entender por qué JAX permite a NumPyro aprovechar la aceleración en hardware, como GPUs, y realizar cálculos matemáticos de manera muy rápida. Esto significa que los modelos probabilísticos pueden ser entrenados y evaluados mucho más rápidamente que con otras herramientas, lo que es particularmente útil en aplicaciones que requieren un alto rendimiento. Además de lo mencionado anteriormente, es un lenguaje más simple que otros como Stan, lo que permite contsruir modelos complejos de manera intuitiva y con menos líneas de código. Esto hace que NumPyro sea más accesible para aquellos que necesitan personalizar los modelos sin lidiar con complejidad adicional. En las siguientes secciones iremos desglosando los diferentes aspectos de NumPyro. A continuación, podemos ver su logotipo en la figura 3. 9 16 Capítulo 3. NumPyro 3.2.3. Manejadores de efectos En esta subsección veremos cómo se utilizan los manejadores de efectos como una abstracción para implementar lenguajes de programación probabilística. Veremos primero en qué consisten y cómo se construyen basándonos principalmente en Pretnar (2015), Plotkin y Pretnar (2009), Phan et al. (2019), Kammar et al. (2013) y Davidson-Pilon (2015). A continuación, investigaremos cómo se usan los manejadores de efectos en NumPyro y describiremos algunos ejemplos de manejadores. Las operaciones con efectos son operaciones que interactúan con alguna parte del código que debe manejarse para ejecutarse, como get yset para el almacenamiento, read yprint para la entrada y salida, o raise para las excepciones. De esta forma surgen los manejadores no solo de los más habituales que son los de excepciones, sino de cualquier otro efecto, dando lugar a un concepto novedoso que, entre otras cosas, puede capturar a dónde se envía un programa, el retroceso, múltiples hilos de ejecución o guardar el estado de un programa para retomarlo posteriormente. Para poder explicar la teoría de los manejadores de efectos nos apoyaremos en el artículo Pretnar (2015), que introduce un lenguaje muy sencillo que abarca un lenguaje de programación general. En este lenguaje tenemos valores fijos y computaciones que pueden verse modificadas por efectos computacionales. El lenguaje es el siguiente: Figura 4: Lenguaje para los manejadores de efectos Este lenguaje se compone de valores que son las partes del código que no cambian durante la evaluación. Los valores pueden ser una variable (x), una constante booleana, una función determinista que tome una entrada xy produzca una salida co un manejador h. El siguiente elemento del lenguaje son los manejadores, que definen cómo se deben procesar ciertas operaciones dentro de una computación. Un manejador tiene varias clausulas, primero encontramos la cláusula return x→cr que delimita cuándo acaba una operación y devuelve un valor x(cr es la computación resultante). El otro tipo de cláusula son las cláusulas de operación donde cada operación op tiene una computación casociada. La operación toma una entrada xy una continuación k, que es lo que debería pasar después de que la operación termina y 3.2. Componentes principales de NumPyro 17 produce una nueva computación c. Es decir, una vez aparece la operación, se realiza la operación y se continua por la computación k. Por último, los distintos tipos de computaciones: return v: Devuelve el valor v. op(v;y.c): Realiza una llamada a una operación op, donde se le pasa como entrada el valor v, se ejecuta la operación y se almacena el valor en y. A continuación, el programa sigue por c, que es la parte del programa que continúa. Un ejemplo sería si queremos leer de memoria, vsería la dirección de memoria, yel valor que hay en la memoria y cpodría ser lo que hacemos luego con el valor que hemos traído de memoria. do x←c1in c2: Secuencia de acciones donde c1se evalúa primero, y su resultado se asigna a xantes de proceder a c2. if vthen c1else c2: Un condicional que evalúa c1si ves true yc2si v es false. v1v2: Aplicación de una función donde v1es la función y v2es el argumento. with vhandle c: Maneja un cálculo cutilizando el manejador v. Por último, cabe destacar que para la semántica, tener la continuación a la operación que se va a realizar es de gran utilidad pero para un programador no lo es porque solo quiere recuperar el valor de la operación. Entonces a la hora de ejecutar podemos definir un efector genérico op_gen que lo que hace es pasar el parámetro de la función y devolver el valor de y, es decir: op_gen ≜fun x7→ op(x;y. return y) Además nos permite recuperar la operación original haciendo, do y←op vin c. A continuación, veremos una semántica operacional de paso pequeño de cómo proceden las computaciones cuando hay manejadores de efectos. Esta semántica se basa en la idea de que las llamadas a computaciones no producen efectos reales sino que actúan como señales que se propagan hasta encontrar un manejador con su cláusula correspondiente. También, una computación que no sea capturada por un manejador se considerará como finalizada, es decir, no pude reducirse más. La semántica operacional sacada de Pretnar (2015) se define de la siguiente forma: 1. c1⇝c′ 1 do x←c1in c2⇝do x←c′ 1in c2 2. do x←return vin c⇝c[v/x] 3. do x←op(v;y.c1)in c2⇝op(v;y.do x←c1in c2) 18 Capítulo 3. NumPyro 4. if true then c1else c2⇝c1 5. if false then c1else c2⇝c2 6. (fun x7→ c)v⇝c[v/x] En las siguientes reglas, definimos h=handler{return x7→ cr, op1(x;k)7→ c1, . . . , opn(x;k)7→ cn}: 7. c⇝c′ with hhandle c⇝with hhandle c′ 8. with hhandle (return v)⇝cr[v/x] 9. with hhandle opi(v;y.c)⇝ci[v/x, (fun y 7→ with hhandle c)/k](1 ≤i≤n) 10. with hhandle op(v;y.c)⇝op(v;y.with hhandle c)(op /∈ {op1, . . . , opn}) Explicaremos a continuación todas las reglas para terminar de aclarar cómo funcionan los manejadores de efectos. En las primeras seis reglas formalizamos el comportamiento de las computaciones. Regla 1: Afirma que si la computación c1se reduce a c′ 1entonces la secuencia do x←c1in c2se reduce a do x←c′ 1in c2. Esto implica que la computación c2espera a que c1termine de evaluarse y si c1se reduce a otro estado el mismo cambio se produce en la secuencia. Regla 2: Si la secuencia c1es simplemente return vin centonces lo que devuelve la secuencia es ccambiando el valor de vpor x. Regla 3: Esta regla es relativa al llamado de operaciones. Lo que dice es que si una operación op(v;y.c1) está en la secuencia de do x←c1in c2se pueden intercambiar, dando lugar, a que se puede realizar la misma operación poniendo como continuación: do x←c1in c2. Regla 4 y regla 5: Estas reglas hacen referencia a la formalización del condicional. Indican que si vevalúa a True, se realiza c1y se ignora c2y si evalúa aFalse entonces se hace justo al revés. 3.2. Componentes principales de NumPyro 19 Regla 6: Esta regla formaliza el comportamiento de las funciones, que indica que aplicar fun x→ca un valor vresulta en realizar la computación c cambiando xpor v. Una vez formalizado adecuadamente el comportamiento de las distintas computaciones veamos qué sucede al aplicar un manejador de efectos h. Regla 7: Si una computación cse reduce a una computación c′, entonces la expresión with hhandle c, también se reduce a with hhandle c′. Es decir, el manejador se aplicará al final de la computación, propagando la aplicación del manejador. Regla 8: En el caso de que la computación manejada sea únicamente return v, el manejador hutiliza su cláusula específica de retorno crsustituyendo xpor v. Es decir, se ejecuta crpero modificando el valor de x. Regla 9: Esta regla se aplica cuando opiestá siendo manejada por h. Aquí, cies la cláusula de operación correspondiente en el manejador h. La expresión cise evalúa con vsustituido por xy la continuación kse reemplaza por fun y→with hhandle c. Esto implica que se vuelve a desplazar el manejador ejecutando la operación y manejando el resultado. Regla 10: Esta regla indica que si una operación op no tiene una cláusula en el manejador hentonces pasa simplemente del manejador realizando la operación y manejando la continuación. Tras haber comprendido cómo funcionan los manejadores de efectos internamente a través de su semántica operacional observaremos unos ejemplos sencillos de qué podría ser un manejador de efectos para luego verlo en NumPyro. Un ejemplo muy simple sería a través de la operación de lectura. handler{read(_, k)→ k"Bob"})es un manejador que hace que cada vez que se llame a la operación de lectura se devuelva "Bob" como cadena de caracteres constante y la continuación k se realizará con "Bob". Otro posible ejemplo es un manejador que recoge todos los procesos de escritura y las devuelve como una única cadena de caracteres. El código sería el siguiente: collect ≜handler       return x7→ return (x, ””) print(s;k)7→ do (x, acc)←k() in return (x, join sacc)        Si la cadena únicamente devuelve un valor xentonces el manejador debe devolver xy una cadena vacía. Si una computación tiene una instrucción de escritura con una cadena s, lo que hacemos es ejecutar la continuación k, ya que el manejador se irá desplazando hasta llegar al caso base donde se devolverá el valor en xy en acc la cadena vacía y luego iremos uniendo hasta llegar a una única cadena de caracteres con todas las salidas. 20 Capítulo 3. NumPyro En NumPyro los manejadores de efectos proporcionan una manera de introducir efectos computacionales en las primitivas de un programa probabilístico como sample oparam. Un ejemplo sería registrar las elecciones aleatorias realizadas en una traza de ejecución. Internamente, los algoritmos de inferencia pueden usar los manejadores de efectos para inspeccionar y modificar el comportamiento del programa. Al igual que hicimos con las primitivas, exploraremos como se utilizan algunos de los manejadores de efectos más importantes de NumPyro. 3.2.3.1. Seed El manejador de efectos seed en NumPyro se utiliza para controlar la aleatoriedad dentro de un programa probabilístico. Esto asegura que las ejecuciones del programa sean reproducibles al controlar la aleatoriedad dentro del modelo probabilístico. Cada llamada a numpyro.sample dentro de la función utiliza una nueva semilla derivada de la semilla inicial pero siguiendo un mismo patrón lo que permite la reproducibilidad. rng_seed: Es un entero, un jnp.ndarray scalar o un jax.random.PRNGKey. Especifica la semilla aleatoria que se usará para inicializar el generador de números aleatorios. fn: La función o modelo probabilístico que se está ejecutando y para el cual se desea controlar la aleatoriedad. hide_types: Una lista opcional de zonas para las que modificar la semilla. Cabe destacar que, a diferencia del lenguaje de programación probabilística Pyro, siempre es necesario envolver la primitiva sample por un manejador de seed. Un ejemplo sería el siguiente. def model(): a=numpyro.sample('a', dist.Normal(0,1)) b=numpyro.sample('b', dist.Bernoulli(0.5)) return numpyro.sample('x', dist.Gamma(2.0,1.0)) seeded_model =seed(model, rng_seed=1) x=seeded_model() 3.2.3.2. Trace El manejador de efectos trace permite registrar la ejecución de un programa probabilístico, capturando detalles sobre las muestras de variables aleatorias y sus distribuciones. Esto es útil para inspeccionar el comportamiento de un modelo y depurar. Los parámetros son los siguientes: fn: La función o modelo cuya ejecución se desea rastrear. 3.2. Componentes principales de NumPyro 21 Además el método get_trace devuelve el rastro completo de la ejecución, incluyendo todas las muestras y distribuciones correspondientes. En el resultado de este método podemos ir mirando variable por variable características como su valor, si es observable o no, la función de distribución o los argumentos. Un ejemplo de su aplicación es el siguiente: def model(): a=numpyro.sample('a', dist.Normal(0,1)) b=numpyro.sample('b', dist.Bernoulli(0.5)) return numpyro.sample('x', dist.Gamma(2.0,1.0)) seeded_model =seed(model, rng_seed=1) exec_trace =trace(seeded_model).get_trace() print(exec_trace['a']) Esta es la traza para la variable ’a’: { 'type':'sample', 'name':'a', 'fn': <numpyro.distributions.continuous.Normal object at 0x000001B23D2DA090>, 'args': (), 'kwargs': { 'rng_key': Array([3819641963, 2025898573], dtype=uint32), 'sample_shape': () }, 'value': Array(-1.1470195, dtype=float32), 'scale': None, 'is_observed': False, 'intermediates': [], 'cond_indep_stack': [], 'infer': {} } 3.2.3.3. Substitute El manejador de efectos substitute reemplaza las variables aleatorias o parámetros en un modelo con valores fijos, lo que es útil para pruebas. Los parámetros son los siguientes: substitute_data: Un diccionario que especifica las variables a sustituir y sus valores o funciones correspondientes. fn: La función o modelo en el que se realizarán las sustituciones. 22 Capítulo 3. NumPyro Un ejemplo de su utilidad se puede ver en este código. def model(): a=numpyro.sample('a', dist.Normal(0,1)) b=numpyro.sample('b', dist.Bernoulli(0.5)) return numpyro.sample('x', dist.Gamma(2.0,1.0)) seeded_model =seed(model, rng_seed=1) substituted_model =numpyro.handlers.substitute(seeded_model, {'a':0.06}) exec_trace =trace(substituted_model).get_trace() print(exec_trace['a']['value']) En este caso, podemos observar que la salida es un valor de 0.06. 3.2.3.4. Condition El manejador de efectos condition se utiliza para fijar ciertos valores observados en el modelo, condicionando variables aleatorias a datos. Esto es útil cuando se quiere restringir el valor de una variable aleatoria a valores específicos o cuando el proceso de inferencia quiere realizar pruebas para buscar un valor óptimo. Es muy parecido a substitute pero solo afecta a sample y cambia la propiedad de que la variable está observada a True. Los parámetros son los siguientes: data: Un diccionario que asigna nombres de variables a valores observados. Este manejador fija las muestras de las variables aleatorias a los valores proporcionados. Por esta razón, era de gran importancia el nombre que se le asigna internamente a cada variable aleatoria. fn: La función probabilística a la que se le aplicará el manejador. Un ejemplo de sus utilidad se puede ver en este código. def model(): a=numpyro.sample('a', dist.Normal(0,1)) b=numpyro.sample('b', dist.Bernoulli(0.5)) return numpyro.sample('x', dist.Gamma(2.0,1.0)) seeded_model =seed(seeded_model, rng_seed=1) conditioned_model =condition(seeded_model, data={'a':-1}) exec_trace =trace(conditioned_model).get_trace() print(exec_trace['a']['value']) print(exec_trace['a']['is_observed']) Por último, la salida de esta programa es el valor de aque está condicionado a -1 y True porque el manejador condition, modifica el parámetro is_observed. 3.2. Componentes principales de NumPyro 23 Por último, de esta sección cabe recalcar que hay varias formas de utilizar los manejadores de efectos. Se pueden utilizar de las siguiente manera: Los decoradores en NumPyro se utilizan para envolver la definición del modelo en NumPyro. @seed(rng_seed =random.PRNGKey(0)) def model(): a=numpyro.sample('a', dist.Normal(0,1)) b=numpyro.sample('b', dist.Bernoulli(0.5)) return numpyro.sample('x', dist.Gamma(2.0,1.0)) seeded_model =model() Otra forma de utilizar los manejadores es como funciones de orden superior, que es cpmo veíamos en los ejemplos descritos arriba. La última forma en la que principalmente podemos encontrar los manejadores son los gestores de contexto. En este caso, se aplica a bloques de código concreto dentro de la función, lo que nos facilita un control preciso sobre el alcance y la duración del manejador. def model(): with seed(rng_seed =random.PRNGKey(0)): a=numpyro.sample('a', dist.Normal(0,1)) with seed(rng_seed =random.PRNGKey(1)): b=numpyro.sample('b', dist.Bernoulli(0.5)) return numpyro.sample('x', dist.Gamma(2.0,1.0)) Como podemos ver en el ejemplo anterior, la diferencia entre los tres casos es: los decoradores aportan una mayor limpieza pero ofrecen menor versatilidad en el código; las funciones de orden superior nos permiten aplicarlo fuera de la definición del modelo, es decir, en la propia ejecución de la función y los gestores de contexto nos permiten tener diferentes manejadores dentro de la definición del modelo aunque no nos permiten modificar la ejecución de la función. 3.2.4. Primer modelo probabilístico en NumPyro En el capítulo 2, veíamos un primer modelo que mostraba si una moneda estaba trucada o no. A continuación veremos el código implementado en NumPyro. def coin_model(data=None): p=numpyro.sample('p', dist.Beta(1,1)) with numpyro.plate('data',len(data) if data is not None else 0): numpyro.sample('obs', dist.Bernoulli(p), obs=data) data =jnp.array([1,0,1,1,0,1,0,1,1,0]) seeded_model =numpyro.handlers.seed(coin_model, rng_seed = 2)(data) 24 Capítulo 3. NumPyro with numpyro.handlers.trace() as exec_trace: seeded_model =numpyro.handlers.seed(coin_model, rng_seed = 2)(data) print(exec_trace['p']) print(exec_trace['obs']) En el modelo vemos primero la definición de la variable ’p’ con una distribución previa Beta(1,1) y luego tenemos un repetidor con la longitud de data, que en este caso es 10. Por tanto, la variable’obs’ como es una variable observable, será 10 lanzamientos de moneda independientes que se han extraído de una bernoulli con parámetro p. Estas variables se utilizarán para actualizar la distribución de ’p’. Sin embargo, si ejecutamos el modelo de la manera que vimos anteriormente con los manejadores de efectos seed ytrace, lo único obtenemos del modelo es que ’p’ tiene un valor aleatorio extraído de la distribución previa y ’obs’ ha recibido el valor de los datos que habíamos observado, sin realizar ningún proceso de actualización de distribuciones lo cuál es algo primordial en la programación probabilística. Eso se debe a que no hemos realizado un proceso de inferencia. Por tanto, pasaremos a explicar cómo se produce la actualización de las distribuciones a través de los métodos de inferencia y cuáles son estos métodos. 3.2.5. Métodos de inferencia En esta subsección introduciremos en qué consisten los métodos de inferencia y para qué se utilizan. Esta sección será breve porque solo queremos introducir la importancia de los métodos de inferencia en la programación probabilística. En caso de querer entrar en detalle sobre los aspectos matemáticos de los métodos de inferencia se podrá consultar mi TFG de Matemáticas. La inferencia estadística consiste en el uso de los datos observables para inferir propiedades o características de la distribución inicial de la que hemos extraído los datos. En nuestro caso estudiaremos un tipo de inferencia estadística particular, la inferencia bayesiana. Supongamos que tenemos una muestra xde una variable aleatoria X, que será nuestra variable observable. Además, tenemos una serie de parámetros representados en el vector θ, que serán nuestras variables inferidas o no observables. Nuestro objetivo es encontrar la función de distribución de Xque denotaremos como F. Para ello, conocemos que la función de distribución Festá totalmente determinado por este vector de parámetros θ. Entonces, tenemos que encontrar cuáles son los valores de estos parámetros. La programación probabilística realiza algo mejor que darnos un valor para esos parámetros: nos ofrece una función de probabilidad para ellos. A nosotros nos interesa ir actualizando esas distribuciones con los datos observados x, para que los valores que puedan tomar los parámetros sean coherentes con los datos observados. Para ello, partimos de una distribución previa del vector de variables no observables, que denotaremos por ρ(θ). Esta distribución previa nos permite introducir nuestras creencias sobre las variables latentes y lo que queremos es ir actualizándola para llegar a la distribución posterior que se denota como ρ(θ|x). Para ellos utilizaremos que conocemos la forma que tiene la función de verosimilitud, denotada por ρ(x|θ)que representa la probabilidad de observar la muestra x dados unos parámetros. Si nos ponemos en el caso del ejemplo de la moneda que 3.2. Componentes principales de NumPyro 25 veíamos con anterioridad, nuestra muestra xsería los lanzamientos de moneda y nuestro parámetro unidimensional θsería lo que denotamos como p. La distribución previa de psería la Beta(1,1) y la función de verosimilitud sería la Bernouilli(p). Para poder obtener la distribución posterior se suele utilizar el Teorema de Bayes de la siguiente forma: ρ(θ|x) = ρ(θ)ρ(x|θ) ρ(x) donde ρ(x)es la función de densidad marginal de x. El gran problema que vamos a tener es el cálculo de esta función de densidad marginal debido a que muchas veces será una integral de gran complejidad. Cuando no es posible encontrar una solución analítica, los métodos de inferencia se utilizan para calcular la distribución posterior de forma numérica. En concreto, nosotros vamos a mencionar dos. 1. Inferencia variacional estocástica (SVI): La idea general de la inferencia variacional es buscar una función de densidad q(θ)dentro de una familia escogida, que esté lo más cerca posible de ρ(θ|x), es decir, queremos plantear un problema de optimización donde se minimice una cierta métrica. En este caso, la métrica utilizada es la cantidad ELBO. Por último, la resolución de este problema de optimización se hace a través del descenso de gradiente estocástico. 2. Métodos de Montecarlo basados en cadenas de Markov: Son un conjunto de métodos de inferencia que construyen una cadena de Markov para nuestra variable latente que tiene como distribución estacionaria la distribución posterior que buscamos. La cadena de Markov se construye partiendo de un estado inicial aleatorio. Después se van proponiendo estados, utilizando únicamente el estado anterior para proponerlo. A continuación, se calcula la probabilidad de aceptación. Si está por encima de un cierto valor se coge como siguiente estado. Este proceso se repite un número nde veces y lo que te queda después de quitar los primeros estados son muestras de la distribución posterior que buscas. El algoritmo más conocido para lograr esto es el algoritmo de Metropolis-Hastings. Sin embargo, han surgido unos nuevos métodos de Montecarlo basados en cadenas de Markov, que en lugar de escoger el nuevo estado moviéndose de forma aleatoria utilizan el gradiente de la distribución posterior para moverse hacia zonas donde la probabilidad de aceptación es mucho mayor. En concreto, nosotros utilizaremos una variación de un método de Montecarlo llamado NUTS(No U-Turn Sampler). En NumPyro, el proceso de inferencia se realiza de la siguiente manera: data_model =jnp.array([1,2,3,4]) rng_key =jax.random.PRNGKey(0) kernel =NUTS(model) mcmc =MCMC(kernel, num_warmup=1000, num_samples=2000, num_chains = 4) mcmc.run(rng_key, data=data_model) 32 Capítulo 3. NumPyro mcmc.print_summary() La figura anterior representa un modelo de efectos aleatorios, que se suele utilizar para eplicar la variabilidad de los datos que no puede explicarse con efectos fijos. En la parte anterior al modelo definimos n,meY_obs, que son los tamaños y las observaciones pertinentes. Hemos supuesto que los datos observados están ordenados por grupo, es decir, primero vienen las 10 observaciones de alpha0, luego las 12 de alpha1y así sucesivamente. El vector conversion lo utilizaremos posteriormente para mapear cada índice con su alphaj, gracias a las propiedades de JAX con los vectores. En el modelo, definimos primero las variables que no dependen de otras. Posteriormente, utilizamos un repetidor para poder muestrear los diferentes alphaj. De ahí sacamos, alpha_vector, que será un vector del mismo tamaño de las observaciones que te dirá que alpha se debe utilizar para cada observación. Posteriormente, condicionamos a las variables observadas con un repetidor. El repetidor condicionará a len(Y)variables pero podemos realizar el condicionamiento sin poner índices gracias a las operaciones vectoriales de JAX y simplemente se utilizará la longitud de los datos para conocer el tamaño del vector. Por último, realizamos el proceso de inferencia que internamente maneja la aleatoriedad con rng_key y utilizamos el algoritmo NUTS visto anteriormente para realizar la inferencia utilizando la diferenciación automática de JAX. Por último, observamos las distintas estadísticas que salen del modelo. Parámetro Mean Std Median 5.0 % 95.0 % n_eff r_hat alpha[0] 1.05 0.27 1.06 0.60 1.47 2566.69 1.00 alpha[1] 1.53 0.25 1.53 1.13 1.96 2606.97 1.00 alpha[2] 1.55 0.22 1.55 1.18 1.90 2847.20 1.00 alpha[3] 1.00 0.31 0.99 0.50 1.49 2293.41 1.00 alpha[4] 1.26 0.27 1.25 0.86 1.73 1960.04 1.00 mu 1.14 0.37 1.16 0.56 1.77 1672.09 1.00 omega 2.35 1.37 2.12 0.21 4.19 2131.91 1.00 tau 1.24 0.24 1.22 0.83 1.58 3199.71 1.00 Tabla 1: Resumen de estadísticas MCMC Además, n_eff es una medida que representa las muestras efectiva. Por último r_hat representa el coeficiente de Gelman-Rubin del que se puede encontrar más información en ?. Que sea 1 indica que todas las cadenas han convergido bien y están muestreando de la misma distribución posterior. Si esta coeficiente es mayor que 1.1, entonces es posibles que se necesite un mayor número de parámetros o un ajuste de distribuciones previas. Cap´ ıtulo 4 Aplicaciones En este capítulo, utilizaremos la programación probabilística para abordar la predicción del precio de cierre de las acciones de la empresa Microsoft. Utilizaremos un modelo probabilístico utilizando redes neuronales bayesianas. Comenzaremos haciendo una introducción a las redes neuronales básicas y a las redes neuronales bayesianas. Posteriormente, introduciremos el experimento que vamos a realizar con un análisis y preprocesamiento de datos. Por último, presentaremos el modelo. 4.1. Introducción a las redes neuronales Antes de introducir las redes neuronales bayesianas es importante conocer la base de las redes neuronales utilizando Goan y Fookes (2020). Para ello usaremos el perceptrón multicapa que sirve como base de las redes neuronales. La figura 8 representa un perceptrón multicapa. Figura 8: Perceptrón multicapa Para una entrada xde dimensión N1y suponiendo que tenemos una única capa 33 34 Capítulo 4. Aplicaciones oculta por simplicidad, la salida de la red fpuede modelarse como: ϕj=σ( N1 X i=1 xi;w1 ij)o de forma matricial Φ=σ(XTW1) f=g( N2 X j=1 ϕj;w2 jk)o de forma matricial F=g(ΦW2). donde σes la función de activación. El parámetro wn ij representa el peso de la salida de la neurona ien la capa n−1 sobre la neurona jde la capa n, siendo 0 la capa de las entradas. La primera ecuación modela la salida de la capa oculta, que tendrá una dimensión N2. La salida de la red es entonces una suma sobre las N2salidas de la capa oculta anterior a la que se le aplica una función gque suele ser la identidad en una regresión y la sigmoide en la clasificación binaria. El esquema anterior puede añadir muchas más capas ocultas, donde la entrada de cada capa es la salida de la capa inmediatamente anterior. Para la capa oculta, cada neurona se comporta de la siguiente manera. Primero se toman como entradas las salidas de la capa anterior y se hace una combinación lineal con la matriz de pesos W. Posteriormente, a este resultado se le aplica una función no lineal(σ) para que la red neuronal pueda predecir situaciones más complejas. A esta función no lineal se la denomina función de activación. Algunas de las funciones de activación más comunes son la función sigmoide (σ(x) = 1 1+exp(−x)), RELU (σ(x) = max(0, x)) o la tangente hiperbólica (σ(x) = exp(x)−exp(−x) exp(x)+exp(−x)). Por último, el objetivo de las redes neuronales es encontrar los pesos que mejor se ajustan a los datos minimizando una función de coste o error (E). El error actual se consigue comparando la salida del modelo con la salida esperada a través de alguna función de coste como puede ser el error cuadrático medio. Después se utiliza el valor de los pesos usando el gradiente de la función del error de la siguiente manera. wij ←wij −α∂E ∂wij Una vez vista la base de las redes neuronales veamos en qué se diferencian las redes neuronales bayesianas. En el enfoque visto anteriormente, los pesos del modelo no se tratan como variables aleatorias sino que se asume que los pesos tienen un valor verdadero que es desconocido. En las redes neuronales bayesianas queremos tratar los pesos como variables aleatorias y aprender una distribución de esos parámetros condicionado a lo que podemos observar en los datos de entrenamiento. Durante el proceso de entrenamiento de las redes neuronales bayesianas los pesos del modelo se infieren en base a los datos que podemos observar, actualizando la distribución posterior a través de métodos de inferencia. Cabe destacar que para la ejecución del modelo se coge una muestra de la distribución de ese momento para poder sacar una salida y luego poder actualizar la distribución. Es decir, el proceso de entrenamiento de la red neuronal, en lugar de realizarse como en un perceptrón multicapa mediante retropropagación, aquí utiliza la salida del modelo como variable observable para actualizar los dife- 4.2. Presentación de los datos y preprocesamiento de datos 35 rentes pesos a través del método de inferencia elegido, ya sea inferencia variacional o métodos de Montecarlo. Posteriormente, durante la fase de test, Se utilizan diferentes valores de los pesos extraídos de la distribución posterior para cada entrada x, construyendo varios modelos y proporcionando como salida final la media de los resultados obtenidos en las diferentes ejecuciones del modelo. Siguiendo lo que hemos visto durante todo el trabajo para un modelo en general, los pesos deberán tener una distribución previa ρ(ω)y tenemos que elegir una función de verosimilitud ρ(ω|x)donde xson los datos observados. Para escoger la función de verosimilitud se suele utilizar la salida de la red neuronal. Un ejemplo muy utilizado es una normal que tenga como media la salida de la red neuronal. La salida de la red neuronal será un valor que calcularemos utilizando una muestra de la distribución de cada peso que tengamos en la red neuronal en ese momento. Una vez obtenidas las distribuciones posteriores para los parámetros, se podrán realizar predicciones con la red neuronal. Estas predicciones no se realizarán utilizando un único conjunto de parámetros fijos como se haría en un perceptrón multicapa sino que se considerarán diferentes conjuntos de parámetros, es decir, ejecutaremos el modelo un número de veces mcon diferentes parámetros y mi predicción será la media de todas esas ejecuciones. Esto se puede expresar de manera matemática de la siguiente forma: Eρ[f] = Zf(ω, x)ρ(ω|x)dω donde f(ω, x)representa la salida de la red neuronal con esos parámetros y ρ(ω|x) es la distribución posterior de los parámetros. Una vez realiza la introducción a las redes neuronales, comenzaremos a realizar nuestro experimento. 4.2. Presentación de los datos y preprocesamiento de datos El dataset de Yahoo Finance para Microsoft contiene información histórica sobre los precios de las acciones y otros aspectos financieros relevantes. A continuación, se presenta una descripción detallada de cada columna en el dataset y su significado. 4.2.1. Descripción de las Columnas Fecha: La fecha en que se registraron los datos. Precio de Cierre: El precio de cierre de la acción al final de la jornada de negociación. Este es el valor más importante para los inversores ya que refleja el último precio con el que se operó ese día. Precio de Apertura: El precio de apertura de la acción al comienzo de la 36 Capítulo 4. Aplicaciones jornada. Este valor indica el primer precio al que se realizó una transacción en el mercado ese día. Máximo: El precio más alto alcanzado durante todo el día. Mínimo: El precio más bajo alcanzado durante la jornada. Volumen: La cantidad de acciones compradas/vendidas durante la jornada. Un alto volumen puede indicar un alto nivel de interés en la acción. Podemos ver un ejemplo de cómo es el dataset en la siguiente tabla. Fecha Valor apertura Maximo Minimo Cierre Volumen 374 2015-12-23 55.700001 55.880001 55.439999 55.820000 27279800 1233 2019-05-24 126.910004 127.419998 125.970001 126.239998 14123400 4.2.2. Análisis y preprocesado de Datos El análisis y preprocesamiento de datos es muy importante para comprender la integridad, limpieza y las relaciones entre los datos. A continuación detallaremos algunos análisis realizados. Comenzaremos descargando los datos del 1 de julio del 2014 al 30 de junio de 2024, que son 10 años de datos. Los datos se establecen de forma diaria pero podemos ver que hay unas 20 entradas por mes. Además comprobamos si hay valores faltantes dentro de los datos proporcionados o días duplicados. Tras ejecutar el código nos encontramos que no tenemos valores faltantes ni duplicados. En las siguientes figuras vemos la gráfica de los precios de cierre y después agrupamos los datos por mes y observamos cuántos datos hay en cada mes. Como podemos observar en las figura 9b el único mes que tiene ligeramente menos de 200 datos es febrero, esto se debe simplemente a que febrero tiene 28 días. Además podemos ver que en total tenemos unos 2400 entradas por lo que no tenemos un gran número de ellas pero añadir más datos puede no ser relevante debido a que no tiene sentido para la predicción de las acciones. Por eso, el enfoque de la programación probabilística puede ser útil en este caso para poder lidiar con la incertidumbre. A continuación mostraremos un análisis descriptivo para entender entre qué valores están nuestras características. Nuestro objetivo es predecir el precio de cierre. 4.2. Presentación de los datos y preprocesamiento de datos 37 (a) Precio de cierre (b) Datos por mes Figura 9: Análisis de las acciones de Microsoft Valor apertura Maximo Minimo Cierre Volumen count 2516 2516 2516 2516 2516 mean 168.26 169.890302 166.576423 168.32 29482345.07 std 113.24 114.29 112.13 113.28 13758461.80 min 40.34 40.74 39.72 40.29 7425600 25 % 62.68 63.05 62.23 62.62 21250600 50 % 134.88 135.75 133.36 134.67 26357100 75 % 258.85 261.37 255.80 258.84 33600175 max 453.07 456.17 451.77 452.85 202522400 Además añadiremos tres columnas nuevas. Una columna con la media con una ventana deslizante de 20 días que lo que te aporta es la media en los últimos 20 días, la exponential_mean de 20 días que calcula la media en función del tiempo de los últimos 20 días ponderando más aquel valor que es más cercano en el tiempo y la diferencia entre el precio máximo y el precio mínimo. Además debido a que tenemos la diferencia entre ambos eliminaremos el precio máximo y el precio mínimo. Ahora podemos observar la matriz de correlación de las variables. Como podemos 38 Capítulo 4. Aplicaciones ver en la figura 10, todas tienen una alta correlación con el precio de cierre excepto el volumen que tiene una correlación muy cercana al 0 por lo que hemos decidido eliminarla debido a que no influirá en el precio de cierre. Figura 10: Matriz de correlación Una vez realizado esto definiremos nuestro modelo de red neuronal bayesiana implementamos la red neuronal bayesiana basándonos en Back y Keith (2019). Una vez que hemos terminado de describir los datos, procedemos a dividir el conjunto en dos partes: las entradas, que llamaremos X, y la salida, que llamaremos Y. Las entradas Xincluyen todas las variables menos el precio de cierre, mientras que la salida Yserá únicamente el precio de cierre. Para mejorar el rendimiento del algoritmo y ayudar a que converja más rápido, normalizamos tanto las entradas como la salida. Esto significa que ajustamos los valores de cada característica por su normal y su desviación típica para que sigan una escala común. A continuación, tomamos los datos de los últimos 100 días de la serie temporal y los reservamos para intentar predecir su comportamiento. El resto de los datos los usaremos para entrenar nuestro modelo. Este enfoque nos permite tener una base sólida para evaluar la precisión de nuestras predicciones. El modelo que utilizamos es un modelo de una capa con una red neuronal bayesiana compuesta por 5 neuronas. Los pesos de esta red neuronal serán nuestras variables latentes, es decir, valores que inicialmente desconocemos pero que aprenderemos durante el entrenamiento. Comenzamos con una distribución previa para estos pesos, que en este caso es una distribución de Laplace con media 0 y desviación estándar 1. A medida que entrenamos el modelo, los pesos van ajustándose basados en los datos, dando lugar a una distribución posterior. Para construir este modelo en NumPyro definiremos los pesos de la siguiente manera: w1 =numpyro.sample("w1", numpyro.distributions.Laplace (jnp.zeros((dim_X, num_neurons)), jnp.ones((dim_X, num_neurons)))) donde dim_Xes el número de columnas de Xynum_neurons es el número de neuronas en cada capa oculta. De esta forma, definiremos el número de pesos que 4.2. Presentación de los datos y preprocesamiento de datos 39 necesitamos con sus respectivas distribuciones previas. Después, multiplicaremos las muestras de w1por la entrada y le aplicaremos la función no lineal para obtener la salida de esta capa. Este proceso lo repetiremos dependiendo del número de capas ocultas para obtener la salida final, en nuestro caso solo 1 vez más. La salida de esta red neuronal será interpretada como la media de una distribución normal, que hemos escogido como función de verosimilitud para nuestros datos. Elegimos la normal por ser una opción común en muchos ejemplos previos, ya que asume que los errores en nuestras predicciones se distribuyen de manera aproximadamente normal. La desviación estándar de esta distribución normal no es fija, sino que la estimamos. Partimos de una distribución previa Gamma(3,1), y la desviación estándar será calculada como el inverso de esta estimación. Este enfoque nos permite ajustar tanto los pesos como la incertidumbre en nuestras predicciones. Luego, utilizamos los datos reales observados para actualizar los pesos de la red neuronal y ajustar la desviación estándar. Este proceso, que se realiza iterativamente, mejora la precisión del modelo al comparar sus predicciones con la realidad. Una vez definido el modelo, realizamos la inferencia utilizando 30,000 muestras, de las cuales las primeras 5,000 se descartan (quemado) para asegurarnos de que el modelo haya realizado lo suficiente, y las 25,000 restantes las utilizamos para la toma de decisiones, sin embargo, iremos cogiendo 1 de cada 10 para que tengan una menor correlación dejándonos con 2,500 muestras. El método de inferencia que usamos es NUTS (No-U-Turn Sampler), un algoritmo que es muy eficiente para este tipo de problemas. Como resultado, obtenemos 2,5000 muestras de las distribuciones posteriores de cada parámetro del modelo. Estas muestras las utilizamos para ejecutar muchas versiones del modelo, en concreto para cada entrada ejecutamos el algoritmo con todas las muestras diferentes. Es decir, no ejecutamos siempre el mismo modelo fijo, sino que permitimos que los pesos cambien aleatoriamente según las distribuciones obtenidas en la inferencia, lo que introduce variabilidad en las predicciones. De todas estas ejecuciones, nuestro resultado final será la media de todas las predicciones realizadas, y también podemos calcular intervalos de confianza extrayendo los cuantiles de esas predicciones. Por ejemplo, si tomamos el cuantil 5 % y el cuantil 95 %, podemos afirmar con un 90% de confianza en qué rango es probable que se encuentre el valor real del precio de cierre. En la figura 11, podemos observar los resultados de nuestras predicciones. La línea naranja nos muestra nuestra predicción que representa la media de las predicciones.Como se puede ver, la predicción es bastante acertada, y el modelo capta bien la incertidumbre, dándonos no solo un valor puntual, sino también un intervalo de confianza del 90 % que refleja la incertidumbre relativa al problema (el sombreado azul). Lo más interesante de este ejercicio no es solo la precisión de la predicción, sino cómo la programación probabilística nos permite modelar la incertidumbre de manera efectiva, brindándonos un intervalo de confianza. En el caso de la predicción de acciones en bolsa, estamos ante un problema extremadamente complejo, y obtener predicciones con un 90 % de precisión es un gran desafío. Aun así, la capacidad de incluir la incertidumbre y generar predicciones con un rango de confianza es una de las principales ventajas de la programación 40 Capítulo 4. Aplicaciones probabilística en contextos tan inciertos como este. Por último, cabe destacar que el tiempo de ejecución es considerablemente elevado para una implementación sencilla como la presentada, la cual tardó 30 minutos. Sin embargo, hemos realizado pruebas en las que, al incrementar mínimamente el número de capas ocultas, neuronas en la red o número de cadenas, el tiempo de ejecución aumentaba exponencialmente. Esto se debe a que un mayor número de parámetros implica más entrenamiento y los métodos de inferencia tienen un alto costo computacional. Figura 11: Resultado de la predicción Cap´ ıtulo 5 Conclusiones y Trabajo Futuro Conclusión A lo largo del proyecto, hemos observado que la programación probabilística se presenta como una herramienta nueva pero de gran utilidad en el campo de la modelización. Su principal ventaja está en la posibilidad de modelar la incertidumbre dentro de cualquier problema, desde los modelos más simples hasta los más complejos. Debido a que la incertidumbre está presente en prácticamente todos los ámbitos, contar con herramientas que la gestionen y modelen de manera eficiente es fundamental. En este sentido, la programación probabilística puede ser la solución, ya que permite representar la variabilidad y las posibles distribuciones de los datos. NumPyro, un lenguaje bastante nuevo en el ámbito de la programación probabilística, destaca especialmente por su integración con JAX y sus manejadores de efectos, lo que lo hace una opción potente y novedosa. Su lenguaje simple y su eficiencia debido al soporte de JAX hacen de NumPyro una opción en desarrollo muy interesante. Además tanto en este trabajo como en el TFG de Matemáticas, donde se realizó un modelo simple de regresión lineal, hemos visto que aporta un nuevo enfoque y funciona bien para modelos sencillos. Sin embargo, para implementar modelos más complejos, como las redes neuronales bayesianas, es necesario profundizar en el estudio y desarrollo de estas herramientas debido a que al aplicar estos modelos probabilísticos a situaciones más complejas se ha observado que la convergencia se vuelve bastante compleja, especialmente cuando el número de parámetros del modelo aumenta o cuando el propio modelo implementado es difícil de ajustar. Durante el entrenamiento realizado, se exploraron diversas estrategias, como aumentar el número de capas ocultas, incrementar el número de muestras o añadir más neuronas a la red y uno de los problemas principales ha sido la dificultad en alcanzar la convergencia y el alto costo computacional requerido para lograrla. A pesar de esto, los resultados obtenidos en la red neuronal bayesiana son esperanzadores para continuar desarrollando esta herramienta y que la programación probabilística, y en particular NumPyro, sean una herramienta con un enorme potencial. 41