Contents

Colored Noise

%pylab inline
%config InlineBackend.figure_format = 'svg'

import arviz as az
import celerite2
from celerite2 import terms
from celerite2.theano import terms as theano_terms
import celerite2.theano
import exoplanet as xo
import glob
import h5py
import os.path as op
import pandas as pd
import pymc3 as pm
import pymc3_ext as pmx
from scipy.interpolate import interp1d
from scipy.signal import welch, decimate
from scipy.signal.windows import tukey
import seaborn as sns
import stan_jupyter as stan
import theano.tensor as T

sns.set_context('notebook')
sns.set_style('ticks')
sns.set_palette('colorblind')

Populating the interactive namespace from numpy and matplotlib

Some random notes from Lecture

Correlated Noise

Sigma = array([[2.0, sqrt(2)*sqrt(3)*0.9], [sqrt(2)*sqrt(3)*0.9, 3]])

xs = np.random.multivariate_normal(zeros(2), Sigma, size=1000)

plot(xs[:,0], xs[:,1], '.')

[<matplotlib.lines.Line2D at 0x1494969a0>]
_images/colored_noise_5_1.svg

Here we made a Gaussian process with kernel $\( k\left( x_i, x_j \right) = \sigma^2 \exp\left( -\frac{1}{2} \left(\frac{x_i - x_j}{\lambda} \right)^2 \right). \)$ For obvious reasons, this is called the “squared exponential kernel.” A good place to learn a lot about Gaussian processes, kernel functions, etc, etc, is [Rasmussen and Williams, 2006], which can be found online here.

xs = linspace(-1, 1, 128)

sigma2 = 10.0
lambda2 = 0.01
def k(xi, xj):
    return sigma2*np.exp(-0.5*np.square(xi-xj)/lambda2)

Sigma = zeros((128,128))
for i in range(128):
    for j in range(128):
        Sigma[i,j] = k(xs[i], xs[j])

imshow(Sigma)

<matplotlib.image.AxesImage at 0x149560eb0>
_images/colored_noise_7_1.svg
for i in range(10):
    plot(xs, np.random.multivariate_normal(zeros_like(xs), Sigma))
_images/colored_noise_8_0.svg
\[ \rho(\tau) = \sigma^2 \exp\left( -\frac{\left| \tau \right|}{\lambda} \right) \]
xs = linspace(-1, 1, 128)

sigma2 = 10.0
tau = 1.0
def k(xi, xj):
    return sigma2*np.exp(-np.abs(xi-xj)/tau)

Sigma = zeros((128,128))
for i in range(128):
    for j in range(128):
        Sigma[i,j] = k(xs[i], xs[j])

imshow(Sigma)

<matplotlib.image.AxesImage at 0x1496c1bb0>
_images/colored_noise_10_1.svg
for i in range(10):
    plot(xs, np.random.multivariate_normal(zeros_like(xs), Sigma))
_images/colored_noise_11_0.svg

An Example from LIGO

I downloaded the Hanford LIGO data for GW150914 from GWOSC. There are a number of different data products there, but I grabbed the 4 kHz sample rate, 32 second long segment surrounding this event. I also grabbed the “parameter estimation” posterior samples from the same place.

with h5py.File('data/H-H1_GWOSC_4KHZ_R1-1126259447-32.hdf5', 'r') as f:
    strain = array(f['strain/Strain'])
    strain_start = array(f['meta/GPSstart'])[()]
    dt = 1.0 / 4096.0
    strain_times = linspace(0, (len(strain)-1)*dt, len(strain))

with h5py.File('data/GW150914_GWTC-1.hdf5', 'r') as f:
    postsamples = array(f['Overall_posterior'])
t_approx = 1126259462.422

plot(strain_times, strain)
axvline(t_approx - strain_start, color='k')
xlabel(r'$t / \mathrm{s}$')
ylabel(r'$h(t)$')

Text(0, 0.5, '$h(t)$')
_images/colored_noise_16_1.svg

We know the 150914 signal ends at low frequency (< 500 Hz), so we will decimate this data to 1024 sample rate (i.e. by a factor of 4):

lowpass_strain = decimate(strain, 4)
lowpass_strain_times = strain_times[::4]
lowpass_dt = dt*4

T = 0.5
seglen = int(round(T/lowpass_dt))
fs, psd = welch(lowpass_strain, fs=1/lowpass_dt, nperseg=seglen)

loglog(fs, psd)
xlabel(r'$f / \mathrm{Hz}$')
ylabel(r'$P(f) / \mathrm{strain}^2 \, \mathrm{Hz}^{-1}$')

Text(0, 0.5, '$P(f) / \\mathrm{strain}^2 \\, \\mathrm{Hz}^{-1}$')
_images/colored_noise_20_1.svg

Let’s whiten the data: the PSD measures the density of the variance in frequency space, so the s.d. in a frequency bin is \(\sqrt{P(f)*\delta f}\):

def whiten(strain, interp_psd, dt):
    Nt = len(strain)
    freqs = np.fft.rfftfreq(Nt, dt)

    # whitening: transform to freq domain, divide by asd, then transform back, 
    # taking care to get normalization right.
    hf = np.fft.rfft(strain)
    white_hf = hf / (np.sqrt(interp_psd(freqs) /dt/2.))
    white_ht = np.fft.irfft(white_hf, n=Nt)
    return white_ht

white_strain = whiten(lowpass_strain, interp1d(fs, psd), lowpass_dt)

plot(lowpass_strain_times, white_strain)

[<matplotlib.lines.Line2D at 0x147bf1c40>]
_images/colored_noise_24_1.svg

