Stan Coin flipping Test

model.stan

data {
  int<lower=0> N;
}

transformed data {
  int<lower=0> n;
  n = N / 2;
}

parameters {
  real<lower=0,upper=1> p;
}

model {
  N ~ binomial(N * 2, p);
}

run.py

import pickle
from hashlib import md5
from pathlib import Path

import numpy as np
from pystan import StanModel


def StanModel_cache(model_code, model_name="anon_model", model_dir="~/.stan"):
    model_dir = Path(model_dir).expanduser()
    if not model_dir.exists():
        model_dir.mkdir()

    code_hash = md5(model_code.encode('ascii')).hexdigest()
    if model_name is None:
        cache_fn = 'cached-model-{}.pkl'.format(code_hash)
    else:
        cache_fn = 'cached-{}-{}.pkl'.format(model_name, code_hash)

    cache_file: Path = model_dir.joinpath(cache_fn)
    if cache_file.exists():
        print("use cached stan model")
        with cache_file.open("rb") as fp:
            sm = pickle.load(fp)
    else:
        print("compile stan model")
        sm = StanModel(model_code=model_code, model_name=model_name)
        with cache_file.open(mode="wb") as fp:
            pickle.dump(sm, fp, pickle.HIGHEST_PROTOCOL)
    return sm


def main():
    model = StanModel_cache(open("model.stan").read())
    for N in (10, 50, 100, 200, 300, 500, 1000):
        res = model.sampling(data={"N": N})
        p = res["p"]
        s = np.sum((0.45 < p) & (p < 0.55)) / len(p)
        print(N, s)
        if s > 0.95:
            break


if __name__ == "__main__":
    main()