<commit_before>import pystan
import numpy as np

# Two compartment model from
# "Stan: A probabilistic programming language for
#  Bayesian inference and optimization" <NAME> (2015)
# http://www.stat.columbia.edu/~gelman/research/published/stan_jebs_2.pdf

a = np.array([0.8, 1.0])
b = np.array([2, 0.1])
sigma = 0.2

x = np.arange(0, 1000, dtype='float')/100
N = len(x)

# The two compartment model we are attempting to fit
y_pred = a[0]*np.exp(-b[0]*x) + a[1]*np.exp(-b[1]*x)

# Include multiplicative noise
y = y_pred * np.exp(np.random.normal(0, sigma, N))

# Compile and sample model
fit = pystan.stan(file='two_compartment.stan',
                  data={'N': N, 'x': x, 'y': y},
                  iter=1000, chains=4)

# Plot parameter estimates of interest
fit.plot(pars=['a', 'b', 'sigma'])

# Print all parameter estimates (limitation of PyStan 2.0)
print(fit)

# Extension: Demonstrate reusing compiled model on new data
y_new = y_pred * np.exp(np.random.normal(0, sigma, N))
fit2 = pystan.stan(fit=fit, data={'N': N, 'x': x, 'y': y_new})
fit2.plot(pars=['a', 'b', 'sigma'])
print(fit2)
<commit_msg>Save and reuse compiled model.  Need to manually delete pkl file to recompile.<commit_after>import pystan
import numpy as np
import os
import pickle


# Two compartment model from
# "Stan: A probabilistic programming language for
#  Bayesian inference and optimization" <NAME> (2015)
# http://www.stat.columbia.edu/~gelman/research/published/stan_jebs_2.pdf

model_file = 'two_compartment.stan'
pkl_file = 'two_compartment.pkl'

a = np.array([0.8, 1.0])
b = np.array([2, 0.1])
sigma = 0.2

x = np.arange(0, 1000, dtype='float')/100
N = len(x)

# The two compartment model we are attempting to fit
y_pred = a[0]*np.exp(-b[0]*x) + a[1]*np.exp(-b[1]*x)

# Include multiplicative noise
y = y_pred * np.exp(np.random.normal(0, sigma, N))

if os.path.isfile(pkl_file):
    # Reuse previously compiled model
    sm = pickle.load(open(pkl_file, 'rb'))
else:
    # Compile and sample model
    sm = pystan.StanModel(file='two_compartment.stan')
    with open(pkl_file, 'wb') as f:
        pickle.dump(sm, f)

fit = sm.sampling(data={'N': N, 'x': x, 'y': y},
                  iter=1000, chains=4)

# Plot parameter estimates of interest
fit.plot(pars=['a', 'b', 'sigma'])

# Print all parameter estimates (limitation of PyStan 2.0)
print(fit)

# Extension: Demonstrate reusing compiled model on new data
y_new = y_pred * np.exp(np.random.normal(0, sigma, N))
fit2 = pystan.stan(fit=fit, data={'N': N, 'x': x, 'y': y_new})
fit2.plot(pars=['a', 'b', 'sigma'])
print(fit2)

{'fraction_non_alphanumeric': 0.16381322957198444, 'fraction_numerical': 0.029571984435797664, 'mean_word_length': 4.152304609218437, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, '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}}
{'dir': 'git-commits-cleaned', 'id': '4967458', 'n_tokens_mistral': 1028, 'n_tokens_neox': 954, 'n_words': 264}