-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
ENH: Use toy example instead of stan MCMC simulation
- Loading branch information
Showing
4 changed files
with
79 additions
and
71 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,27 +1,80 @@ | ||
# Example of Surgical institutional ranking from | ||
# http://www.openbugs.net/Examples/Surgical.html | ||
|
||
#!/usr/bin/env python | ||
"""Model fit example program.""" | ||
import argparse | ||
import collections | ||
import itertools | ||
import os | ||
import time | ||
|
||
|
||
RESULTS_DIR = 'results' | ||
PARAMETERS = collections.OrderedDict(( | ||
('x', [0.1 * i for i in range(1, 11)]), | ||
('y', [-0.1, 0.1]), | ||
('z', ['control', 'positive', 'negative']))) | ||
COMBINATIONS = list(itertools.product(*PARAMETERS.values())) | ||
|
||
|
||
def n_remaining(): | ||
"""Return number of completed combinations.""" | ||
total = len(COMBINATIONS) | ||
if not os.path.exists(RESULTS_DIR): | ||
completed = 0 | ||
else: | ||
completed = len([name | ||
for name in os.listdir(RESULTS_DIR) | ||
if os.path.isfile(os.path.join(RESULTS_DIR, name))]) | ||
return total - completed | ||
|
||
|
||
def run(index): | ||
"""Run example model fit.""" | ||
prefix = '{}: '.format(index) | ||
record = collections.OrderedDict( | ||
zip(PARAMETERS.keys(), | ||
COMBINATIONS[index - 1])) | ||
print(prefix + 'Fitting model to parameters: {} ...'.format( | ||
", ".join(["{} = {}".format(param, value) | ||
for param, value in record.items()]))) | ||
value = {'positive': 0.8, 'negative': 0.1, 'control': 0.3} | ||
result = record['x'] ** record['y'] / value[record['z']] | ||
time.sleep(1) | ||
path_result = os.path.join(RESULTS_DIR, '{:02d}.dat'.format(index)) | ||
if not os.path.exists(RESULTS_DIR): | ||
os.mkdir(RESULTS_DIR) | ||
with open(path_result, 'w') as handle: | ||
handle.writelines([str(result)]) | ||
print(prefix + '... done! Saved result {:6.3f} to {}'.format(result, path_result)) | ||
|
||
|
||
def parse_args(): | ||
"""Parse command-line arguments.""" | ||
parser = argparse.ArgumentParser( | ||
description='Model fit example program.') | ||
parser.add_argument('--sim', action='store_true', | ||
help='Number of simulations to run') | ||
parser.add_argument('--remaining', action='store_true', | ||
help='Remaining simulations to run') | ||
parser.add_argument('indices', metavar='I', nargs='*', type=int, | ||
help='Model indices to run') | ||
return parser | ||
|
||
|
||
def main(argv=None): | ||
"""Entry point for model.py""" | ||
parser = parse_args() | ||
args = parser.parse_args(argv) | ||
if args.sim: | ||
print(len(COMBINATIONS)) | ||
return | ||
if args.remaining: | ||
print(n_remaining()) | ||
return | ||
if not args.indices: | ||
parser.error('Need at least 1 index to run!') | ||
for index in args.indices: | ||
run(index) | ||
|
||
|
||
import pystan | ||
|
||
# This example considers mortality rates in 12 hospitals performing | ||
# cardiac surgery in babies. The data are shown below. | ||
# | ||
DATA = { | ||
# Number of hospitals. | ||
'N': 12, | ||
# Number of operations. | ||
'n': [ 47, 148, 119, 810, 211, 196, 148, 215, 207, 97, 256, 360], # noqa: E201, E501 | ||
# Number of deaths. | ||
'r': [ 0, 18, 8, 46, 8, 13, 9, 31, 14, 8, 29, 24] # noqa: E201, E501 | ||
} | ||
|
||
if __name__ == "__main__": | ||
# Disable threading. | ||
os.environ['STAN_NUM_THREADS'] = "1" | ||
|
||
fixed_effects = pystan.StanModel(file="model_indep.stan") | ||
random_effects = pystan.StanModel(file="model_dep.stan") | ||
burn_in = fixed_effects.sampling(data=DATA, warmup=1000, iter=0) | ||
fit = random_effects.sampling(fit=burn_in, data=DATA, warmup=0, iter=10000) | ||
if __name__ == '__main__': | ||
main() |
This file was deleted.
Oops, something went wrong.
This file was deleted.
Oops, something went wrong.