<jupyter_start><jupyter_text>_Bayesian Framework & Probability Distribution_<jupyter_code>%matplotlib inline
from IPython.core.pylabtools import figsize
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
import scipy.stats as stats
# Example : coin flip
figsize(11, 9)
dist = stats.beta
n_trials = [0, 1, 2, 3, 4, 5, 8, 15, 50, 500]
data = stats.bernoulli.rvs(0.5, size=n_trials[-1])
x = np.linspace(0, 1, 100)

for k, N in enumerate(n_trials):
    sx = plt.subplot(len(n_trials) / 2, 2, k + 1)
    plt.setp(sx.get_yticklabels(), visible=False)
    heads = data[:N].sum()
    y = dist.pdf(x, 1 + heads, 1 + N - heads)
    plt.plot(x, y, label="%d flip count, \n%d head count" % (N, heads))
    
    plt.fill_between(x, 0, y, color="#348ABD", alpha=0.4)
    plt.vlines(0.5, 0, 4, color="k", linestyles="--", lw=1)
    
    leg = plt.legend()
    leg.get_frame().set_alpha(0.4)
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>_Poisson-distribution_<jupyter_code>figsize(12.5, 4)
a = np.arange(16)
poi = stats.poisson
lambda_ = [1.5, 4.25]
colours = ["#348ABD", "#A60628"]

plt.subplot(1, 3, 1)
plt.bar(a, poi.pmf(a, lambda_[0]), color=colours[0])

plt.subplot(1, 3, 2)
plt.bar(a, poi.pmf(a, lambda_[1]), color=colours[1])

plt.subplot(1, 3, 3)
plt.bar(a, poi.pmf(a, lambda_[0]), color=colours[0], label="$\lambda = %.1f$" % lambda_[0], \
        alpha=0.60)
plt.bar(a, poi.pmf(a, lambda_[1]), color=colours[1], label="$\lambda = %.1f$" % lambda_[0], \
        alpha=0.60)

plt.xticks(a + 0.4, a)
plt.legend()<jupyter_output><empty_output><jupyter_text>_exponential distribution_<jupyter_code>a = np.linspace(0, 4, 100)
expo = stats.expon
lambda_ = [0.5, 1]
for l, c in zip(lambda_, colours):
    plt.plot(a, expo.pdf(a, scale=1. / l))<jupyter_output><empty_output><jupyter_text>_Bayesian Inference Using Computer_- Example : Action Inference using text messages.<jupyter_code># 데이터 가져오기
# 변수 선언

figsize(12.5, 3.5)

count_data = np.loadtxt(r"01_txtdata.csv")
n_count_data = len(count_data)
n_count_data, count_data
plt.bar(np.arange(n_count_data), count_data)
plt.xlabel("Time(days)", fontsize=13)
plt.ylabel("Numbers of messages")
plt.xlim(0, n_count_data)
# 모수의 분포 구하기

import pymc as pm
# pymc 변수 선언

alpha = 1.0/count_data.mean()

lambda_1 = pm.Exponential("lambda_1", alpha)
lambda_2 = pm.Exponential("lambda_2", alpha)

tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data)
@pm.deterministic
def lambda_(tau=tau, lambda_1=lambda_1, lambda_2=lambda_2):
    out = np.zeros(n_count_data)
    out[:tau] = lambda_1
    out[tau:] = lambda_2
    return out
observation = pm.Poisson("obs", lambda_, value=count_data, observed=True)

model = pm.Model([observation, lambda_1, lambda_2, tau])
mcmc = pm.MCMC(model)
mcmc.sample(40000, 10000)
np.mean(lambda_1_samples)
lambda_1_samples = mcmc.trace('lambda_1')[:]
lambda_2_samples = mcmc.trace('lambda_2')[:]
tau_samples = mcmc.trace('tau')[:]
lambda_1_samples == lambda_2_samples
np.ones_like(tau_samples)

# w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples)
# w.shape
figsize(14.5, 10)

ax = plt.subplot(311)
ax.set_autoscaley_on(False)

plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85, \
         label="$\lambda_1$'s posterior pd'", color="#A60628", normed=True)
plt.xlim([10, 30])

ax = plt.subplot(312)
ax.set_autoscaley_on(False)
plt.hist(lambda_2_samples, histtype='stepfilled', bins=30, alpha=0.85, \
         label="$\lambda_2$'s posterior pd'", color="#7A68A6", normed=True,)
plt.xlim([10, 30])

plt.subplot(313)
w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples) # scaling
plt.hist(tau_samples, bins=n_count_data, alpha=1, color="#467821", weights=w, rwidth=2.)
plt.xticks(np.arange(n_count_data))
plt.ylim([0, .75])
plt.xlim([35, len(count_data) - 20])
tau_samples
# 해석

figsize(12.5, 5)

N = tau_samples.shape[0]
expected_texts_per_day = np.zeros(n_count_data)