Funky things happen at the end of the segment. Why is that? (Periodicity.) It’s better to whiten a windowed version of the data that enforces periodicity. Here we use a Tukey window (which tapers off to zero on either end over one-second intervals):

strain_window = tukey(len(lowpass_strain), alpha=2/32)
plot(lowpass_strain_times, strain_window)

[<matplotlib.lines.Line2D at 0x147c3d940>]
_images/colored_noise_26_1.svg
white_strain = whiten(lowpass_strain*strain_window, interp1d(fs, psd), lowpass_dt)

plot(lowpass_strain_times, white_strain)

[<matplotlib.lines.Line2D at 0x147ca8910>]
_images/colored_noise_28_1.svg

And here’s our merger:

imerg = argmin(abs(lowpass_strain_times - (t_approx - strain_start)))
rng = range(int(round(imerg-0.2/lowpass_dt)),int(round(imerg+0.05/lowpass_dt)))
plot(lowpass_strain_times[rng]-(t_approx-strain_start), white_strain[rng])
xlabel(r'$t - t_c$ $(\mathrm{s})$')
ylabel(r'$h(t)$ (whitened)')

Text(0, 0.5, '$h(t)$ (whitened)')
_images/colored_noise_30_1.svg

A good place to go for discussions like what we are about to have is the “GW150914 basic physics” paper [Abbott et al., 2016]; for “real” post Newtonian expansions, see [Blanchet, 2014]. The basic physics paper gives the evolution of the inspiral phase with time: $\( \frac{1}{2\pi} \frac{\mathrm{d} \Phi}{\mathrm{d} t} = f_\mathrm{GW} = \frac{5^{3/8}}{\left( 8 \pi \right)} \left( \frac{G M_c}{c^3} \right)^{-5/8} \left( t - t_c \right)^{-3/8} \)$

ts = lowpass_strain_times[rng]-(t_approx-strain_start)
fs = (1100.0*(4.9e-6*20.0)**(5.0/3.0)*(-ts))**(-3.0/8.0)
plot(ts, fs)
xlabel(r'$t / \mathrm{s}$')
ylabel(r'$f_{\mathrm{GW}} / \mathrm{Hz}$')

/var/folders/35/vcq24mtj2_g5wk96mbck0cw400018s/T/ipykernel_6604/3102124229.py:2: RuntimeWarning: invalid value encountered in power
  fs = (1100.0*(4.9e-6*20.0)**(5.0/3.0)*(-ts))**(-3.0/8.0)

Text(0, 0.5, '$f_{\\mathrm{GW}} / \\mathrm{Hz}$')
_images/colored_noise_32_2.svg

Because the inspiral frequency reaches \(\infty\) in finite time, we will have to do something. Here’s something:

  • To truncate the frequency at some \(f_\mathrm{max}\) (a parameter).

  • To have an exponentially-damped amplitude decay once the frequency reaches \(f_\mathrm{max}\).

  • To use a power-law amplitude with unknown index.

So $\( h(t) = A_{100} \left( \frac{f(t)}{100 \, \mathrm{Hz}} \right)^\alpha \cos \left( \phi_0 + \Phi(t) \right) \frac{1}{1 + \exp\left( \frac{t - t\left(f_\mathrm{max}\right)}{\tau} \right)} \)\( with \)\( f(t) = \begin{cases} \left( \frac{\left( 8 \pi \right)^{8/3}}{5} \left( \frac{G \mathcal{M}}{c^3} \right)^{5/3} \left( t_c - t \right) \right)^{-3/8} & f(t) < f_\mathrm{max} \\ f_\mathrm{max} & \mathrm{otherwise} \end{cases} \)$

def low_pass_filter_amplitude(fs, frolloff, fwidth):
    return 1.0/(1.0 + exp((fs - frolloff)/fwidth))
def basic_waveform(ts, A100, dlnAdlnf, Mc, tc, phic, max_freq, tau):
    fs = (1100.0*(4.9e-6*Mc)**(5.0/3.0)*(tc-ts))**(-3.0/8.0)

    tmax = tc - max_freq**(-8.0/3.0)/(1100.0*(4.9e-6*Mc)**(5.0/3.0))
    
    Phiinsp = 2.0*pi*8.0/5.0*(1100.0*(4.9e-6*Mc)**(5.0/3.0))**(-3.0/8.0)*(tc - ts)**(5.0/8.0)
    Phiconst = 2.0*pi*8.0/5.0*(1100.0*(4.9e-6*Mc)**(5.0/3.0))**(-3.0/8.0)*(tc - tmax)**(5.0/8.0) + 2.0*pi*max_freq*(tmax-ts)
    
    Phi = where(ts < tmax, Phiinsp, Phiconst)
        
    fs[(fs >= max_freq) | (ts >= tc)] = max_freq
    
    amps = A100*(fs/100)**dlnAdlnf
    
    return (1.0-low_pass_filter_amplitude(fs, 30.0, 2.0))*amps*cos(phic - Phi)*low_pass_filter_amplitude(ts, tmax, tau)

plot(ts, basic_waveform(ts, 2.0, 1.5, 30.0, 0.0, 0.0, 250.0, 0.01))
plot(ts, white_strain[rng])
xlabel(r'$t-t_\mathrm{GW150914}$ ($\mathrm{s}$)')
ylabel(r'Whitened Strain')

/var/folders/35/vcq24mtj2_g5wk96mbck0cw400018s/T/ipykernel_6604/4150586030.py:4: RuntimeWarning: invalid value encountered in power
  fs = (1100.0*(4.9e-6*Mc)**(5.0/3.0)*(tc-ts))**(-3.0/8.0)
