Skip to content

Instantly share code, notes, and snippets.

@johnnyopao
Last active July 21, 2018 17:25
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save johnnyopao/2c830efb1fe91b28c31905b5939e8e1f to your computer and use it in GitHub Desktop.
Save johnnyopao/2c830efb1fe91b28c31905b5939e8e1f to your computer and use it in GitHub Desktop.
import os
import sys
CWD = os.path.dirname(os.path.realpath(__file__))
sys.path.insert(0, os.path.join(CWD, "lib"))
# The example below was taken and modified from the pyjags documentation
# https://pyjags.readthedocs.io/en/latest/getting_started.html
import pyjags
import numpy as np
np.random.seed(0)
np.set_printoptions(precision=1)
N = 500
a = 70
b = 4
sigma = 50
x = np.random.uniform(0, 100, size=N)
y = np.random.normal(a + x*b, sigma, size=N)
code = '''
model {
for (i in 1:N) {
y[i] ~ dnorm(alpha + beta * x[i], tau)
}
alpha ~ dunif(-1e3, 1e3)
beta ~ dunif(-1e3, 1e3)
tau <- 1 / sigma^2
sigma ~ dgamma(1e-4, 1e-4)
}
'''
model = pyjags.Model(code, data=dict(x=x, y=y, N=N), chains=4)
samples = model.sample(500, vars=['alpha', 'beta', 'sigma'])
def summary(samples, varname, p=95):
values = samples[varname]
ci = np.percentile(values, [100-p, p])
return '{:<6} mean = {:>5.1f}, {}% credible interval [{:>4.1f} {:>4.1f}]'.format(
varname, np.mean(values), p, *ci)
def main(event, context):
print('Did stuff with pyjags')
return {
'Summary': summary(samples, 'alpha')
}
if __name__ == "__main__":
main('', '')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment