Dependencia de la observación en los modelos pymc3

Question

Tengo un modelo, que se puede simplificar conceptualmente a:

$$ a \sim TruncNormal(\mu = 1.0, \sigma=0.01, min = 0.9, max = 1.1)$$

$$y = a \cdot sin(b) $$

Puedo hacer observaciones sobre $y$ pero todas estas observaciones tienen su propia y única $b$ ángulo. Denotemos una muestra de $a$ como $a[i]$ . $$y_{measured}[1] = a[1] \cdot sin(b_1)$$ $$y_{measured}[2] = a[2] \cdot sin(b_2)$$ $$y_{measured}[2] = a[3] \cdot sin(b_3)$$ $$ ... $$ $$y_{measured}[n] = a[n] \cdot sin(b_n)$$

Estos ángulos están en el rango de $(-\pi,\pi)$ y se eligen explícitamente antes de realizar una medición. Más explícitamente, podemos elegir medir en

$b_1 = \frac{-pi}{2}$ entonces podemos observar $y_{measured} = 1.01$

podemos hacer otra medición en

$b_2 = \frac{-pi}{4}$ entonces podemos observar $y_{measured} = 0.717$ .

Obviamente, si no considero el derecho $b_i$ ángulo para el derecho $y_{measured}$ al inferir un ello llevaría a un sinsentido.

Estoy tratando de implementar esto en python usando pymc3.

import pymc3 as pm
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import math

a_data = 1.042
b_data = np.arange(-math.pi, math.pi, 0.1)
exact_measurement_data = a_data * np.sin(b_data)
measurement = exact_measurement_data + np.random.normal(exact_measurement_data, 0.01, b_data.size)

with pm.Model() as mdl:
   a = pm.TruncatedNormal('a', mu = 1.0, sd = 0.01, lower = 0.9, upper = 1.1)
   b = b_data # ??? Is this correct?
   y = pm.Deterministic('y', a * pm.math.sin(b))
   y_measured = pm.Normal('y_measured', mu = y, sd = 0.01, observed = measurement)

Como se puede ver, aún no está completo. Mi pregunta sería cómo definir $b$ y cómo inferir $a$ .

$b$ puede ser un valor exacto, o puedo definir una prioridad para $b$ . Esta prioridad sería, por ejemplo, una distribución normal en torno a todos y cada uno de los elementos de la lista, y equivaldría a un error de medición del ángulo.

Mi problema se debe al hecho de que mi observación para $y_{measured}[i]$ depende de $b[i]$ . Más explícitamente, la medida para $y[i]$ está en función de $b$ Puedo medir un valor para $y$ sólo para un determinado valor de $b[i]$ . La gama para $b[i]$ es en el peor de los casos $(-\pi,\pi)$ pero me parece bien no tener $b = \frac{\pi}{2}$ o valores problemáticos similares en las observaciones. Mi objetivo con el modelo es inferir $a$ , en base a las mediciones. ¿Se puede hacer esto en pymc3?

Answer 1

Lo que describes difiere del código que has proporcionado. Afortunadamente, has proporcionado un ejemplo reproducible, así que lo trataré como una descripción inequívoca de tu problema.

Suponiendo que sus datos sean como los siguientes

import math
import numpy as np

np.random.seed(42)

a_data = 1.042
b_data = np.arange(-math.pi, math.pi, 0.1)
exact_measurement_data = a_data * np.sin(b_data)
measurement = exact_measurement_data + np.random.normal(exact_measurement_data, 0.01, b_data.size)

En este caso, se puede describir utilizando la notación matemática como

$$\begin{align} \mu_i &= a \sin(b_i) \\ y_i &\sim \mathcal{N}(\mu_i,\, 0.01) \end{align}$$

Obsérvese que aquí he suprimido los sufijos "data" y "measured" para simplificar la notación, $b$ es b_data y $y$ es measurement en su ejemplo. Donde tienes pares de $(y_i, b_i)$ puntos de datos, y un único $a$ parámetro. En su descripción, ha utilizado el ambiguo $x_i$ y $x[i]$ subíndices que no coinciden con el ejemplo, por lo que es difícil comentarlos.

La descripción matemática se traduce directamente en el código de PyMC3. Al igual que con otros lenguajes de programación, su objetivo es ajustarse a la notación lo más posible dadas las restricciones relacionadas con los lenguajes de programación en los que se basan.

import pymc3 as pm

with pm.Model() as mdl:
   a = pm.TruncatedNormal('a', mu = 1.0, sd = 0.01, lower = 0.9, upper = 1.1)
   mu = pm.Deterministic('y', a * pm.math.sin(b_data))
   y_measured = pm.Normal('y_measured', mu = mu, sd = 0.01, observed = measurement)