/var/folders/35/vcq24mtj2_g5wk96mbck0cw400018s/T/ipykernel_6604/4150586030.py:8: RuntimeWarning: invalid value encountered in power
  Phiinsp = 2.0*pi*8.0/5.0*(1100.0*(4.9e-6*Mc)**(5.0/3.0))**(-3.0/8.0)*(tc - ts)**(5.0/8.0)

Text(0, 0.5, 'Whitened Strain')
_images/colored_noise_35_2.svg
data = {
    'nt': len(ts),
    'ts': ts,
    'h_obs': white_strain[rng]
}
with open('stan/phenom_model.stan', 'r') as f:
    code = f.read()
model = stan.build(program_code=code, data=data)

Building: found in cache, done.
Messages from stanc:
Warning in '/var/folders/35/vcq24mtj2_g5wk96mbck0cw400018s/T/httpstan_zxfqm9oh/model_qsspgd2o.stan', line 138, column 19: Argument 0.01 suggests there may be parameters that are not unit scale; consider rescaling with a multiplier (see manual section 22.12).
Warning in '/var/folders/35/vcq24mtj2_g5wk96mbck0cw400018s/T/httpstan_zxfqm9oh/model_qsspgd2o.stan', line 16, column 29: The variable fs may not have been assigned a value before its use.
Warning: The parameter phi_x has no priors.

fit = az.from_pystan(model.sample(num_chains=4))

Sampling:   0%
Sampling:   0% (1/8000)
Sampling:   0% (2/8000)
Sampling:   0% (3/8000)
Sampling:   0% (4/8000)
Sampling:   1% (103/8000)
Sampling:   3% (203/8000)
Sampling:   4% (302/8000)
Sampling:   5% (401/8000)
Sampling:   6% (501/8000)
Sampling:   8% (600/8000)
Sampling:   9% (700/8000)
Sampling:  10% (800/8000)
Sampling:  11% (900/8000)
Sampling:  12% (1000/8000)
Sampling:  14% (1100/8000)
Sampling:  15% (1200/8000)
Sampling:  16% (1300/8000)
Sampling:  18% (1400/8000)
Sampling:  19% (1500/8000)
Sampling:  20% (1600/8000)
Sampling:  21% (1700/8000)
Sampling:  22% (1800/8000)
Sampling:  24% (1900/8000)
Sampling:  25% (2000/8000)
Sampling:  26% (2100/8000)
Sampling:  28% (2200/8000)
Sampling:  29% (2300/8000)
Sampling:  30% (2400/8000)
Sampling:  31% (2500/8000)
Sampling:  32% (2600/8000)
Sampling:  34% (2700/8000)
Sampling:  35% (2800/8000)
Sampling:  36% (2900/8000)
Sampling:  38% (3000/8000)
Sampling:  39% (3100/8000)
Sampling:  40% (3201/8000)
Sampling:  41% (3301/8000)
Sampling:  43% (3402/8000)
Sampling:  44% (3501/8000)
Sampling:  45% (3601/8000)
Sampling:  46% (3700/8000)
Sampling:  48% (3800/8000)
Sampling:  49% (3901/8000)
Sampling:  50% (4001/8000)
Sampling:  51% (4101/8000)
Sampling:  53% (4201/8000)
Sampling:  54% (4300/8000)
Sampling:  55% (4400/8000)
Sampling:  56% (4500/8000)
Sampling:  58% (4600/8000)
Sampling:  59% (4700/8000)
Sampling:  60% (4800/8000)
Sampling:  61% (4900/8000)
Sampling:  62% (5000/8000)
Sampling:  64% (5100/8000)
Sampling:  65% (5200/8000)
Sampling:  66% (5300/8000)
Sampling:  68% (5400/8000)
Sampling:  69% (5500/8000)
Sampling:  70% (5600/8000)
Sampling:  71% (5700/8000)
Sampling:  72% (5800/8000)
Sampling:  74% (5900/8000)
Sampling:  75% (6000/8000)
Sampling:  76% (6100/8000)
Sampling:  78% (6200/8000)
Sampling:  79% (6301/8000)
Sampling:  80% (6401/8000)
Sampling:  81% (6501/8000)
Sampling:  83% (6601/8000)
Sampling:  84% (6700/8000)
Sampling:  85% (6800/8000)
Sampling:  86% (6900/8000)
Sampling:  88% (7000/8000)
Sampling:  89% (7100/8000)
Sampling:  90% (7200/8000)
Sampling:  91% (7300/8000)
Sampling:  92% (7400/8000)
Sampling:  94% (7500/8000)
Sampling:  95% (7600/8000)
Sampling:  96% (7700/8000)
Sampling:  98% (7800/8000)
Sampling:  99% (7900/8000)
Sampling: 100% (8000/8000)
Sampling: 100% (8000/8000), done.
Messages received during sampling:
  Gradient evaluation took 0.000112 seconds
  1000 transitions using 10 leapfrog steps per transition would take 1.12 seconds.
  Adjust your expectations accordingly!
  Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
  Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/var/folders/35/vcq24mtj2_g5wk96mbck0cw400018s/T/httpstan_60ciys2i/model_qsspgd2o.stan', line 131, column 2 to column 33)
  If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
  but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
  Gradient evaluation took 0.000211 seconds
  1000 transitions using 10 leapfrog steps per transition would take 2.11 seconds.
  Adjust your expectations accordingly!
  Gradient evaluation took 0.000205 seconds
  1000 transitions using 10 leapfrog steps per transition would take 2.05 seconds.
  Adjust your expectations accordingly!
  Gradient evaluation took 0.000205 seconds
  1000 transitions using 10 leapfrog steps per transition would take 2.05 seconds.
  Adjust your expectations accordingly!

az.plot_trace(fit, var_names=['A100', 'dlnAdlnf', 'phi', 'tc', 'Mc', 'max_freq', 'tau', 'sigma'])

array([[<AxesSubplot:title={'center':'A100'}>,
        <AxesSubplot:title={'center':'A100'}>],
       [<AxesSubplot:title={'center':'dlnAdlnf'}>,
        <AxesSubplot:title={'center':'dlnAdlnf'}>],
       [<AxesSubplot:title={'center':'phi'}>,
        <AxesSubplot:title={'center':'phi'}>],
       [<AxesSubplot:title={'center':'tc'}>,
        <AxesSubplot:title={'center':'tc'}>],
       [<AxesSubplot:title={'center':'Mc'}>,
        <AxesSubplot:title={'center':'Mc'}>],
       [<AxesSubplot:title={'center':'max_freq'}>,
        <AxesSubplot:title={'center':'max_freq'}>],
       [<AxesSubplot:title={'center':'tau'}>,
        <AxesSubplot:title={'center':'tau'}>],
       [<AxesSubplot:title={'center':'sigma'}>,
        <AxesSubplot:title={'center':'sigma'}>]], dtype=object)
_images/colored_noise_38_1.svg
plot(ts, fit.posterior.h_model.mean(axis=(0,1)))
plot(ts, white_strain[rng])
xlabel(r'$t-t_\mathrm{GW150914}$ ($\mathrm{s}$)')
ylabel(r'Whitened Strain')

Text(0, 0.5, 'Whitened Strain')
_images/colored_noise_39_1.svg

Note: the “real” value of \(M_c\) is \(\sim 30 \, M_\odot\).

az.plot_posterior(fit, var_names='Mc')

<AxesSubplot:title={'center':'Mc'}>
_images/colored_noise_41_1.svg

State Space Models

N = 1024
xs = linspace(-1, 1, N)

sigma2 = 10.0
tau = 1.0
def k(xi, xj):
    return sigma2*np.exp(-np.abs(xi-xj)/tau)

Sigma = zeros((N,N))
for i in range(N):
    for j in range(N):
        Sigma[i,j] = k(xs[i], xs[j])

imshow(Sigma)

<matplotlib.image.AxesImage at 0x14ae52520>
_images/colored_noise_43_1.svg
for i in range(10):
    plot(xs, np.random.multivariate_normal(zeros_like(xs), Sigma))
_images/colored_noise_44_0.svg

Real term corresponds to $\( \rho(t) = a \exp\left( c t \right) \)$

\[ \frac{\mathrm{d} y}{\mathrm{d} t} + \frac{y}{\tau} = n \]
\[ \left\langle n \right\rangle = 0 \]
\[ \left\langle n(t) n(t') \right\rangle = \sigma^2 \delta(t - t') \]
\[ P_y(f) \sim \frac{\sigma^2 \tau}{1 + \left( 2 \pi f \tau\right)^2} \]
\[ y(t) = y\left( t_0 \right) \exp\left( -\frac{t-t_0}{\tau} \right) + \int_0^t \mathrm{d} s \, \exp\left( -\frac{t-s}{\tau} \right) n(s) \]
\[ \left\langle y(t) y(t') \right\rangle = \frac{\sigma^2}{\tau} \exp\left( -\frac{\left| t - t' \right|}{\tau} \right) \]

Finite Difference analog:

\[ y_{i+1} - a y_i = n_{i+1} \]
a_true = 3
c_true = 3
kernel = terms.RealTerm(a=a_true, c=c_true)
gp = celerite2.GaussianProcess(kernel)

ts = cumsum(rand(1024))/512
plot(ts)

[<matplotlib.lines.Line2D at 0x14ada1b50>]
_images/colored_noise_49_1.svg
gp.compute(ts)

for _ in range(10):
    plot(ts, gp.sample())
_images/colored_noise_51_0.svg

Likelihood: $\( p\left( \left\{ y_i \mid i = 1, \ldots N \right\} \mid a, c \right) = p\left( y_1 \mid a, c \right) \times p\left( y_2 \mid y_1, a, c \right) \times p\left( y_3 \mid y_2, y_1, a, c\right) \ldots \)$

ys = gp.sample()

with pm.Model() as celerite_model:
    a = pm.Lognormal('a', mu=log(1), sigma=1)
    c = pm.Lognormal('c', mu=log(1), sigma=1)
    
    kernel = theano_terms.RealTerm(a=a, c=c)
    gp = celerite2.theano.GaussianProcess(kernel)
    gp.compute(ts)
    _ = gp.marginal('likelihood', observed=ys)

with celerite_model:
    trace = pmx.sample()

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/pymc3_ext/sampling/sampling.py:97: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.
  return pm.sample(draws=draws, tune=tune, model=model, step=step, **kwargs)
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [c, a]
100.00% [4000/4000 00:13<00:00 Sampling 2 chains, 0 divergences] 
Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 29 seconds.
ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

The acceptance probability does not match the target. It is 0.9599042519267924, but should be close to 0.9. Try to increase the number of tuning steps.

with celerite_model:
    pm.plot_trace(trace, lines=[('a', {}, a_true), ('c', {}, c_true)])

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'

ERROR (theano.graph.opt): Optimization failure due to: local_join_empty
ERROR (theano.graph.opt): node: Join(TensorConstant{0}, Elemwise{exp,no_inplace}.0, TensorConstant{[]}, TensorConstant{[]})
ERROR (theano.graph.opt): TRACEBACK:
ERROR (theano.graph.opt): Traceback (most recent call last):
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 2017, in process_node
    replacements = lopt.transform(fgraph, node)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4181, in local_join_empty
    ret = tt.patternbroadcast(ret, node.outputs[0].broadcastable)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/basic.py", line 4181, in patternbroadcast
    return theano.tensor.opt.apply_rebroadcast_opt(rval)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4103, in apply_rebroadcast_opt
    rval2 = local_rebroadcast_lift.transform(None, rval.owner)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/graph/opt.py", line 1209, in transform
    return self.fn(*args, **kwargs)
  File "/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/theano/tensor/opt.py", line 4046, in local_rebroadcast_lift
    if len(fgraph.clients[input]) == 1:
