書記が数学やるだけ#783 GLMモデルのベイズ統計

import pymc3 as pm
import math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import load_breast_cancer
%matplotlib inline

data_breast_cancer = load_breast_cancer()
df_target = pd.DataFrame(data_breast_cancer["target"], columns=["target"])
df_data = pd.DataFrame(data_breast_cancer["data"], columns=data_breast_cancer["feature_names"])
df = pd.concat([df_target, df_data], axis=1)
df1 = df[["mean radius","target"]]
df1 = df1.rename(columns={'mean radius': 'mean_radius'})

with pm.Model() as model_0:
    alpha = pm.Normal('alpha', mu=0, sd=10)
    beta = pm.Normal('beta', mu=0, sd=10)
    mu = alpha + pm.math.dot(x, beta)
    theta = pm.Deterministic('theta', 1 / (1 + pm.math.exp(-mu)))
    #bd:θ=0.5のときのmean_radius
    bd = pm.Deterministic('bd', -alpha/beta)
    yl = pm.Bernoulli('yl', p=theta, observed=y)

    chain = pm.sample(1000)

varnames = ['alpha', 'beta', 'bd']
pm.traceplot(chain, varnames)
pm.summary(trace)
plt.figure()

df = sns.load_dataset("tips")
x = df['total_bill']
y = df['tip']

with pm.Model() as model:
    alpha = pm.Normal('alpha', mu=0, sd=10)
    beta = pm.Normal('beta', mu=0, sd=10)
    theta = pm.Deterministic('theta', alpha + beta * x)
    y_pred = pm.Poisson('y_pred', mu=np.exp(theta), observed=y)
    trace = pm.sample(1000)

pm.traceplot(trace)
pm.summary(trace)
plt.figure()