for day in range(0, n_count_data):
    ix = day < tau_samples
    expected_texts_per_day[day] = (lambda_1_samples[ix].sum() + lambda_2_samples[~ix].sum()) / N
    
plt.plot(range(n_count_data), expected_texts_per_day, lw=4, color="#E24A33", \
         label="message's expectation")
plt.bar(np.arange(len(count_data)), count_data, color="#348ABD", alpha=0.65)<jupyter_output><empty_output><jupyter_text>*부록*- 두 람다가 정말 다른지 통계적으로 알 수 있을까?- 변환점 두 개로 확장하기<jupyter_code># 두 람다가 정말 다른지 통계적으로 알 수 있을까?
print((lambda_1_samples < lambda_2_samples).sum())
print(lambda_1_samples.shape[0], lambda_2_samples.shape[0])
print((lambda_1_samples < lambda_2_samples).mean())
for d in [1,2,3,4,5,6,7,8,9,10]:
    v = (abs(lambda_1_samples - lambda_2_samples) >= d).mean()
    print("Larger than %d? %.2f" % (d,v))
# 변환점 두 개로 확장하기
lambda_1 = pm.Exponential("lambda_1", alpha)
lambda_2 = pm.Exponential("lambda_2", alpha)
lambda_3 = pm.Exponential("lambda_3", alpha)

tau_1 = pm.DiscreteUniform("tau_1", lower=0, upper=n_count_data-1)
tau_2 = pm.DiscreteUniform("tau_2", lower=tau_1, upper=n_count_data)

@pm.deterministic
def lambda_(tau_1=tau_1, tau_2=tau_2, lambda_1=lambda_1, lambda_2=lambda_2, lambda_3=lambda_3):
    out = np.zeros(n_count_data)
    out[:tau_1] = lambda_1
    out[tau_1:tau_2] = lambda_2
    out[tau_2:] = lambda_3
    
    return out
observation = pm.Poisson("obs", lambda_, value=count_data, observed=True)
model = pm.Model([observation, lambda_1, lambda_2, lambda_3, tau_1, tau_2])
mcmc = pm.MCMC(model)
mcmc.sample(40000, 10000)
lambda_1_samples = mcmc.trace('lambda_1')[:]
lambda_2_samples = mcmc.trace('lambda_2')[:]
lambda_3_samples = mcmc.trace('lambda_3')[:]
tau_1_samples = mcmc.trace('tau_1')[:]
tau_2_samples = mcmc.trace('tau_2')[:]
figsize(14.5, 10)

ax = plt.subplot(511)
ax.set_autoscaley_on(False)
plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85, \
         label="$\lambda_1$'s posterior pd'", color="#A60628", normed=True)

ax = plt.subplot(512)
ax.set_autoscaley_on(False)
plt.hist(lambda_2_samples, histtype='stepfilled', bins=30, alpha=0.85, \
         label="$\lambda_2$'s posterior pd'", color="#7A68A6", normed=True,)

ax = plt.subplot(513)
ax.set_autoscaley_on(False)
plt.hist(lambda_3_samples, histtype='stepfilled', bins=30, alpha=0.85, \
         label="$\lambda_3$'s posterior pd'", color="#7A68A6", normed=True,)

plt.subplot(514)
w = 1.0 / tau_1_samples.shape[0] * np.ones_like(tau_1_samples)
plt.hist(tau_1_samples, bins=n_count_data, alpha=1, color="#467821", weights=w, rwidth=2.)
plt.xticks(np.arange(n_count_data))
plt.ylim([0, .75])
plt.xlim([35, len(count_data) - 20])

plt.subplot(515)
w = 1.0 / tau_2_samples.shape[0] * np.ones_like(tau_2_samples)
plt.hist(tau_2_samples, bins=n_count_data, alpha=1, color="#467821", weights=w, rwidth=2.)
plt.xticks(np.arange(n_count_data))
plt.ylim([0, .75])
plt.xlim([35, len(count_data) - 20])<jupyter_output>/Users/quartz/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.
  warnings.warn("The 'normed' kwarg is deprecated, and has been "
/Users/quartz/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.
  warnings.warn("The 'normed' kwarg is deprecated, and has been "
/Users/quartz/anaconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py:6462: UserWarning: The 'normed' kwarg is deprecated, and has been replaced by the 'density' kwarg.
  warnings.warn("The 'normed' kwarg is deprecated, and has been "
{'fraction_non_alphanumeric': 0.16715347849176845, 'fraction_numerical': 0.05589484864577801, 'mean_word_length': 5.505181347150259, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 21, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '363029', 'license': 'MIT', 'n_tokens_mistral': 3435, 'n_tokens_neox': 3261, 'n_words': 577, 'path': '01_bayesian_interference.ipynb', 'repo_name': 'seokyeongheo/bayesian-estimation-study', 'size': '7532'}