AttributeError: 'NoneType' object has no attribute 'clients'
_images/colored_noise_56_1.svg
\[ \prod_{i=1}^p \left[ \frac{\mathrm{d}}{\mathrm{d} t} - r_i \right] y(t) = \prod_{j=1}^q \left[ \frac{\mathrm{d} }{\mathrm{d} t} - b_i \right] n(t) \]

Recall: $\( q < p \)$ for finite variance.

\[ P_y(f) = \sigma^2 \frac{\prod_{j=1}^q \left| \left(2\pi i f\right) - b_j \right|^2}{\prod_{i=1}^p \left| \left(2 \pi i f \right) - r_i \right|^2} \]
\[ P_y \to \frac{1}{f^{2p - 2q}} \]

as $\( f \to \infty \)$

CARMA(\(p\), \(q\))

\[ r_i = \gamma + 2\pi i f \]
\[ Q = \frac{\pi f}{\gamma} \]
omegas = [2*pi, 2*pi*sqrt(2)]
Qs = [10, 20]
variances = [1.0, 1.0]
kernel = sum([terms.SHOTerm(Q=Qs[i], sigma=sqrt(variances[i]), w0=omegas[i]) for i in range(2)])
gp = celerite2.GaussianProcess(kernel)

dts = 0.5*rand(1024)
ts = cumsum(dts)
gp.compute(ts)
ys = gp.sample()

plot(ts, ys)

[<matplotlib.lines.Line2D at 0x14afb3d00>]
_images/colored_noise_61_1.svg
\[ p\left( y_1 \mid \vec{A}, \vec{Q}, \vec{\omega}_0 \right) \times p\left( y_2 \mid y_1, \vec{A}, \vec{Q}, \vec{\omega}_0 \right) \times p\left( y_3 \mid y_2, y_1, \vec{A}, \vec{Q}, \vec{\omega}_0 \right) \times \ldots \]
\[ N\left[ 0, \Sigma\left(\vec{A}, \vec{Q}, \vec{\omega}_0 \right) \right]\left( y_1 \right) \times N\left[ \mu\left( t_2 \mid y_1 \right), \Sigma\left(t_2 - t_1 \mid \vec{A}, \vec{Q}, \vec{\omega}_0 \right) \right]\left( y_2 \right) \times N\left[ \mu\left( t_3 \mid y_2, y_1 \right), \Sigma\left(t_3 - t_2 \mid \vec{A}, \vec{Q}, \vec{\omega}_0 \right) \right]\left( y_3 \right) \times \ldots \]

Kelly, et al. (2014) or Foreman-Mackey, et al. (2017)

fs = linspace(0, 2, 1024)
plot(fs, kernel.get_psd(2*pi*fs)*2*pi)
yscale('log')
xlabel(r'$f$')
ylabel(r'$P_y(f)$')

Text(0, 0.5, '$P_y(f)$')
_images/colored_noise_63_1.svg
var(ys), sum(variances)

(1.9372381767218398, 2.0)

Nmode = 2
with pm.Model() as sho_model:
    sigmas = pm.Lognormal('sigmas', mu=log(1), sigma=1, shape=Nmode)
    fs = pm.Lognormal('fs', mu=log(1), sigma=1, shape=Nmode)
    Qs = pm.Lognormal('Qs', mu=log(10), sigma=1, shape=Nmode)
    
    gp = celerite2.theano.GaussianProcess(sum([theano_terms.SHOTerm(sigma=sigmas[i], Q=Qs[i], w0=2*pi*fs[i]) for i in range(Nmode)]))
    gp.compute(ts)
    _ = gp.marginal('likelihood', observed=ys)

with sho_model:
    sho_trace = pmx.sample()

Dan Simpson GP Tutorial: https://dansblog.netlify.app/posts/2021-11-03-yes-but-what-is-a-gaussian-process-or-once-twice-three-times-a-definition-or-a-descent-into-madness/

Keplerian and Asteroseismology in RV

Here we will reproduce the results of [Farr et al., 2018]. The data for that paper live here; I have downloaded them and present them here:

dsnames = dict([("table3.dat", "CFHT"), ("table4.dat", "DAO"), ("table5.dat", "McDonald 2.1m"), 
                ("table6.dat", "McDonald 2.7m"), ("table7.dat", "McDonald 2.7m + Tull"),
                ("table8.dat", "TLS"), ("table9.dat", "BOAO"), ("tablesong.dat", "SONG")])

fnames = []
dataframes = []
for fn in glob.glob("data/Aldebaran/*.dat"):
    arr = np.loadtxt(fn)
    fnames.append(op.basename(fn))
    dataframes.append(pd.DataFrame(arr, columns=['time', 'rv', 'rv_uncert']))

Observations of Aldebaran come from several different instruments over several tens of years:

Ns = len(dataframes)
with sns.color_palette('husl', n_colors=Ns):
    for fn, df in zip(fnames, dataframes):
        errorbar(df['time'], df['rv']-mean(df['rv']), yerr=df['rv_uncert'], fmt='.', label=dsnames[fn])
    xlabel(r'BJD (d)')
    ylabel(r'$v_r$ ($\mathrm{m} \, \mathrm{s}^{-1}$)')
    legend(loc='lower center')
_images/colored_noise_73_0.svg

