<jupyter_start><jupyter_text>bap_logistic_regression> API details.<jupyter_code>#hide
from nbdev.showdoc import *
%matplotlib inline
import numpy as np
import pandas as pd
import pymc3 as pm
import theano.tensor as tt
import matplotlib.pyplot as plt
import scipy.stats as stats
import arviz as az
from scipy.special import expit as logistic
import seaborn as sns
sns.set_context('talk')
sns.set_style('white')

RANDOM_SEED = 123<jupyter_output>WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.<jupyter_text>1. Logistic regression  The IRIS dataset<jupyter_code># iris = sns.load_dataset("iris")
iris = pd.read_csv("data/iris.csv")
iris.head()
iris.shape
sns.set(rc={'figure.figsize':(12,8)})
sns.stripplot(x="species", y="sepal_length", data=iris, jitter=True)
sns.pairplot(iris, hue='species', diag_kind='kde')<jupyter_output><empty_output><jupyter_text>The logistic model<jupyter_code># Data preprocessing
df = iris.query("species == ('setosa', 'versicolor')")
y_0 = pd.Categorical(df['species']).codes
x_label = 'sepal_length'
x_0 = df[x_label].values
x_center = x_0 - x_0.mean()
df.shape
# model
with pm.Model() as logit_0:
    α = pm.Normal('α', mu=0, sd=10)
    β = pm.Normal('β', mu=0, sd=10)
    
    μ = α + pm.math.dot(x_center, β)
    #θ = pm.Deterministic('θ', pm.math.sigmoid(μ))
    θ = pm.Deterministic('θ', 1 / (1 + pm.math.exp(-μ)))
    bd = pm.Deterministic('bd', -α/β)
    
    yl = pm.Bernoulli('yl', p=θ, observed=y_0)
    
    trace_0 = pm.sample(2000)
az.summary(trace_0, var_names=['α', 'β', 'bd'])
pm.traceplot(trace_0, ['α', 'β', 'bd'])
# Visual posterior dist and observed data
theta = trace_0['θ'].mean(axis=0)
idx = np.argsort(x_center)
plt.plot(x_center[idx], theta[idx], color='C2', lw=3)
plt.vlines(trace_0['bd'].mean(), 0, 1, color='k')
bd_hpd = az.hpd(trace_0['bd'])
plt.fill_betweenx([0, 1], bd_hpd[0], bd_hpd[1], color='k', alpha=0.5)

plt.scatter(x_center, np.random.normal(y_0, 0.02),
            marker='.', color=[f'C{x}' for x in y_0])
az.plot_hpd(x_center, trace_0['θ'], color='C2')

plt.xlabel(x_label)
plt.ylabel('θ', rotation=0)
# use original scale for xticks
locs, _ = plt.xticks()
plt.xticks(locs, np.round(locs + x_0.mean(), 1))<jupyter_output>/opt/conda/lib/python3.7/site-packages/arviz/stats/stats.py:338: UserWarning: hpd will be deprecated Please replace hdi
  warnings.warn(("hpd will be deprecated " "Please replace hdi"),)
/opt/conda/lib/python3.7/site-packages/arviz/data/base.py:146: UserWarning: More chains (8000) than draws (100). Passed array should have shape (chains, draws, *shape)
  UserWarning,<jupyter_text>Prediction<jupyter_code>def classify(n, threshold):
    n = np.array(n)
    mu = trace_0['α'].mean() + trace_0['β'].mean() * n
    prob = 1 / (1 + np.exp(-mu))
    return prob, prob > threshold
# Test
classify([5, 5.5, 6], 0.4)
y_prob, y_pred = classify(x_center, 0.5)
accuracy = (y_pred == (y_0 == 1)).mean()
accuracy<jupyter_output><empty_output><jupyter_text>2. Multiple logistic regression  Data preparation<jupyter_code>df = iris.query("species == ('setosa', 'versicolor')")
y_1 = pd.Categorical(df['species']).codes
x_n = ['sepal_length', 'sepal_width']
x_1 = df[x_n].values<jupyter_output><empty_output><jupyter_text>Model<jupyter_code>with pm.Model() as multiple_logit:
    α = pm.Normal('α', mu=0, sd=10)
    β = pm.Normal('β', mu=0, sd=2, shape=len(x_n))
    
    μ = α + pm.math.dot(x_1, β)
    θ = pm.Deterministic('θ', 1 / (1 + pm.math.exp(-μ)))
    bd = pm.Deterministic('bd', -α/β[1] - β[0]/β[1] * x_1[:,0])
    
    yl = pm.Bernoulli('yl', p=θ, observed=y_1)
    
    trace_1 = pm.sample(2000)
pm.traceplot(trace_1, var_names=['α', 'β', 'bd'])
var_names = ['α', 'β']
az.plot_forest(trace_1, var_names=var_names)
az.summary(trace_1, var_names=var_names)
# Plot posterior dist with observed data
idx = np.argsort(x_1[:,0]) 
bd = trace_1['bd'].mean(0)[idx] 
plt.scatter(x_1[:,0], x_1[:,1], c=[f'C{x}' for x in y_0]) 
plt.plot(x_1[:,0][idx], bd, color='k'); 

az.plot_hpd(x_1[:,0], trace_1['bd'], color='k')
 
plt.xlabel(x_n[0]) 
plt.ylabel(x_n[1])
plt.savefig('B11197_04_05.png', dpi=300);<jupyter_output>/opt/conda/lib/python3.7/site-packages/arviz/data/base.py:146: UserWarning: More chains (8000) than draws (100). Passed array should have shape (chains, draws, *shape)
  UserWarning,<jupyter_text>Prediction<jupyter_code>def classify(x_1, threshold):
    x_1 = np.array(x_1)
    mu = trace_1['α'].mean() + trace_1['β'][:,0].mean() * x_1[:,0] + trace_1['β'][:,1].mean() * x_1[:,1]
    prob = 1 / (1 + np.exp(-mu))
    return prob, prob > threshold
def classify2(x_1, threshold):
    x_1 = np.array(x_1)
    mu = -0.14344304026621013 + 0.92119197 * x_1[:,0] + -1.55931326 * x_1[:,1]
    prob = 1 / (1 + np.exp(-mu))
    return prob, prob > threshold
y_prob, y_pred = classify(x_1, 0.5)
accuracy = (y_pred == (y_0 == 1)).mean()
accuracy
y_prob, y_pred = classify2(x_1, 0.5)
accuracy = (y_pred == (y_0 == 1)).mean()
accuracy<jupyter_output><empty_output><jupyter_text>Interpreting the coefficients of a logistic regression<jupyter_code>probability = np.linspace(0.01, 1, 100)
odds = probability / (1 - probability)

_, ax1 = plt.subplots()
ax2 = ax1.twinx()
ax1.plot(probability, odds, 'C0')
ax2.plot(probability, np.log(odds), 'C1')

ax1.set_xlabel('probability')
ax1.set_ylabel('odds', color='C0')
ax2.set_ylabel('log-odds', color='C1')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.187175685693106, 'fraction_numerical': 0.04447739065974796, 'mean_word_length': 5.057239057239057, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 34, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '446300', 'license': 'Apache-2.0', 'n_tokens_mistral': 2401, 'n_tokens_neox': 2273, 'n_words': 436, 'path': '22_bap_logistic_regression.ipynb', 'repo_name': 'danhphan/bayesian_inferences', 'size': '5396'}