I am trying to model a process, which looks a bit like a function, with heavy-tailed noise, and some values are also multiplied by a random value between 0 and 1. Consider the following example data.
rng = np.random.default_rng(4) x = np.linspace(-3, 3, 50) y = 5 + np.sin(x) noise = 0.2 * rng.standard_t(df=1.5, size=x.shape) # multiply 20% of the values by a random number from 0-1 eff = np.minimum(5 * rng.random(size=x.shape), 1) obs = (y + noise) * eff
I can leverage
StudentT to capture the underlying function quite well. As per the following:
with pm.Model() as model: # Underlying GP ls = pm.Gamma("ls", mu=3, sigma=1) scale = pm.Gamma("scale", mu=1, sigma=0.5) cov = scale ** 2 * pm.gp.cov.ExpQuad(input_dim=1, ls=ls) gp = pm.gp.Latent(cov_func=cov) f = gp.prior("f", x[:, None]) # Noise offset = pm.Normal("offset", mu=5, sigma=2) nu = 1 + pm.Gamma("noise_nu", mu=1, sigma=1) sigma = pm.HalfNormal("noise_sigma", sigma=1) pm.StudentT("noisy_f", nu=nu, mu=offset + f, sigma=sigma, observed=obs) # Sample trace = pm.sample() x_new = np.linspace(-4, 4, 100) with model: f_pred = gp.conditional("f_pred", x_new[:, None]) pm.StudentT("noisy_f_pred", nu=nu, mu=offset + f_pred, sigma=sigma) ppc = pm.sample_posterior_predictive(trace, var_names=["f_pred", "noisy_f_pred"])
However, though this approach captures the underlying function well, it overestimates the noise variance by quite a bit and obviously doesn’t capture the random drops in “efficiency” (i.e. the multiplications by 0-1)
How can I model this process more accurately?
I’ve tried a couple of things but they always lead to strange or completely wrong results. I’m interested in creating a generative model so I can sample other similar-looking data.
SIDE-NOTE: I am loving the new
HSGP class and the ability to leverage
prior_linearized to treat GPs like a regular linear model and get rapid inference…but didn’t want to complicate my example above by including it.