It is customary to give asteroseismic frequencies in micro-Hz; \(1 \, \mu\mathrm{Hz} \simeq 1 / \left( 12 \, \mathrm{d} \right)\).

mu_Hz_day = 1.0 / (3600.0 * 24.0 / 1e6)
mu_Hz_day

11.574074074074073

Here we build the model that we discussed in class. The RV measured at instrument \(i\) at time \(t_j\) is modeled as a sum of a Keplerian RV from a high-mass planet orbiting Aldebaran, a GP with two terms in the kernel—an exponential term for granulation noise on very long timescales, and an SHO term representing an un-resolved superposition of asteroseismic modes whose amplitude peaks at some frequency \(\nu_\mathrm{max}\)—for the stellar surface motion, a site-dependent constant RV offset, and a measurement noise:

\[ v_{r,i}\left( t_j \right) = v_{r,\mathrm{Kepler}}\left( t_j \right) + v_{r, \mathrm{stellar}}\left( t_j \right) + v_{0,i} + \epsilon_{ij} \]

where

\[ \epsilon_{ij} \sim N\left( 0, f_i \sigma_i\left( t_j \right) \right) \]

is the measurement noise term, assumed to be a zero-mean Gaussian with s.d. given by a site-dependent scale factor \(f_i\) times the reported measurement uncertainty for site \(i\) at time \(j\).

Doing a minor bit of algebra, we have

\[ v_{r, \mathrm{stellar}}\left( t_j \right) + \epsilon_{ij} = v_{r,i}\left( t_j \right) - v_{r,\mathrm{Kepler}}\left( t_j \right) - v_{0,i}, \]

which isolates the “stochastic” terms on the left from the “deterministic” terms on the right. celerite2 is capable of handling the combined covariance matrix for the LHS, which consists of the sum of an exponential kernel, an SHO kernel, and a diagonal (a “white” kernel) covariance \(\sim \epsilon^2\).

As usual, we re-center the measurement time series on the median sampled time for better stability in the fit (otherwise we are inferring the time of the pericentre passage closest to Julian day 0 (Monday, 1 January 4713 BC!), which is not a thing that is well-constrained, generating horrible degeneracies with the period—tiny changes in the period couple to tiny shifts in the time of percentre passage observed at Stonehenge near the end of the neolithic era to produce the same RV curve we see today!).

Because our RV measurements are a mix of “intrinsic” and “site-specific” effects, we compute some index gymnastics to keep track of which site is associated to each RV measurement below; then we feed a (sorted) array of times, RVs, measurement uncertainty, and site indices to the model.

Ninst = len(dataframes)

times = np.concatenate([df['time'] for df in dataframes])
indexes = np.concatenate([ind*ones_like(df['time'], dtype=int) for (ind, df) in enumerate(dataframes)])
raw_rvs = np.concatenate([df['rv'] for df in dataframes])
raw_rv_uncerts = np.concatenate([df['rv_uncert'] for df in dataframes])

time_sort_index = argsort(times)

tmid = median(times)

sorted_times = times[time_sort_index] - tmid
sorted_raw_rvs = raw_rvs[time_sort_index]
sorted_raw_rv_uncerts = raw_rv_uncerts[time_sort_index]
sorted_indexes = indexes[time_sort_index]

instrument_mean_rvs = array([mean(df['rv']) for df in dataframes])
instrument_std_rvs = array([std(df['rv']) for df in dataframes])

ts_pred = np.arange(sorted_times[0]-1, sorted_times[-1]+1, mu_Hz_day / 5)

# uniform in log from ~1/100 of nu_max to 100*nu_max
fs = exp(linspace(log(0.01), log(100.0), 1024))

