<jupyter_start><jupyter_text>Trigonometric time series model<jupyter_code>%pylab inline
import pymc3 as pm
import theano.tensor as tt
t = np.linspace(0., 10., 1000)
# 3.5Hz sine wave
func = lambda a, b, omega: a*np.sin(2*np.pi*omega*t)+b*np.cos(2*np.pi*omega*t)
y = func(1., 0., 3.5)
data = y + np.random.normal(size = t.shape[0])

_, ax = plt.subplots(1, 1, figsize=(15, 3))
ax.plot(t, data)
ax.plot(t, y);
with pm.Model() as model:
    a = pm.Normal("a", mu=1, sd=2)
    b = pm.Normal("b", mu=1, sd=2)
    omega = pm.Gamma("omega", 1., 1.)
    regression = a * tt.sin(2 * np.pi * omega * t) +\
                 b * tt.cos(2 * np.pi * omega * t)
    sd = pm.HalfCauchy("sd", 0.5)
    observed = pm.Normal("observed",
                         mu=regression, 
                         sd=sd, 
                         observed=data)
with model:
    step1 = pm.Metropolis([omega])
    trace = pm.sample(10000, tune=5000, step=step1)
pm.traceplot(trace);
omegapost = trace.get_values(varname='omega', combine=False)
apost = trace.get_values(varname='a', combine=False)
bpost = trace.get_values(varname='b', combine=False)
_, ax = plt.subplots(1, 1, figsize=(15, 3))
ax.plot(t, data, alpha=.5)
ax.plot(t, y, color='k', alpha=.5)
for ia, ib, io in zip(apost, bpost, omegapost):
    y_ = func(ia.mean(), ib.mean(), io.mean())
    ax.plot(t, y_,  alpha=.5);
with model:
    trace2 = pm.sample(1000, tune=1000)
pm.traceplot(trace2);
pm.pairplot(trace2, alpha=.05);
model.free_RVs
logp_dlogp_cond = model.logp_dlogp_function([model.free_RVs[0], 
                                             model.free_RVs[2]])
pt = model.test_point
pt
pt['b'] = np.array(0.)
pt['sd_log__'] = np.log(0.5)
logp_dlogp_cond.set_extra_values(pt)
pt['a'] = np.array(2.)
pt
logp_dlogp_cond.dict_to_array(pt)
pt['b'] = np.array(0.)
pt['sd_log__'] = np.log(1.)
logp_dlogp_cond.set_extra_values(pt)
# grid search
omega_ = np.linspace(0, 4, 200)
a_ = np.linspace(-4, 4, 200)
ov_, av_ = np.meshgrid(omega_, a_)
logvec = np.asarray([logp_dlogp_cond(np.asarray([a__, np.log(o__)]))[0]
                    for o__, a__ in zip(ov_.flatten(), av_.flatten())])
ll = logvec.reshape(av_.shape)

_, ax = plt.subplots(1, 1, figsize=(10, 5))
ax.imshow(ll.T, cmap='viridis',
         extent=[a_[0], a_[-1], omega_[0], omega_[-1]], origin='lower'
         )
ax.set_xlabel('a')
ax.set_ylabel('omega');
pt['b'] = np.array(0.)
pt['sd_log__'] = np.log(1.)
logp_dlogp_cond.set_extra_values(pt)
# grid search
omega_ = np.linspace(3, 4, 200)
a_ = np.linspace(-4, 4, 200)
ov_, av_ = np.meshgrid(omega_, a_)
logvec = np.asarray([logp_dlogp_cond(np.asarray([a__, np.log(o__)]))[0]
                    for o__, a__ in zip(ov_.flatten(), av_.flatten())])
ll = logvec.reshape(av_.shape)

_, ax = plt.subplots(1, 1, figsize=(10, 5))
ax.imshow(ll.T, cmap='viridis',
         extent=[a_[0], a_[-1], omega_[0], omega_[-1]], aspect=4, origin='lower'
         )
