From 2fa11d48a8d1215b143a9e401820ab153fecd1fe Mon Sep 17 00:00:00 2001 From: Pariksheet Nanda Date: Mon, 13 May 2019 18:22:29 -0400 Subject: [PATCH] ENH: Use toy example instead of stan MCMC simulation --- README.md | 4 +- examples/model.py | 103 +++++++++++++++++++++++++++++--------- examples/model_dep.stan | 27 ---------- examples/model_indep.stan | 16 ------ 4 files changed, 79 insertions(+), 71 deletions(-) mode change 100644 => 100755 examples/model.py delete mode 100644 examples/model_dep.stan delete mode 100644 examples/model_indep.stan diff --git a/README.md b/README.md index 3bc3fd4..3f81348 100644 --- a/README.md +++ b/README.md @@ -139,10 +139,8 @@ task workers complete records and become available. ```sh # From the command-line -python2.7 -m pip install --user --upgrade pip -python2.7 -m pip install --user --upgrade pystan module purge -module load gcc/4.9.3 # Needed for pystan 2.19.0 +module load python/2.7.6-gcc-unicode cd ~/parallel-slurm/examples rm -rf joblog submit.out results/ results.csv for i in {1..5}; do sbatch 03-submit-param-sweep.slurm; done diff --git a/examples/model.py b/examples/model.py old mode 100644 new mode 100755 index 1413db4..edd1691 --- a/examples/model.py +++ b/examples/model.py @@ -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() diff --git a/examples/model_dep.stan b/examples/model_dep.stan deleted file mode 100644 index d953256..0000000 --- a/examples/model_dep.stan +++ /dev/null @@ -1,27 +0,0 @@ -data { - int N; - int n[N]; - int r[N]; -} - -parameters { - real mu; - real tau; - real b[N]; - real p[N]; - real pop_mean; - real sigma; -} - -model { - for (i in 1:N) { - b[i] ~ normal_lpdf(mu, tau); - r[i] ~ binomial_lpmf(n[i], p[i]); - p[i] ~ logit(b[i]); - } - - pop_mean = exp(mu) / (1 + exp(mu)); - mu ~ normal(0.0, 1.0E-6); - sigma = 1 / sqrt(tau); - tau ~ gamma(0.001, 0.001); -} diff --git a/examples/model_indep.stan b/examples/model_indep.stan deleted file mode 100644 index 3031202..0000000 --- a/examples/model_indep.stan +++ /dev/null @@ -1,16 +0,0 @@ -data { - int N; - int n[N]; - int r[N]; -} - -parameters { - real p[N]; -} - -model { - for (i in 1:N) { - p[i] ~ beta_lpdf(1.0, 1.0); - r[i] ~ binomial_lpmf(n[i], p[i]); - } -}