with pm.Model() as aldeb_model:
    # Aosc in m/s is the RMS amplitude of the mode superposition.
    Aosc = pm.Lognormal('Aosc', mu=log(100), sigma=1)
    Qosc = pm.Bound(pm.Lognormal, lower=1, upper=10)('Qosc', mu=log(2), sigma=1)
    # fmax (== nu_max) in muHz
    fmax = pm.Lognormal('fmax', mu=log(2), sigma=0.2)
    
    # Agran in m/s is the RMS amplitude of the granulation signal
    Agran = pm.Lognormal('Agran', mu=log(100), sigma=1)
    # In days!  (For no good reason---just because I was sloppy in class.)
    taugran = pm.Lognormal('taugran', mu=log(30), sigma=1)
    
    # RV semi-amplitude in m/s
    K = pm.Lognormal('K', mu=log(150), sigma=1)
    # Period in days
    P = pm.Normal('P', mu=630, sigma=10)
    
    # Eccentricity vector gives e = |A|, omega = arg(A):
    A = pmx.UnitDisk("A", testval=[0.01, 0.01])
    e = pm.Deterministic("e", T.sqrt(A[0]*A[0] + A[1]*A[1]))
    omega = pm.Deterministic("omega", T.arctan2(A[1], A[0]))

    # The "phase" of t0.  This is a trick: if we choose a 2D vector with N(0, 1) 
    # prior on each component, then its angles will be uniformly distributed, but 
    # with the correct topology: angles crossing from -pi to +pi will be near each 
    # other in this parameter space, while the obvious definition 
    # phi0 = pm.Uniform('phi0', 0, 2*pi) has a discontinuity there.
    phi0_vec = pm.Normal('phi0_vec', mu=0, sigma=1, shape=2, testval=[0.1, 0.1])
    phi0 = pm.Deterministic('phi0', T.arctan2(phi0_vec[1], phi0_vec[0]))
    t0 = pm.Deterministic("t0", P*phi0/(2*pi))
    
    # RV offset at each site, in units of that site's RV s.d., about the mean of the observed RVs.
    dvr0 = pm.Normal('dvr0', mu=0, sigma=1, shape=Ninst)
    # Actual RV offset in m/s
    vr0 = pm.Deterministic('vr0', instrument_mean_rvs + instrument_std_rvs*dvr0)
    
    # RV measurement error scaling factor for each site.
    funcert = pm.Lognormal('funcert', mu=log(1), sigma=2, shape=Ninst)
    
    # Construct the orbit.
    orbit = xo.orbits.KeplerianOrbit(period=P, t0=t0, ecc=e, omega=omega)
    
    # The Kepler part of the RV signal.
    rv_kepler = pm.Deterministic("v_r_kepler", orbit.get_radial_velocity(sorted_times, K=K))

    # Construct the GP with two terms in the kernel.
    gp = celerite2.theano.GaussianProcess(theano_terms.RealTerm(a=Agran*Agran, c=1/taugran) +
                                          theano_terms.SHOTerm(sigma=Aosc, Q=Qosc, w0=2*pi*fmax/mu_Hz_day))
    
    # Factorize the GP covariance matrix, including the diagonal measurement error term (scaled by the site-specific factor)
    gp.compute(sorted_times, yerr=sorted_raw_rv_uncerts*funcert[sorted_indexes])
    
    # This is the "stochastic" part of the RV signal; note the index tricks to apply the constant RV offset
    rv_residual = sorted_raw_rvs - rv_kepler - vr0[sorted_indexes]
    
    rv_stellar_mean = pm.Deterministic('v_r_stellar_mean', gp.predict(rv_residual))
    
    # Here is the likelihood function:
    _ = gp.marginal('stellar_rv_likelihood', observed=rv_residual)
    
    # Now we output a bunch of derived quantities.
    # First the PSD:
    psd = pm.Deterministic('psd', gp.kernel.get_psd(2*pi*fs/mu_Hz_day)*2*pi/mu_Hz_day)
    
    # We also record the keplerian and expected stellar component on a dense grid of times for plotting later
    rv_kepler_pred = pm.Deterministic('v_r_kepler_pred', orbit.get_radial_velocity(ts_pred, K=K))
    rv_stellar_pred_mean = pm.Deterministic('v_r_stellar_pred_mean', gp.predict(rv_residual, t=ts_pred))

with aldeb_model:
    aldeb_trace = pmx.sample()

/Users/wfarr/miniconda3/envs/Astrostatistics688/lib/python3.9/site-packages/pymc3_ext/sampling/sampling.py:97: FutureWarning: In v4.0, pm.sample will return an `arviz.InferenceData` object instead of a `MultiTrace` by default. You can pass return_inferencedata=True or return_inferencedata=False to be safe and silence this warning.
  return pm.sample(draws=draws, tune=tune, model=model, step=step, **kwargs)
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [funcert, dvr0, phi0_vec, A, P, K, taugran, Agran, fmax, Qosc, Aosc]
100.00% [4000/4000 01:57<00:00 Sampling 2 chains, 62 divergences] 
Sampling 2 chains for 1_000 tune and 1_000 draw iterations (2_000 + 2_000 draws total) took 140 seconds.
There were 24 divergences after tuning. Increase `target_accept` or reparameterize.
There were 38 divergences after tuning. Increase `target_accept` or reparameterize.
The number of effective samples is smaller than 25% for some parameters.

with aldeb_model:
    summary = pm.summary(aldeb_trace, var_names=['P', 'dvr0', 'Aosc', 'Qosc', 'fmax', 'Agran', 'K', 'e', 'omega', 't0', 'funcert'])
summary
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
P 628.754 1.883 625.329 632.269 0.057 0.040 1111.0 967.0 1.00
dvr0[0] -0.059 0.161 -0.353 0.243 0.005 0.004 1240.0 754.0 1.00
dvr0[1] -0.066 0.123 -0.291 0.174 0.003 0.003 1305.0 870.0 1.00
dvr0[2] -0.029 0.251 -0.536 0.440 0.007 0.006 1401.0 985.0 1.00
dvr0[3] -0.030 0.156 -0.326 0.245 0.005 0.004 1056.0 680.0 1.00
dvr0[4] -0.860 0.665 -2.078 0.406 0.019 0.013 1299.0 959.0 1.01
dvr0[5] -0.065 0.148 -0.340 0.215 0.004 0.004 1337.0 737.0 1.00
dvr0[6] 0.471 0.200 0.088 0.825 0.006 0.005 1053.0 680.0 1.00
dvr0[7] -0.121 0.230 -0.525 0.326 0.006 0.005 1330.0 739.0 1.00
Aosc 56.166 2.877 50.887 61.653 0.079 0.056 1298.0 970.0 1.00
Qosc 2.970 0.578 1.882 3.984 0.015 0.011 1363.0 1020.0 1.00
fmax 2.272 0.082 2.123 2.426 0.002 0.002 1195.0 1184.0 1.00
Agran 73.348 9.298 57.518 90.331 0.282 0.206 1212.0 1117.0 1.00
K 125.074 13.101 98.087 147.738 0.407 0.288 1033.0 1232.0 1.00
e 0.166 0.070 0.036 0.292 0.002 0.002 1201.0 803.0 1.00
omega -1.367 0.849 -2.919 -0.200 0.028 0.023 742.0 627.0 1.00
t0 146.445 27.314 96.073 200.094 1.121 0.825 608.0 565.0 1.00
funcert[0] 1.760 1.084 0.326 3.936 0.031 0.022 959.0 1113.0 1.00
funcert[1] 0.911 0.654 0.002 2.047 0.020 0.014 811.0 804.0 1.00
funcert[2] 0.821 0.881 0.002 2.540 0.023 0.016 1129.0 1048.0 1.00
funcert[3] 0.031 0.025 0.000 0.074 0.001 0.001 957.0 735.0 1.00
funcert[4] 0.644 0.032 0.587 0.706 0.001 0.001 1283.0 1029.0 1.00
funcert[5] 1.793 1.602 0.003 4.885 0.057 0.040 554.0 185.0 1.00
funcert[6] 0.020 0.015 0.000 0.047 0.001 0.000 1054.0 469.0 1.00
funcert[7] 1.140 0.774 0.004 2.393 0.034 0.024 417.0 364.0 1.00

