<jupyter_start><jupyter_text>https://github.com/pymc-devs/pymc3/pull/2876<jupyter_code>%matplotlib inline
import pandas as pd
import numpy as np
import pymc3 as pm
import theano.tensor as tt

import matplotlib.pylab as plt
p = np.random.rand(10, 5)
pm.Categorical.dist(p=p).random()
k = np.shape(p)[-1]
[np.random.choice(k, p=p_, size=None)
                    for p_ in p]<jupyter_output><empty_output><jupyter_text>Cherry picked old commit<jupyter_code>from pymc3.model import modelcontext, get_named_nodes_and_relations
from pymc3.util import is_transformed_name, get_untransformed_name
from tqdm import tqdm


def sample_prior(samples=500, model=None, vars=None, size=None,
                 random_seed=None, progressbar=True):
    """Generate samples from the prior of a model.

    Parameters
    ----------
    samples : int
        Number of samples from the prior to generate. Defaults to 500.
    model : Model (optional if in `with` context)
    vars : iterable
        Variables for which to compute the posterior predictive samples.
        Defaults to `model.named_vars`.
    size : int
        The number of random draws from the distribution specified by the
        parameters in each sample of the trace.
    random_seed : int
        Seed for the random number generator.
    progressbar : bool
        Whether or not to display a progress bar in the command line.

    Returns
    -------
    dict
        Dictionary with the variables as keys. The values are arrays of prior samples.
    """

    model = modelcontext(model)

    if vars is None:
        vars = set(model.named_vars.keys())

    if random_seed is not None:
        np.random.seed(random_seed)

    if progressbar:
        indices = tqdm(range(samples))

    def step(model):
        point = {}
        for var in model.basic_RVs:
            var_name = var.name
            if hasattr(var, 'distribution'):
                if is_transformed_name(var_name):
                    val = var.distribution.dist.random(
                        point=point, size=size)
                    var_name_untranf = get_untransformed_name(var_name)
                    point[var_name_untranf] = val
                    val = var.distribution.transform_used.forward_val(val)
                else:
                    val = var.distribution.random(point=point, size=size)
            else:
                nn, _, _ = get_named_nodes_and_relations(var)
                val = var.eval({model.named_vars[v]: point[v]
                                for v in nn})

            point[var_name] = val
        return point

    try:
        prior = {var: [] for var in vars}
        for _ in indices:
            point = step(model)
            for var_name, var in model.named_vars.items():
                if var_name in vars:
                    if var_name in point:
                        val = point[var_name]
                    else:
                        nn, _, _ = get_named_nodes_and_relations(var)
                        val = var.eval({model.named_vars[v]: point[v]
                                        for v in nn})
                    prior[var_name].append(val)

    except KeyboardInterrupt:
        pass

    finally:
        if progressbar:
            indices.close()

    return {k: np.asarray(v) for k, v in prior.items()}
observed = np.random.normal(10, 1, size=200)
with pm.Model() as m: 
    # Use a prior that's way off to show we're actually sampling from it 
    mu = pm.Normal('mu', mu=-10, sd=1) 
    positive_mu = pm.Deterministic('positive_mu', np.abs(mu)) 
    pm.Normal('x_obs', mu=mu, sd=1, observed=observed) 
    prior = sample_prior() 

assert (prior['mu'] < 0).all() 
assert (prior['positive_mu'] > 0).all() 
assert (prior['x_obs'] < 0).all()<jupyter_output>100%|██████████| 500/500 [00:00<00:00, 1499.66it/s]<jupyter_text>More complex example<jupyter_code>data = pd.read_table(pm.get_data('efron-morris-75-data.tsv'), sep="\t")
at_bats, hits = data[['At-Bats', 'Hits']].values.T
N = len(hits)<jupyter_output><empty_output><jupyter_text>Reference numpy implementation<jupyter_code>import scipy.stats as stats
def model_prior_sample(nsample):
    phi = stats.beta.rvs(1., 1., size=nsample)
    logkappa = np.random.exponential(1/5., nsample)
    kappa = np.exp(logkappa)
    thetas = stats.beta.rvs(phi*kappa, (1.-phi)*kappa, size=(N, nsample))
    y = stats.binom.rvs(n=at_bats[:, None], p=thetas)
    
    psample = {}
    psample['phi'] = phi
    psample['logkappa'] = logkappa
    psample['kappa'] = kappa
    psample['thetas'] = thetas.T
    psample['y'] = y.T
    return psample
