<jupyter_start><jupyter_text>Regression and Other Stories: SimulationSimulation of probability models. See Chapter 5 in Regression and Other Stories.<jupyter_code>import arviz as az
from bambi import Model
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats<jupyter_output>WARNING (theano.configdefaults): install mkl with `conda install mkl-service`: No module named 'mkl'
WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.<jupyter_text>Simulate how many girls in 400 births?<jupyter_code>n_girls = stats.binom(400, .488).rvs(1)
print(n_girls)<jupyter_output>[193]<jupyter_text>Repeat simulation 1000 times<jupyter_code># With scipy.stats we do not need a loop to conduct 1000 simulations, but instead can specify it by changing the argument
# to rvs
n_sims = 1000
n_girls = stats.binom(400, .488).rvs(n_sims)<jupyter_output><empty_output><jupyter_text>Plot<jupyter_code>az.plot_dist(n_girls);<jupyter_output><empty_output><jupyter_text>Accounting for twins<jupyter_code>birth_types = pd.Series(["fraternal twin","identical twin","single birth"]) \
.sample(n=400,replace=True, weights=(1/125, 1/300, 1-1/125-1/300))

birth_types.value_counts()<jupyter_output><empty_output><jupyter_text>Plot<jupyter_code># For this smaller dataset we'll use pandas apply. For seriously large datasets consider ufuncs or other optimizations
# outside of the scope of this

def birthtype(birth_type):
    if birth_type == "single birth":
        counts = stats.binom(1, 0.488).rvs(1)
    elif birth_type == "identical twin":
        counts = 2*stats.binom(1, 0.495).rvs(1)
    elif birth_type == "fraternal twin":
        counts = stats.binom(2, 0.495).rvs(1)
    return counts.squeeze()

girls = birth_types.apply(birthtype)
girls.head()
az.plot_dist(n_girls, kind="hist");<jupyter_output><empty_output><jupyter_text>Simulation of continuous and mixed discrete/continuous models<jupyter_code>n_sims = 1000
y1 = stats.norm(3, .5).rvs(n_sims)
y2 = np.exp(y1)
y3 = stats.binom(20, .6).rvs(n_sims)
y4 = stats.poisson(5).rvs(n_sims)<jupyter_output><empty_output><jupyter_text>Plot<jupyter_code>fig, axes = plt.subplots(2,2, figsize=(8,8), constrained_layout=True)
axes[0][0].hist(y1, bins = np.arange(y1.min(), y1.max() + .2, .2))
axes[0][0].set_title("1000 draws from a normal dist. \n with mean 3, sd .5")

axes[0][1].hist(y2, bins = np.arange(0, y2.max() + 5, 5))
axes[0][1].set_title("1000 draws from corresponding \n lognormal dist")


axes[1][0].hist(y3, bins = np.arange(-.5, 20.5, 1))
axes[1][0].set_title("1000 draws from binomial distt with \n 20 tries probability .6")

axes[1][1].hist(y4, bins = np.arange(-.5, y4.max() + 1, 1))
axes[1][1].set_title("1000 draws from the Poisson dist \n with mean 5");<jupyter_output><empty_output><jupyter_text>Generate the height of one randomly chosen adult<jupyter_code>male = stats.binom(1, .48).rvs(1)
height = np.where(male ==1, stats.norm(69.1, 2.9).rvs(1), stats.norm(63.7, 2.7).rvs(1))

male, height<jupyter_output><empty_output><jupyter_text>Select 10 adults at random<jupyter_code>N = 10

male = stats.binom(1, .48).rvs(N)
height = np.where(male==1, stats.norm(69.1, 2.9).rvs(N), stats.norm(63.7, 2.7).rvs(N))
avg_height = height.mean()
avg_height<jupyter_output><empty_output><jupyter_text>Repeat the simulation 1000 times<jupyter_code>n_sims = 1000
avg_height = np.empty(n_sims)

for i in range(n_sims):
    N=10
    male = stats.binom(1, .48).rvs(N)
    height = np.where(male==1, stats.norm(69.1, 2.9).rvs(N), stats.norm(63.7, 2.7).rvs(N))
    avg_height[i] = height.mean()
    
ax = az.plot_dist(avg_height, kind="hist")
ax.set_title("Dist of average height of 10 adults");<jupyter_output><empty_output><jupyter_text>The maximum height of the 10 people<jupyter_code>def height_sim(n):
    male = stats.binom(1, .48).rvs(N)
    height = np.where(male==1, stats.norm(69.1, 2.9).rvs(N), stats.norm(63.7, 2.7).rvs(N))
    max_height = height.max()
    return max_height
max_height = np.array([height_sim(10) for _ in range(1000)])
az.plot_dist(max_height, kind="hist");<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1501093028904542, 'fraction_numerical': 0.05635171241195045, 'mean_word_length': 5.728758169934641, 'pattern_counts': {'":': 3, '<': 47, '<?xml version=': 0, '>': 47, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '55331', 'license': 'MIT', 'n_tokens_mistral': 1763, 'n_tokens_neox': 1599, 'n_words': 343, 'path': 'ROS/ProbabilitySimulation/ProbabilitySimulation.ipynb', 'repo_name': 'farhanreynaldo/Bambi_resources', 'size': '4117'}