Here is the traceplot; the model seems to have a lot of “divergences” (places where the ball rolling around the bowl has somehow shot off to “infinity” because it hit a bump or the sampler otherwise failed to accurately integrate the trajectory). These seem (?) to be associated to small funcert (for site 3 and site 6, which are the McDonald 2.7m and 2.1m, the measurement uncertainty scaling factors are \(\lesssim 0.1\), indicating something very strange about those measurements—they are “too good?”)

with aldeb_model:
    pm.plot_trace(aldeb_trace, var_names=['P', 'dvr0', 'Aosc', 'Qosc', 'fmax', 'Agran', 'K', 'e', 'omega', 't0', 'funcert'])
_images/colored_noise_81_0.svg

Normally we would make a residual plot to consider whether the model has correctly “whitened” the residuals. But for performance reasons, we have only output the mean GP prediction for the stellar RV timeseries, which is not the same thing as a fair draw from the GP—it will be more “centrally concentrated,” so the residuals will look too good. Instead, we will show the mean of the means GP (whew!) and the mean Keplerian signal below. Note the “regression to the mean” of the GP between observations—without observations to inform the state space model, the GP mean decays away to zero:

plot(ts_pred, aldeb_trace['v_r_kepler_pred'].mean(axis=0), label='Keplerian')
plot(ts_pred, aldeb_trace['v_r_stellar_pred_mean'].mean(axis=0), label='Stellar GP Mean')
xlabel(r'$t - t_\mathrm{ref}$ ($\mathrm{BJD}$)')
ylabel(r'$v_r$ ($\mathrm{m} \, \mathrm{s}^{-1}$)')
legend(loc='best')

<matplotlib.legend.Legend at 0x14fc75a60>
_images/colored_noise_83_1.svg

Here we plot the Keplerian plus the GP mean against the data for one random draw from the chain (remember—the GP mean will pass “too close” to the data compared to a fair draw from the GP). Every time this cell is executed a different curve will appear.

Nt = aldeb_trace['funcert'].shape[0]
Ns = aldeb_trace['funcert'].shape[1]
with sns.color_palette('husl', n_colors=Ns):
    fig = plt.figure(figsize=(15,5))
    i = randint(Nt)
    for df, v0, fu in zip(dataframes, aldeb_trace['vr0'][i,:], aldeb_trace['funcert'][i,:]):
        errorbar(df['time']-tmid, df['rv']-v0, yerr=df['rv_uncert']*fu, fmt='.')
    plot(ts_pred, aldeb_trace['v_r_stellar_pred_mean'][i,:] + aldeb_trace['v_r_kepler_pred'][i,:], color='k')
    xlabel(r'$t - t_\mathrm{med}$ ($\mathrm{BDJ}$)')
    ylabel(r'$v_r$ ($\mathrm{m} \, \mathrm{s}^{-1}$)')
_images/colored_noise_85_0.svg

Here is the PSD of the stellar GP, clearly showing a peak around \(2 \, \mu\mathrm{Hz}\)

l, = plot(fs, fs*median(aldeb_trace['psd'],axis=0))
fill_between(fs, fs*quantile(aldeb_trace['psd'], 0.84, axis=0), fs*quantile(aldeb_trace['psd'], 0.16, axis=0), color=l.get_color(), alpha=0.25)
fill_between(fs, fs*quantile(aldeb_trace['psd'], 0.975, axis=0), fs*quantile(aldeb_trace['psd'], 0.025, axis=0), color=l.get_color(), alpha=0.25)
xscale('log')
yscale('log')
xlabel(r'$f$ ($\mu\mathrm{Hz}$)')
ylabel(r'$f P(f)$ ($\mathrm{m}^2 \, \mathrm{s}^{-2}$)')

Text(0, 0.5, '$f P(f)$ ($\\mathrm{m}^2 \\, \\mathrm{s}^{-2}$)')
_images/colored_noise_87_1.svg
271-207

64

And the inferred peak value: \(\nu_\mathrm{max} = 2.28^{+0.102}_{-0.064} \, \mu\mathrm{Hz}\).

with aldeb_model:
    pm.plot_posterior(aldeb_trace, var_names='fmax', point_estimate='median', round_to=4, hdi_prob=0.68)
    title('')
    xlabel(r'$\nu_\mathrm{max}$ ($\mu\mathrm{Hz}$)')
_images/colored_noise_90_0.svg

Something very strange is going on with the errorbars in the McDonald instrument—our model would prefer then to be something like 1% of the quoted uncertainty…. We also see that the Song telescope appears to be a bit pessisimistic about its uncertainty (inflated by ~30% relative to what the model would prefer).

with sns.color_palette('husl', n_colors=len(fnames)):
    for i, fn in enumerate(fnames):
        label = dsnames[fn]
        sns.kdeplot(aldeb_trace['funcert'][:,i], log_scale=True, label=label)
    legend(loc='best')
    xlabel(r'$f_\mathrm{uncert}$')
_images/colored_noise_92_0.svg

previous

Regression and Curve Fitting

next

Populations and Catalogs