ax.set_xlabel('a')
ax.set_ylabel('omega');
pt['b'] = np.array(0.)
pt['sd_log__'] = np.log(1.)
logp_dlogp_cond.set_extra_values(pt)
# grid search
omega_ = np.linspace(0, .25, 200)
a_ = np.linspace(-4, 4, 200)
ov_, av_ = np.meshgrid(omega_, a_)
logvec = np.asarray([logp_dlogp_cond(np.asarray([a__, np.log(o__)]))[0]
                    for o__, a__ in zip(ov_.flatten(), av_.flatten())])
ll = logvec.reshape(av_.shape)

_, ax = plt.subplots(1, 1, figsize=(10, 5))
ax.imshow(ll.T, cmap='viridis',
         extent=[a_[0], a_[-1], omega_[0], omega_[-1]], aspect=16, origin='lower'
         )
ax.set_xlabel('a')
ax.set_ylabel('omega');
pt['b'] = np.array(0.)
pt['sd_log__'] = np.log(1.)
logp_dlogp_cond.set_extra_values(pt)
# grid search
omega_ = np.linspace(0, .25, 200)
a_ = np.linspace(-10, 10, 200)
ov_, av_ = np.meshgrid(omega_, a_)
logvec = np.asarray([logp_dlogp_cond(np.asarray([a__, np.log(o__)]))[0]
                    for o__, a__ in zip(ov_.flatten(), av_.flatten())])
ll = logvec.reshape(av_.shape)

_, ax = plt.subplots(1, 1, figsize=(10, 5))
ax.imshow(ll.T, cmap='viridis',
         extent=[a_[0], a_[-1], omega_[0], omega_[-1]], aspect=16*2, origin='lower'
         )
ax.set_xlabel('a')
ax.set_ylabel('omega');<jupyter_output><empty_output><jupyter_text>GP<jupyter_code>with pm.Model() as model:
    η_per = pm.HalfCauchy("η_per", beta=2, testval=1.0)
    period = pm.Normal("period", mu=1 / 3.5, sd=1.)
    periodic = η_per**2 * pm.gp.cov.Cosine(1, period)

    gp_periodic = pm.gp.Marginal(cov_func=periodic)

    noise = pm.HalfCauchy("noise", beta=2, testval=1.0)
    y_ = gp_periodic.marginal_likelihood(
        "y", X=t[:, np.newaxis], y=data, noise=noise)
mp = model.test_point
mp
with model:
    mp = pm.find_MAP()
mp
mp2 = mp
# mp2['noise_log__'] = np.log(1)
# mp2['η_per_log__'] = np.log(1)
# mp2['period'] = np.array(1/3.5)

mu_pred, cov_pred = gp_periodic.predict(t[:, np.newaxis], point=mp2)
_, ax = plt.subplots(1, 1, figsize=(15, 3))
ax.plot(t, data, alpha=.5)
ax.plot(t, y, color='k', alpha=.5)
ax.plot(t, mu_pred)
model.free_RVs
logp_dlogp_cond = model.logp_dlogp_function([model.free_RVs[0], model.free_RVs[1]])
pt = model.test_point
pt
logp_dlogp_cond.set_extra_values(mp)
logp_dlogp_cond.dict_to_array(pt)
# grid search
omega_ = np.linspace(0, 4, 100)
a_ = np.linspace(0, 1, 100)
ov_, av_ = np.meshgrid(omega_, a_)
logvec = np.asarray([logp_dlogp_cond(np.asarray([np.log(a__), 1/o__]))[0]
                    for o__, a__ in zip(ov_.flatten(), av_.flatten())])
ll = logvec.reshape(av_.shape)

_, ax = plt.subplots(1, 1, figsize=(10, 5))
ax.imshow(ll.T, cmap='viridis',
         extent=[a_[0], a_[-1], omega_[0], omega_[-1]], aspect=1/16*2, origin='lower'
         )
ax.set_xlabel('a')
ax.set_ylabel('omega');<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1996822033898305, 'fraction_numerical': 0.03954802259887006, 'mean_word_length': 3.7685185185185186, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '26960', 'license': 'MIT', 'n_tokens_mistral': 2576, 'n_tokens_neox': 2438, 'n_words': 439, 'path': 'Notebooks/Timeserie_model.ipynb', 'repo_name': 'junpenglao/All-that-likelihood-with-PyMC3', 'size': '5664'}