%%time

psample0 = model_prior_sample(10000)
plt.style.use('seaborn-darkgrid')
x = np.linspace(0, 1, 200)
a = 1.
b = 1.
pdf = stats.beta.pdf(x, a, b)
plt.plot(x, pdf, label=r'$\alpha$ = {}, $\beta$ = {}'.format(a, b))
plt.hist(psample0['phi'], 500, normed=True);
plt.xlabel('x', fontsize=12)
plt.ylabel('f(x)', fontsize=12)
plt.ylim(0, 4.5)
plt.legend(loc=9)
plt.show()
plt.style.use('seaborn-darkgrid')
x = np.linspace(0, 3, 100)
lam = 5.
pdf = stats.expon.pdf(x, scale=1.0/lam)
plt.plot(x, pdf, label=r'$\lambda$ = {}'.format(lam))
plt.hist(psample0['logkappa'], 500, normed=True);
plt.xlabel('x', fontsize=12)
plt.ylabel('f(x)', fontsize=12)
plt.ylim(0, 4.5)
plt.legend(loc=9)
plt.show()
plt.hist(psample0['y']);
plt.hist(psample0['thetas']);
psample0['thetas'].mean(axis=0) * at_bats
psample0['y'].mean(axis=0)<jupyter_output><empty_output><jupyter_text>Check old implementation<jupyter_code>with pm.Model() as model:
    phi = pm.Beta('phi', alpha=1., beta=1.)

    kappa_log = pm.Exponential('logkappa', lam=5.)
    kappa = pm.Deterministic('kappa', tt.exp(kappa_log))

    thetas = pm.Beta('thetas', alpha=phi*kappa, beta=(1.0-phi)*kappa, shape=N)

    y = pm.Binomial('y', n=at_bats, p=thetas, shape=N, observed=hits)
y.distribution.random(point=model.test_point, size=1000)
np.random.binomial(45, thetas.tag.test_value, size=(1000, N))
# names = ['phi', 'kappa_log', 'kappa', 'thetas', 'y']
# values = pm.distributions.draw_values([model[name] for name in names], size=100)
# p = {k: v for k, v in zip(names, values)}
%%time
psample1 = sample_prior(10000, model=model)
plt.style.use('seaborn-darkgrid')
x = np.linspace(0, 1, 200)
a = 1.
b = 1.
pdf = stats.beta.pdf(x, a, b)
plt.plot(x, pdf, label=r'$\alpha$ = {}, $\beta$ = {}'.format(a, b))
plt.hist(psample1['phi'], 500, normed=True);
plt.xlabel('x', fontsize=12)
plt.ylabel('f(x)', fontsize=12)
plt.ylim(0, 4.5)
plt.legend(loc=9)
plt.show()
np.min(psample1['phi']), np.max(psample1['phi'])
plt.style.use('seaborn-darkgrid')
x = np.linspace(0, 3, 100)
lam = 5.
pdf = stats.expon.pdf(x, scale=1.0/lam)
plt.plot(x, pdf, label=r'$\lambda$ = {}'.format(lam))
plt.hist(psample1['logkappa'], 500, normed=True);
plt.xlabel('x', fontsize=12)
plt.ylabel('f(x)', fontsize=12)
plt.ylim(0, 4.5)
plt.legend(loc=9)
plt.show()
plt.hist(psample1['y']);
plt.hist(psample1['thetas']);
psample1['thetas'].mean(axis=0) * at_bats
psample1['y'].mean(axis=0)<jupyter_output><empty_output><jupyter_text>Check new implementation<jupyter_code>from pymc3.distributions.distribution import draw_values

