Data
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
from sklearn.datasets import load_breast_cancer
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression
import pystan
cancer = load_breast_cancer()
indices = np.random.choice(len(cancer.data), 100, replace=False)
data = cancer.data[indices, :5]
target = cancer.target[indices]
print (cancer.feature_names[:5])
malignant_count = len(np.where(target==0)[0])
benign_count = len(np.where(target==1)[0])
print('# of 0 (malignant): ', malignant_count)
print('# of 1 (benign): ', benign_count)
x_train, x_test, y_train, y_test = train_test_split(data, target,
test_size=0.2, random_state=0)
print ('x_train shape: ', x_train.shape)
print ('y_train shape: ', y_train.shape)
print ('x_test shape: ', x_test.shape)
print ('y_test shape: ', y_test.shape)
Logistic Regression
clf = LogisticRegression()
clf.fit(x_train, y_train)
print ('Accuracy train: ', clf.score(x_train, y_train))
print ('Accuracy validation: ', clf.score(x_test, y_test))
Bayesian Logistic Regression
cancer_data = {'N': x_train.shape[0], 'M': x_train.shape[1], 'X': x_train, 'y': y_train}
lr_code = """
data {
int N;
int M;
real X[N, M];
int<lower=0, upper=1> y[N];
}
parameters {
real beta0;
real beta[M];
}
model {
for (i in 1:N)
y[i] ~ bernoulli(inv_logit (beta0 + dot_product(X[i] , beta)));
}
"""
%%time
stm = pystan.StanModel(model_code=lr_code)
%%time
n_itr = 2000
n_warmup = 500
chains = 3
fit = stm.sampling(data=cancer_data, iter=n_itr, chains=chains, n_jobs=-1,
warmup=n_warmup, algorithm="NUTS", verbose=False)
fit
la = fit.extract(permuted=True)
names = fit.model_pars
n_param = np.sum([1 if len(x) == 0 else x[0] for x in fit.par_dims])
mean_list = np.array(fit.summary()['summary'])[:,0]
f, axes = plt.subplots(n_param, 2, figsize=(10, 4*n_param))
cnt = 0
for name in names:
dat = la[name]
if dat.ndim == 2:
for j in range(dat.shape[1]):
d = dat[:,j]
sns.distplot(d, hist=False, rug=True, ax=axes[cnt, 0])
sns.tsplot(d, alpha=0.8, lw=.5, ax=axes[cnt, 1])
cnt += 1
else:
# Intercept
sns.distplot(dat, hist=False, rug=True, ax=axes[cnt, 0])
sns.tsplot(dat, alpha=0.8, lw=.5, ax=axes[cnt, 1])
cnt += 1
name_list = []
for name in names:
if la[name].ndim == 2:
for i in range(dat.shape[1]):
name_list.append("{}{}".format(name,i+1))
else:
name_list.append(name)
for i in range(2):
for j, t in enumerate(name_list):
axes[j, i].set_title(t)
plt.show()
def logistic(x, beta):
tmp = [1]
tmp.extend(x)
x = tmp
return (1+np.exp(-np.dot(x, beta)))**(-1)
def check_accuracy(data, target, param, threshold = 0.5):
ans_list = []
for i in range(len(data)):
res = logistic(data[i], param)
ans = 1 if res > threshold else 0
ans_list.append(ans == target[i])
return np.mean(ans_list)
param = mean_list[0:6]
print ('Accuracy train: ', check_accuracy(x_train, y_train, param))
print ('Accuracy test: ', check_accuracy(x_test, y_test, param))