# Simple GP model giving large MAP values for covariance timescale

**URL:** <https://discourse.pymc.io/t/simple-gp-model-giving-large-map-values-for-covariance-timescale/7844>\
**Category:** Questions\
**Created:** [August 4, 2021, 11:17pm UTC](https://discourse.pymc.io/t/simple-gp-model-giving-large-map-values-for-covariance-timescale/7844 "2021-08-04T23:17:01Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![burke86](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/burke86/32/4323_2.png) [@burke86](https://discourse.pymc.io/u/burke86)\
**Post date:** [August 4, 2021, 11:17pm UTC](https://discourse.pymc.io/t/simple-gp-model-giving-large-map-values-for-covariance-timescale/7844/1 "2021-08-04T23:17:01Z")

</div>

I’m trying to fit a simple GP model with exponential covariance (damped random walk model) with timescale `l`. However, the MAP values do not seem reasonable, and fail to recover the input values from simulated time series.

The code for the model is:

```auto
def fit_drw(t, y, yerr, cadence, baseline, amplitude, precision):
        
        import pymc3 as pm
        import numpy as np
        import matplotlib.pyplot as plt
        
        with pm.Model() as model:
            
            # damped random walk
            l = pm.Uniform("l", lower=np.sqrt(1/2), upper=np.sqrt(1e8*baseline/2))
            # 2l^2 = tau
            # l = sqrt(tau/2)
            sigma_drw = pm.Uniform("sigma_drw", lower=0.1*precision, upper=10*amplitude)
            cov = 2*sigma_drw**2 * pm.gp.cov.Exponential(1, l)
            gp_drw = pm.gp.Marginal(cov_func=cov)

            # The Gaussian process is a sum of these three components
            gp = gp_drw 

            # Since the normal noise model and the GP are conjugates, we use `Marginal` with the `.marginal_likelihood` method
            X = t[:, None]
            y_ = gp.marginal_likelihood("y", X=X, y=y, noise=yerr)
            mp = pm.find_MAP() #start={'l': 100, 'sigma_drw': amplitude})
            
            # Predict
            tpred = np.linspace(np.min(t), np.max(t)+400, 1000)
            Xpred = tpred[:, None]
            
            mu, var = gp.predict(Xpred, point=mp, diag=True)
            sd = np.sqrt(var)
            
            plt.figure(figsize=(8,4))
            plt.plot(tpred, mu, "dodgerblue", lw=3)
            plt.fill_between(tpred, mu-sd, mu+sd, color="dodgerblue", alpha=0.2)
            plt.errorbar(t, y, yerr=yerr, color="k", linestyle='none', ms=3, alpha=1)
            plt.show()
                        
            print(mp)
            
        return

```

The MAP prediction looks reasonable when plotted below. However, the MAP value of `l` is huge! Why?

```auto
{'l_interval__': array(12.60262633), 'sigma_drw_interval__': array(1.12011405), 'l': array(499998.31842261), 'sigma_drw': array(8.61098201)}

```

 ![Screen Shot 2021-08-04 at 6.09.25 PM](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/e/e9e728c779d12fbd171264a8c6e67d2865c9e680.jpeg)

For completion, the coed to simulate the time series is:

```auto
from astroML.time_series import generate_damped_RW

dt = 50
baseline = 5000
t = np.arange(0, baseline, dt)
y = generate_damped_RW(t, tau=250, xmean=20, SFinf=0.3, z=0.0)

fit_drw(t, y, 0.005*y, cadence=dt, baseline=baseline, amplitude=np.max(y)-np.min(y), precision=0.005)

```

The details of the generate\_damped\_RW function should not be a concern, other than the input timescale is given by `2l^2 = tau`.

Thanks!

---

<div class="post-metadata">

**Author:** ![burke86](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/burke86/32/4323_2.png) [@burke86](https://discourse.pymc.io/u/burke86)\
**Post date:** [August 4, 2021, 11:42pm UTC](https://discourse.pymc.io/t/simple-gp-model-giving-large-map-values-for-covariance-timescale/7844/2 "2021-08-04T23:42:33Z")

</div>

Note I had to make the upper limit on the `l` prior very large to get a reasonable prediction (otherwise the MAP value of `l` was just giving the upper-limit).

---

<div class="post-metadata">

**Author:** ![jlindbloom](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jlindbloom/32/5016_2.png) [@jlindbloom](https://discourse.pymc.io/u/jlindbloom)\
**Post date:** [August 5, 2021, 2:00am UTC](https://discourse.pymc.io/t/simple-gp-model-giving-large-map-values-for-covariance-timescale/7844/3 "2021-08-05T02:00:13Z")

</div>

I messed around with your code and found the same thing. I have a feeling that you’re experiencing the same issue that Michael Betancourt talks about in [this blog post](https://betanalpha.github.io/assets/case_studies/gp_part3/part3.html#2_inferring_the_hyperparameters_with_a_naive_bayesian_fit). See [this thread](https://discourse.pymc.io/t/how-to-put-a-restriction-on-the-value-of-hyper-parameters-in-gaussian-process/6072/2) for some advice. Long story short, priors are important (but especially with GPs)!

---

<div class="post-metadata">

**Author:** ![burke86](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/burke86/32/4323_2.png) [@burke86](https://discourse.pymc.io/u/burke86)\
**Post date:** [August 5, 2021, 2:47am UTC](https://discourse.pymc.io/t/simple-gp-model-giving-large-map-values-for-covariance-timescale/7844/4 "2021-08-05T02:47:29Z")

</div>

Thanks @jlindbloom,

I tried the InverseGamma prior (using the minimum sampling rate and time series length as the lower and upper length scales) and Normal prior for the amplitude term. Now we are getting somewhere, but the recovered `l` (or `tau`) seem systematically larger (by a value of ~1e5), and much larger than the time series length! Do you have any ideas what could cause this difference by a multiplicative factor?

Below I show the input timescale vs. recovered MAP timescale after dividing by 1e5:

 ![Screen Shot 2021-08-04 at 9.46.12 PM](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/b/b2b147c84a62e7d4b0ef0469e8e860bd9c1043bc.png)

---

<div class="post-metadata">

**Author:** ![jlindbloom](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/jlindbloom/32/5016_2.png) [@jlindbloom](https://discourse.pymc.io/u/jlindbloom)\
**Post date:** [August 5, 2021, 2:56am UTC](https://discourse.pymc.io/t/simple-gp-model-giving-large-map-values-for-covariance-timescale/7844/5 "2021-08-05T02:56:49Z")

</div>

What precisely do you mean by “input” \tau? Do you mean the true timescale used to generate your training data?

Can you share your code for this?

---

<div class="post-metadata">

**Author:** ![burke86](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/burke86/32/4323_2.png) [@burke86](https://discourse.pymc.io/u/burke86)\
**Post date:** [August 5, 2021, 3:01am UTC](https://discourse.pymc.io/t/simple-gp-model-giving-large-map-values-for-covariance-timescale/7844/6 "2021-08-05T03:01:01Z")

</div>

Yes, input \tau is the the timescale of the simulated time series. Code below:

```auto
def fit_drw(t, y, yerr, cadence, baseline, amplitude, precision):
        
        import pymc3 as pm
        import numpy as np
        import matplotlib.pyplot as plt
                
        l_lower = cadence
        l_upper = baseline
        l_sigma = (l_upper-l_lower)/6
        l_mu = l_lower + 3*l_sigma
        
        with pm.Model() as model:
            
            # damped random walk
            tau = pm.InverseGamma("tau_drw", mu=l_mu, sigma=l_sigma)
            # 2l^2 = tau_drw
            # l = sqrt(tau/2)
            sigma_drw = pm.Normal("sigma_drw", mu=amplitude, sigma=precision)
            cov = 2*sigma_drw**2 * pm.gp.cov.Exponential(1, np.sqrt(tau/2))
            gp_drw = pm.gp.Marginal(cov_func=cov)

            # white noise
            sigma_n = pm.Normal("sigma_n", mu=precision, sigma=precision)
            cov = pm.gp.cov.WhiteNoise(sigma_n)
            gp_wn = pm.gp.Marginal(cov_func=cov)

            # The Gaussian process is a sum of these two components
            gp = gp_drw + gp_wn

            # Since the normal noise model and the GP are conjugates, we use `Marginal` with the `.marginal_likelihood` method
            X = t[:, None]
            y_ = gp.marginal_likelihood("y", X=X, y=y, noise=0)
            mp = pm.find_MAP(start={'tau_drw': 0.25*baseline, 'sigma_drw': amplitude, 'sigma_n':precision})
            
            # Predict
            tpred = np.linspace(np.min(t), np.max(t), 1000)
            Xpred = tpred[:, None]
            
            mu, var = gp.predict(Xpred, point=mp, diag=True)
            sd = np.sqrt(var)
            
            plt.figure(figsize=(12,4))
            plt.plot(tpred, mu, "dodgerblue", lw=3)
            plt.fill_between(tpred, mu-sd, mu+sd, color="dodgerblue", alpha=0.2)
            plt.errorbar(t, y, yerr=yerr, color="k", linestyle='none', ms=3, alpha=1)
            plt.show()
        
        return mp

```

```auto
from astroML.time_series import generate_damped_RW

dt = 50
baseline = 5000
t = np.arange(0, baseline, dt)

# Simulate time series at varying tau
tau_outs = []
tau_ins = np.linspace(50, 500, 20)
for tau_in in tau_ins:

    y = generate_damped_RW(t, tau=tau_in, xmean=20, SFinf=0.3, z=0.0)

    noise_level = 0.005
    y += np.random.normal(noise_level)

    mp = fit_drw(t, y, noise_level*y, cadence=dt, baseline=baseline, amplitude=np.max(y)-np.min(y), precision=noise_level)
    tau_outs.append(mp['tau_drw'])

```

```auto
# Plot input versus MAP tau
plt.scatter(tau_ins, np.array(tau_outs)/1e5)
plt.plot([0, 500], [0, 500], color='k', lw=2)
plt.xlim([0, 500])
plt.ylim([0, 500])
plt.ylabel(r'MAP $\tau / 10^5$', fontsize=16)
plt.xlabel(r'Input $\tau$', fontsize=16)

```

---

<div class="post-metadata">

**Author:** ![burke86](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/burke86/32/4323_2.png) [@burke86](https://discourse.pymc.io/u/burke86)\
**Post date:** [August 5, 2021, 3:25am UTC](https://discourse.pymc.io/t/simple-gp-model-giving-large-map-values-for-covariance-timescale/7844/7 "2021-08-05T03:25:39Z")

</div>

Update: The multiplicative “correction” factor seems to be dependent on the time series baseline (length). ~20 times the baseline seems to work. Strange. I am still investigating.

---

<div class="post-metadata">

**Author:** ![BioGoertz](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/biogoertz/32/4489_2.png) [@BioGoertz](https://discourse.pymc.io/u/BioGoertz)\
**Post date:** [August 5, 2021, 11:02am UTC](https://discourse.pymc.io/t/simple-gp-model-giving-large-map-values-for-covariance-timescale/7844/8 "2021-08-05T11:02:55Z")

</div>

With GPs, I would strongly recommend standardizing your input and output variables (subtract their mean and divide by their standard deviation) prior to fitting. I find this makes it easier to specify priors (largely because you can just use off-the-shelf suggestions rather than rescaling things) and to interpret parameter posteriors. I _think_ there’s also something computationally advantageous about picking priors that are closer to “standard” values (e.g., N(0,1) rather than N(1e5,1e4)), something about the tuning, but I could be completely off-base with that. But I do think part of your issue may arise from the fact that your input and output variables lie on vastly different scales, increments of ~100 vs 0.01. That could be why your MAP value for `sigma_drw` actually _also_ seems 1-2 orders of magnitude larger than I might expect. The interaction between GP variance and lengthscale is complex and somewhat unintuitive; technically speaking, I’m pretty sure variance and lengthscale are unidentifiable (though their ratio is identifiable), meaning there are ~infinite combinations that give identical behavior, so I wouldn’t worry about trying to reproduce their values or draw too much interpretation from their values.

Also why take the square root of tau for the lengthscale? That’s just going to force tau to be larger than necessary, and sort of defeats the purpose of Betancourt’s principled InverseGamma prior.

I’m also suspicious that what you’re seeing is partly to do with the way you are introducing noise to your simulations and specifying noise in your model. It’s not common to account for noise in the GP with a WhiteNoise kernel; typically you would specify a prior for the noise argument to `gp.marginal_likelihood`, e.g. `σ = pm.Exponential('σ', lam=1) ... gp.marginal_likelihood('ml', X=X, y=y, noise=σ)`. Those may be equivalent formulations, but I’m not sure. I noticed that in your code earlier in the thread, you specified yerr for this argument. This essentially acts as weights to the mean of the GP, and probably isn’t what you want. Additionally, you added `noise_level` to y in order to jitter those values, but you _multiplied_ it by y to specify `yerr`, which doesn’t make sense to me.

You’ll want to specify a strictly positive prior for the GP variance (your `sigma_drw`, more commonly called _eta_ , _η_, for a GP) as well as noise; for standardized data, Gamma(2,1) works well for variance, Exponential(1) for noise. Because you’re squaring it, a Normal prior will end up having a bimodal posterior around 0, which will just make it unnecessarily difficult for your sampler to converge.

* * *

tl;dr: try scaling/standardizing your inputs and outputs, which should allow you to use more standard priors, and specify noise as noise, not a separate kernel.

---

<div class="post-metadata">

**Author:** ![burke86](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/burke86/32/4323_2.png) [@burke86](https://discourse.pymc.io/u/burke86)\
**Post date:** [August 5, 2021, 11:28pm UTC](https://discourse.pymc.io/t/simple-gp-model-giving-large-map-values-for-covariance-timescale/7844/9 "2021-08-05T23:28:23Z")

</div>

Thanks for your reply @BioGoertz.

Standardizing the inputs seems like a good idea.

The motivation for taking the square root of the lengthscale is so the inferred lengthscale is the same units as the input times. However, there is nothing wrong with doing the conversion between `tau` and `l` afterwords and adopting the standard priors. It also makes the analytic PSD easy to write down. See link below if you are curious about this convention in my field:

> **[Is Quasar Optical Variability a Damped Random Walk?](https://arxiv.org/abs/1202.3783)**
>
> The damped random walk (DRW) model is increasingly used to model the
> variability in quasar optical light curves, but it is still uncertain whether
> the DRW model provides an adequate description of quasar optical variability
> across all time scales....

 ![Screen Shot 2021-08-05 at 6.17.57 PM](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/0/0fa9e6babaced1e529ceceb58c33daf35399e17d.png)

The lengthscale is not completely unidentifiable/degenerate with the variance, although the covariance between the two parameters is strong. As my figure above shows, the MAP lengthscale is proportional to the simulated input lengthscale.

Thanks for noticing the noise typo and suggestions on the priors. I’ll keep looking into this.

---

<div class="post-metadata">

**Author:** ![burke86](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/burke86/32/4323_2.png) [@burke86](https://discourse.pymc.io/u/burke86)\
**Post date:** [August 6, 2021, 4:17am UTC](https://discourse.pymc.io/t/simple-gp-model-giving-large-map-values-for-covariance-timescale/7844/10 "2021-08-06T04:17:14Z")

</div>

Going to mark this as solved. Indeed something was wrong with my priors initially. MAP timescale still seems off by a factor of a few, but that could be a bias in the MAP estimate perhaps. Also once the timescale gets too long compared to its length (input tau ~ 10 length; gray lines), the time series is not stationary and the MAP timescales saturate, which is a well-known effect. I’m not sure if there is good solution for this.

 ![Screen Shot 2021-08-05 at 11.12.44 PM](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/0/0da39f31c4e3c8a40278190958b3c5bdc3cb86c1.png)