1
$\begingroup$

I have a model, which can be simplified conceptually to:

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

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

I can make observations about $y$, but these observations all have their own, unique $b$ angle. Let us denote a sample from $a$ as $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)$$

These angles are in the range from $(-\pi,\pi)$ and are explicitly chosen before taking a measurement. More explicitly, we can choose to measure at

$b_1 = \frac{-pi}{2}$ then we can observe $y_{measured} = 1.01$

we can make another measurement at

$b_2 = \frac{-pi}{4}$ then we can observe $y_{measured} = 0.717$.

Obviously if I do not consider the right $b_i$ angle for the right $y_{measured}$ when inferring a it would lead to a nonsense.

I am trying to implement this in python using 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)

As it can be seen, it is not yet complete. My question would be how do define $b$ and how to infer $a$.

$b$ can be an exact value, or I can define a prior for $b$. This prior would be e.g. a normal distribution around each and every element in the list.It would be the equivalent of a measurement error of the angle.

My problem is due to the fact, that my observation for $y_{measured}[i]$ is dependent on $b[i]$. More explicitly, the measurement for $y[i]$ is in function of $b$, I can measure a value for $y$ only for a certain value of $b[i]$. The range for $b[i]$ is in the worst case $(-\pi,\pi)$, but I am OK with not having $b = \frac{\pi}{2}$ or similar problematic values in the observations. My goal with the model is to infer $a$, based on the measurements. Can this be done in pymc3?

$\endgroup$
6
  • $\begingroup$ It’s not clear what is the problem? What is wrong with just setting a prior for $b$? $\endgroup$
    – Tim
    Jun 27 at 12:14
  • $\begingroup$ I do not know how to correlate between observations for y and b. I have basically two arrays of values. Should i just add a uniform prior for b for the whole range? Does the index in the observations for b and y habe an intrinsic correlation in pymc3? $\endgroup$
    – 50k4
    Jun 27 at 12:23
  • $\begingroup$ Correlate what? From your description it sounds as just another parameter in the model. $\endgroup$
    – Tim
    Jun 27 at 12:40
  • $\begingroup$ It is. But every y in my observations is dependent on a b angle. If I just add a list of b observations as a list, does pymc „pair“ the observations for b with the observations for y based on their index in the list? $\endgroup$
    – 50k4
    Jun 27 at 12:43
  • $\begingroup$ Pair what exactly? In your code $b$ is a parameter, $y$ is a vector of data points. In a frequentist setting, $b$ would be a single number while $y$ would be multiple numbers, so there is nothing to "pair". $\endgroup$
    – Tim
    Jun 27 at 20:28
0
$\begingroup$

What you describe differs from the code that you provided. Hopefully, you provided a reproducible example, so I'll treat it as an unambiguous description of your problem.

Assuming that your data is something like below

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)

In such a case, this can be described using mathematical notation as

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

Notice that I dropped here the "data" and "measured" suffixes to simplify the notation, $b$ is b_data and $y$ is measurement in your example. Where you have pairs of $(y_i, b_i)$ datapoints, and a single $a$parameter. In your description, you used the ambiguous $x_i$ and $x[i]$ subscripts that do not match the example, so it is hard to comment on them.

The math description translates directly to the PyMC3 code. As with other programming languages, their aim is to match the notation as close as possible given the restrictions related to the programming languages they are based on.

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)
$\endgroup$
3
  • $\begingroup$ This pretty much looks like my code, with less unnecessary assignments, like b = b_data. Do I infer a with pm.sample()? $\endgroup$
    – 50k4
    Jun 28 at 7:06
  • $\begingroup$ @50k4 pm.sample if you want to sample from the posterior. I'd encourage you to start with the tutorials or a book $\endgroup$
    – Tim
    Jun 28 at 7:08
  • $\begingroup$ Thanks! books and tutorials got me to the code I posted. $\endgroup$
    – 50k4
    Jun 28 at 7:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.