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[0])
# multiply 20% of the values by a random number from 0-1
eff = np.minimum(5 * rng.random(size=x.shape[0]), 1)
obs = (y + noise) * eff
```

I can leverage `gp.Latent`

and `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.