<jupyter_start><jupyter_text>now fit an AR(1) model to the timeseries. We do not assume the time series is centred on zero (because in production you won't know the centre). Hence the timeseries is shifted to zero by the 'center' parameter. Fitting the AR(1) model has the consequence of also finding where the timeseries is centred, which gives us an HPD for the timeseries mean.<jupyter_code>with pm.Model() as ar1:
    rho_ = pm.Uniform('rho',-1,1) #we assume process is stationary, so -1<k_<1
    tau_ = pm.Gamma('tau',mu=1,sd=1)
    #center = pm.Normal('center', mu=timeseries.mean(), sigma=5) #set the prior for the true mean initially centred on the population mean
    likelihood = pm.AR1('likelihood', k=rho_, tau_e=tau_, observed=timeseries)#-center)
    trace = pm.sample(target_accept=0.9)
    
kdeplot(trace['tau'])
l =ar1['likelihood']
pm.summary(trace)
np.sqrt(1 / ((1-0.99**2)*1))
import theano.tensor as tt
tt.sqrt(1 / ((1-0.99**2)*1))<jupyter_output>INFO (theano.gof.compilelock): Refreshing lock /Users/lmar3213/.theano/compiledir_Darwin-17.7.0-x86_64-i386-64bit-i386-3.7.6-64/lock_dir/lock<jupyter_text>my own version of ar1() using correct variance:<jupyter_code>import theano.tensor as tt

class fixedAR1(pm.distributions.Continuous):
    """
    Autoregressive process with 1 lag.
    Parameters
    ----------
    k: tensor
       effect of lagged value on current value
    tau_e: tensor
       precision for innovations
    """

    def __init__(self, k, tau_e, *args, **kwargs):
        super().__init__(*args, **kwargs)
        self.k = k = tt.as_tensor_variable(k)
        self.tau_e = tau_e = tt.as_tensor_variable(tau_e)
        self.tau = tau_e * (1 - k ** 2)
        self.mode = tt.as_tensor_variable(0.)

    def logp(self, x):
        """
        Calculate log-probability of AR1 distribution at specified value.
        Parameters
        ----------
        x: numeric
            Value for which log-probability is calculated.
        Returns
        -------
        TensorVariable
        """
        k = self.k
        tau_e = self.tau_e

        x_im1 = x[:-1]
        x_i = x[1:]
        
        ##removed this:
        #boundary = Normal.dist(0., tau=tau_e).logp
        ##added this:
        var_ar1 = 1 / ((1-k**2)*tau_e)
        sd_ar1 = tt.sqrt(var_ar1)
        boundary = pm.Normal.dist(0., sigma=sd_ar1).logp

        innov_like = pm.Normal.dist(k * x_im1, tau=tau_e).logp(x_i)
        return boundary(x[0]) + tt.sum(innov_like)

np.random.seed(3) 
def gen_correlated_curve(rho, num=100):
    num_with_runup = num+5000
    y = np.zeros((num_with_runup,))
    for i in range(1,num_with_runup):
        y[i] = rho * y[i-1] + np.random.normal()
    return y[-num:]

timeseries = gen_correlated_curve(0.8, 30)
print(timeseries.mean())
plt.plot(timeseries)
with pm.Model() as ar1:
    rho_ = pm.Uniform('rho',-1,1) #we assume process is stationary, so -1<k_<1
    tau_ = pm.Gamma('tau',mu=1,sd=1)
    center = pm.Normal('center', mu=timeseries.mean(), sigma=5) #set the prior for the true mean initially centred on the population mean
    likelihood = fixedAR1('likelihood', k=rho_, tau_e=tau_, observed=timeseries-center)
    trace = pm.sample(target_accept=0.9)
    
print(pm.stats.hpd(trace['center']))

from statsmodels.tsa.ar_model import AutoReg
fit_result = AutoReg(timeseries-timeseries.mean(), lags = [1]).fit()
estimated_rho = fit_result.params[1]
estimated_rho
pm.stats.hpd(trace['center'])
pm.summary(trace)<jupyter_output><empty_output><jupyter_text>Trying a custom AR(1):See: https://stats.stackexchange.com/questions/346747/pymc3-using-ar1-process-as-prior-causes-bad-initial-energy-error<jupyter_code>from scipy import stats
np.random.seed(2) 
def simple_ar(n_samples, k, sig, mu):
    n = n_samples+5000
    out = np.zeros((n,), dtype='float64')
    out[0] = mu

    for i in range(1, n):
        out[i] = k * (out[i-1] - mu) + stats.norm.rvs(scale=sig, size=(1,)) + mu

    return out[-n_samples:]

out = simple_ar(3000, 0.99, 1,0)
plt.plot(out)
kdeplot(stats.norm.rvs(scale=1, size=(1000,)))
kdeplot([np.random.normal() for _ in range(1000)])
np.random.seed(2) 
def simple_ar(n_samples, k, sig, mu):
    n = n_samples+5000
    out = np.zeros((n,), dtype='float64')
    out[0] = mu

    for i in range(1, n):
        out[i] = k * (out[i-1] - mu) + np.random.normal() + mu

    return out[-n_samples:]

out = simple_ar(3000, 0.99, 0, 0)
print(np.mean(out))
import theano as T
import theano.tensor as tt
import pymc3 as pm
from pymc3.distributions.distribution import generate_samples, draw_values

class NonzeroMeanAR1(pm.distributions.continuous.PositiveContinuous):
    """
    Autoregressive process with 1 lag.
    Parameters
    ----------
    k : tensor
       effect of lagged value on current value
    tau_e : tensor
       precision for innovations
    """

    def __init__(self, k, tau_e, mu, *args, **kwargs):
        super(NonzeroMeanAR1, self).__init__(*args, **kwargs)
        self.k = k = tt.as_tensor_variable(k)
        self.tau_e = tau_e = tt.as_tensor_variable(tau_e)
        self.tau = tau_e * (1 - k ** 2)
        self.mu = tt.as_tensor_variable(mu)
        #self.mode = tt.as_tensor_variable(0.)
        self.mean = self.mu

    def logp(self, x):
        k = self.k
        tau_e = self.tau_e
        mu = self.mu

        x_im1 = x[:-1]
        x_i = x[1:]
        boundary = pm.Normal.dist(mu=mu, tau=tau_e).logp(x[0])

        innov_like = pm.Normal.dist(mu=k * (x_im1 - mu) + mu, tau=tau_e).logp(x_i)
        return boundary + tt.sum(innov_like)

    def random(self, point=None, size=None):
        tau = draw_values([self.tau], point=point)[0]
        mu = draw_values([self.mu], point=point)[0]
        k = draw_values([self.k], point=point)[0]
        eps = generate_samples(sps.norm.rvs, loc=0, scale=tau**-0.5,
                               dist_shape=self.shape,
                               size=size)

        y, _ = T.scan(lambda err, prev, k, mu: k*(prev - mu) + mu + err,
                      outputs_info=mu*tt.ones((1,), dtype='float64'),
                      sequences=[eps],
                      non_sequences=[k, mu])

        return y
import scipy.stats as sps
def simple_ar(n_samples, k, sig, mu):
    out = np.zeros((n_samples,), dtype='float64')
    out[0] = mu

    for i in range(1, n_samples):
        out[i] = k * (out[i-1] - mu) + sps.norm.rvs(scale=sig, size=(1,)) + mu

    return out

y_obs = simple_ar(500, 0.9, 0.2, 100)

with pm.Model() as ar1:
    rho_ = pm.Uniform('rho',-1,1) #we assume process is stationary, so -1<k_<1
    tau_ = pm.Gamma('tau',mu=1,sd=1)
    center = pm.Normal('center', mu=y_obs.mean(), sigma=5) #set the prior for the true mean initially centred on the population mean
    likelihood = NonzeroMeanAR1('likelihood', k=rho_, tau_e=tau_, mu = center, observed=y_obs)
    trace = pm.sample(target_accept=0.9)
    
pm.stats.hpd(trace['center'])
with pm.Model() as ar_model:
    k = pm.HalfFlat('k')
    tau = pm.HalfNormal('tau', sd=100)
    mu = pm.HalfFlat('mu')
    y = pm.AR1('y', k=k, tau_e=tau, observed=y_obs-mu)
    trace = pm.sample(1000, tune=500)
    
print(pm.stats.hpd(trace['center']))

ts = timeseries+10
with pm.Model() as ar1:
    rho_ = pm.Uniform('rho',-1,1) #we assume process is stationary, so -1<k_<1
    tau_ = pm.Gamma('tau',mu=1,sd=1)
    center = pm.Normal('center', mu=ts.mean(), sigma=5) #set the prior for the true mean initially centred on the population mean
    likelihood = pm.AR1('likelihood', k=rho_, tau_e=tau_, observed=ts-center)
    trace = pm.sample(target_accept=0.9)
    
print(pm.stats.hpd(trace['center']))
kdeplot(trace['center'])<jupyter_output>INFO (theano.gof.compilelock): Refreshing lock /Users/lmar3213/.theano/compiledir_Darwin-17.7.0-x86_64-i386-64bit-i386-3.7.6-64/lock_dir/lock
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [center, tau, rho]
Sampling 2 chains, 0 divergences: 100%|██████████| 2000/2000 [00:01<00:00, 1173.06draws/s]<jupyter_text>Analysis:<jupyter_code>#HPD is way off:
pm.stats.hpd(trace['center'])
#Trajectory looks OK:
plt.plot(trace['center'])
#it's definitely got lost:
kdeplot(trace['center'])<jupyter_output><empty_output><jupyter_text>Comparison to other measures:<jupyter_code>#Now we can compare to another technique, based on the autocorrelation function. This is from 
#Use of the Weighted Histogram Analysis Method for the Analysis of Simulated and Parallel Tempering Simulations
#First it calculates the autocorrelation time constant, then calculates SEM based on 'effective number of samples'
#then just apply the 95 CI to that SEM. 
z95 = 1.959963984540054

#this is the Chodera lab 'statistical inefficiency' measure.
def statistical_inefficiency(corr, mintime=3):
    N = corr.size
    C_t = sm.tsa.stattools.acf(corr, fft=True, unbiased=True, nlags=N)
    t_grid = np.arange(N).astype('float')
    g_t = 2.0 * C_t * (1.0 - t_grid / float(N))
    ind = np.where((C_t <= 0) & (t_grid > mintime))[0][0]
    g = 1.0 + g_t[1:ind].sum()
    return max(1.0, g)

 #and these actually return the values for SEM:
def ci_from_chodera(timeseries):
    autocorrelation_time = statistical_inefficiency(timeseries)
    #original n
    n = len(timeseries)
    #calculate SEM and return CI using effective n:
    sem = np.std(timeseries) / np.sqrt(n/autocorrelation_time)
    #return CI:
    return timeseries.mean()-sem*z95, timeseries.mean()+sem*z95

print(sem_from_chodera(timeseries))<jupyter_output>(-5.056946809414023, 1.8577243818480016)
{'fraction_non_alphanumeric': 0.14499738082765845, 'fraction_numerical': 0.03719224724986904, 'mean_word_length': 3.248331108144192, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 19, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '547458', 'license': 'MIT', 'n_tokens_mistral': 3732, 'n_tokens_neox': 3477, 'n_words': 814, 'path': 'scrapbook_code/autocorr_for_discourse.ipynb', 'repo_name': 'ljmartin/estimating_sem_timeseries', 'size': '9545'}