with pm.Model() as m:
    p = pm.Beta('p', 1., 1.)
draw_values([m['p']], size=10)
import theano

X = theano.shared(np.arange(3))
with pm.Model() as m:
    ind = pm.Categorical('i', np.ones(3)/3)
    x = pm.Deterministic('X', X[ind])
prior = pm.sample_prior_predictive(10, model=m)
prior
with pm.Model() as m:
    p = pm.Beta('p', 1., 1.)
    n = pm.DiscreteUniform('n', 5, 20)
    y = pm.Binomial('y', n, p, observed=5)
    prior = pm.sample_prior_predictive(10)
prior
with pm.Model() as model:
    phi = pm.Beta('phi', alpha=1., beta=1.)

    kappa_log = pm.Exponential('logkappa', lam=5.)
    kappa = pm.Deterministic('kappa', tt.exp(kappa_log))

    thetas = pm.Beta('thetas', alpha=phi*kappa, beta=(1.0-phi)*kappa, shape=N)

    y = pm.Binomial('y', n=at_bats, p=thetas, shape=N, observed=hits)
%%time
psample1 = pm.sample_prior_predictive(10000, model=model)
plt.style.use('seaborn-darkgrid')
x = np.linspace(0, 1, 200)
a = 1.
b = 1.
pdf = stats.beta.pdf(x, a, b)
plt.plot(x, pdf, label=r'$\alpha$ = {}, $\beta$ = {}'.format(a, b))
plt.hist(psample1['phi'], 500, normed=True);
plt.xlabel('x', fontsize=12)
plt.ylabel('f(x)', fontsize=12)
plt.ylim(0, 4.5)
plt.legend(loc=9)
plt.show()
np.min(psample1['phi']), np.max(psample1['phi'])
plt.style.use('seaborn-darkgrid')
x = np.linspace(0, 3, 100)
lam = 5.
pdf = stats.expon.pdf(x, scale=1.0/lam)
plt.plot(x, pdf, label=r'$\lambda$ = {}'.format(lam))
plt.hist(psample1['logkappa'], 500, normed=True);
plt.xlabel('x', fontsize=12)
plt.ylabel('f(x)', fontsize=12)
plt.ylim(0, 4.5)
plt.legend(loc=9)
plt.show()
plt.hist(psample1['y']);
plt.hist(psample1['thetas']);
psample1['thetas'].mean(axis=0) * at_bats
psample1['y'].mean(axis=0)
vars = set(model.named_vars.keys())
names = pm.util.get_default_varnames(model.named_vars, include_transformed=False)
# draw_values fails with auto-transformed variables. transform them later!
print([name for name in names])
values = draw_values([model[name] for name in names], size=1000)
draw_values([model['phi']], size=1000)<jupyter_output><empty_output><jupyter_text>sample from prior with sampler<jupyter_code>with pm.Model() as model_:
    phi = pm.Beta('phi', alpha=1., beta=1.)

    kappa_log = pm.Exponential('logkappa', lam=5.)
    kappa = pm.Deterministic('kappa', tt.exp(kappa_log))

    thetas = pm.Beta('thetas', alpha=phi*kappa, beta=(1.0-phi)*kappa, shape=N)

    y = pm.Binomial('y', n=at_bats, p=thetas, shape=N)
%%time
with model_:
    psample2 = pm.sample(10000, tune=5000, chains=1, step=pm.Metropolis())
plt.hist(psample2['y']);
plt.hist(psample2['thetas']);
psample1['thetas'].mean(axis=0) * at_bats
psample1['y'].mean(axis=0)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1673815359477124, 'fraction_numerical': 0.029513888888888888, 'mean_word_length': 3.492201834862385, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 29, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '541047', 'license': 'MIT', 'n_tokens_mistral': 3957, 'n_tokens_neox': 3632, 'n_words': 759, 'path': 'Miscellaneous/Test_sample_prior.ipynb', 'repo_name': 'junpenglao/Planet_Sakaar_Data_Science', 'size': '9792'}