＜学習シリーズ＞線形モデル編：ロジスティック回帰

[IN]
%matplotlib qt
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import os
import glob

from sklearn import datasets


#データセットの読み込み
iris = datasets.load_iris() #irisデータセットの読み込み
data, target = iris.data, iris.target  #データとラベルを取得
data_all = np.concatenate([data, target[:, np.newaxis]], axis=1) #データとラベルを結合 axis=1:列方向
df = pd.DataFrame(data_all, columns=iris.feature_names+['label']) #データフレームに変換
id2label = {i:name for i, name in enumerate(iris.target_names)} #ラベルの辞書

print(df.shape)
print(id2label)

#データ可視化：2次元
fig, axs = plt.subplots(1, 2, figsize=(10, 6), facecolor='w')

df_setosa, df_versicolor, df_virginica = df[df['label']==0], df[df['label']==1], df[df['label']==2]
x1, x2, x3 = df_setosa['sepal length (cm)'], df_versicolor['sepal length (cm)'], df_virginica['sepal length (cm)']
y1, y2, y3 = df_setosa['label'], df_versicolor['label'], df_virginica['label']

#setosaとversicolorのみを抽出
axs[0].scatter(x1, y1, label='setosa', color='r', marker='o')
axs[0].scatter(x2, y2, label='versicolor', color='b', marker='x')
axs[0].set(xlabel='sepal length (cm)', ylabel='label', title='setosa and versicolor')
axs[0].grid(), axs[0].legend()

#3種類のデータをプロット
axs[1].scatter(x1, y1, label='setosa', color='r', marker='o')
axs[1].scatter(x2, y2, label='versicolor', color='b', marker='x')
axs[1].scatter(x3, y3, label='virginica', color='g', marker='^')
axs[1].set(xlabel='sepal length (cm)', ylabel='label', title='Plot of Iris dataset')
axs[1].grid(), axs[1].legend()

#全体のタイトル
fig.suptitle('Iris datasetの2次元プロット', fontsize=20)

#画像の保存
if not os.path.exists('output'): os.mkdir('images') #outputディレクトリがなければ作成
fig.savefig('output/iris_dataset_2d.png')

[OUT]
(150, 5)
{0: 'setosa', 1: 'versicolor', 2: 'virginica'}

[IN]
#データ可視化：3次元
z1, z2, z3 = df_setosa['label'], df_versicolor['label'], df_virginica['label']
fig = plt.figure(figsize=(10, 10), facecolor='w')
ax = Axes3D(fig)
ax.scatter(x1, y1, z1, label='setosa', color='r', marker='o')
ax.scatter(x2, y2, z2, label='versicolor', color='b', marker='x')
ax.scatter(x3, y3, z3, label='virginica', color='g', marker='^')
ax.set(xlabel='sepal length (cm)', ylabel='sepal width (cm)', zlabel='label', title='Iris datasetの3次元プロット')
plt.show()

[OUT]

[Terminal]
pip install ipympl

[IN]
def plotgraph(X, Y, xlabel, ylabel, title, xticks, yticks, verbose=False):
    fig, ax = plt.subplots(figsize=(8, 6), facecolor='white')
    ax.plot(X, Y, label='data')
    ax.set(xlabel=xlabel, ylabel=ylabel, title=title, xticks=xticks, yticks=yticks)
    if verbose:
        x_min, x_max, y_min, y_max = X.min(), X.max(), Y.min(), Y.max()
        ax.plot((x_min, x_max), (0, 0), color='black', lw=1) #x軸
        ax.plot((0, 0), (y_min, y_max), color='black', lw=1) #y軸
        
    ax.grid(), ax.legend()
    plt.show()

def logit(p):
    return np.log(p / (1 - p))

x = np.linspace(0.01, 0.99, 100) 
y = logit(x) #logit関数
plotgraph(x, y, '確率p', 'logit(p)', 'logit関数', np.arange(0, 1.1, 0.1), np.arange(-5, 6, 1), verbose=True)

[OUT]

[IN]
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LinearRegression, LogisticRegression

X = np.array([1, 2, 3, 4, 5])
Y = np.array([0, 0, 1, 1, 1])
# 最小二乗法による回帰分析
model1 = LinearRegression().fit(X.reshape(-1,1), Y)
y_pred = model1.predict(X.reshape(-1,1))
diffs = Y - y_pred # 残差

# 可視化
fig, ax = plt.subplots(1,1, figsize=(7,7), facecolor='w')
sns.set(style="darkgrid")
plt.scatter(X, Y)
# ロジスティック回帰による分類の可視化
ax.plot(X, Y, 'o', label="data", c='red')
ax.plot(X, y_pred, '-', label="logistic regression", c='blue', lw=1.0)
ax.set(xticks=X.flatten(), xlabel="x", ylabel="y", title="Logistic Regression")
#残差を両矢印で表示
for idx, diff in enumerate(diffs):
    ax.annotate("", xy=(X[idx], y_pred[idx]), xytext=(X[idx], Y[idx]), arrowprops=dict(arrowstyle="<->", color='blue'))
    ax.text(X[idx]+0.2, (Y[idx]+y_pred[idx])/2, f'{diff:.2f}', ha='center', va='center', color='blue')

ax.legend()
plt.show()

[OUT]

[IN]
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

#p(y=1|z)の尤度関数
def prob1(z):
    return 1 / (1 + np.exp(-z)) #z:線形予測値(w^TX)
#p(y=0|z)の尤度関数
def prob0(z):
    return 1 - prob1(z) #1-p(y=1|z)=exp(-z)/(1+exp(-z))

def likelihood(z, y):
    return prob1(z)**y * prob0(z)**(1-y)  #尤度関数 
#線形予測子z=w^TXを想定
z = np.linspace(-10, 10, 100)
#尤度関数の計算
likelihood_y0 = likelihood(z, 0) #y=0の尤度関数
likelihood_y1 = likelihood(z, 1) #y=1の尤度関数

#可視化
fig = plt.figure(figsize=(8, 6), facecolor='w')
plt.plot(z, likelihood_y0, label=r'$L(z|y=0)=\frac{e^{-z}}{1+e^{-z}}$', c='red')
plt.plot(z, likelihood_y1, label=r'$L(z|y=1)=\frac{1}{1+e^{-z}}$', c='blue')
plt.xlabel('線形予測子z', fontsize=12), plt.ylabel('尤度関数L', fontsize=12)
plt.title(r'ロジスティック回帰の尤度関数：$L(z|y)=(\frac{1}{1+e^{-z}})^y (\frac{e^{-z}}{1+e^{-z}})^{1-y}$', fontsize=16, y=-0.18)
plt.grid()
plt.legend(fontsize=14)
plt.show()

[OUT]

[IN]
%matplotlib qt
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
from mpl_toolkits.mplot3d import Axes3D
from sklearn.linear_model import LinearRegression
from scipy.optimize import minimize

# データの読み込み
data = np.array([[70, 80, 1],
                 [75, 65, 1],
                 [80, 90, 0],
                 [60, 40, 0],
                 [45, 55, 0],
                 [55, 70, 1],
                 [90, 95, 1],
                 [85, 80, 1],
                 [40, 30, 0],
                 [50, 50, 0]])

df = pd.DataFrame(data, columns=['英語', '数学', '合否'])
display(df.T)

X = data[:, :-1]  # 説明変数
y = data[:, -1]   # 目的変数

#重回帰で表現
linear = LinearRegression() #インスタンス化
linear.fit(X, y) #学習
y_pred_linear = linear.predict(X) #予測


#重回帰を可視化
fig = plt.figure(figsize=(10, 10), facecolor='w')
ax = Axes3D(fig)
ax.scatter(X[:, 0], X[:, 1], y, c='r', marker='o', label='ラベル', s=50)

x_min, x_max, y_min, y_max = X[:, 0].min(), X[:, 0].max(), X[:, 1].min(), X[:, 1].max()
X_mesh, y_mesh = np.meshgrid(np.linspace(x_min, x_max, 100), np.linspace(y_min, y_max, 100))
Z = linear.predict(np.c_[X_mesh.ravel(), y_mesh.ravel()])

ax.plot_surface(X_mesh, y_mesh, Z.reshape(X_mesh.shape), alpha=0.3, color='b')
ax.set(xlabel='英語', ylabel='数学', zlabel='合否')
plt.grid(), plt.legend()
plt.show()

[OUT]
スペースの都合上、転置して表示

[IN]
%matplotlib inline
p=np.linspace(0.01, 1, 100) #0.01~1の範囲で100個のデータを作成
y0, y1 = 0, 1 #ラベル値

plt.figure(figsize=(8, 6), facecolor='w')
plt.plot(p, y0*np.log(p) , label='ylog(p), y=0', c='red')
plt.plot(p, (1-y0)*np.log(1-p) , label='(1-y)log(1-p), y=0', c='red', ls='--')
plt.plot(p, y1*np.log(p) , label='ylog(p), y=1', c='blue')
plt.plot(p, (1-y1)*np.log(1-p) , label='(1-y)log(1-p), y=1', c='blue', ls='--')
plt.xlabel(r'確率$p$', fontsize=12)
plt.ylabel(r'$ylog(p)$+$(1-y)log(1-p)$', fontsize=12)
plt.grid(), plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0, fontsize=12)
plt.title('対数尤度関数', fontsize=12, y=-0.16)
plt.show()

[OUT]

[IN]
from sklearn.linear_model import LogisticRegression

log = LogisticRegression()
log.fit(X, y)
print(log.predict(X))
print(f'Weight\n{log.coef_}', f'bias\n{log.intercept_}', sep='\n')

[OUT]
[1 1 1 0 0 0 1 1 0 0]
Weight
[[0.04551016 0.05192708]]
bias
[-6.37705906]

[IN]
#ロジスティック回帰で表現
weights = np.array([-6.37705906, 0.04551016, 0.05192708])
def sigmoid(x):
    return 1 / (1 + np.exp(-x))
X_with_bias = np.c_[np.ones((X_mesh.ravel().shape[0], 1)), X_mesh.ravel(), y_mesh.ravel()]
_Z_proba = np.dot(X_with_bias, weights)
Z_proba = sigmoid(_Z_proba)

#重回帰を可視化
fig = plt.figure(figsize=(10, 10), facecolor='w')
ax = Axes3D(fig)
ax.scatter(X[:, 0], X[:, 1], y, c='r', marker='o', label='ラベル', s=50)

x_min, x_max, y_min, y_max = X[:, 0].min(), X[:, 0].max(), X[:, 1].min(), X[:, 1].max()
X_mesh, y_mesh = np.meshgrid(np.linspace(x_min, x_max, 100), np.linspace(y_min, y_max, 100))
ax.plot_surface(X_mesh, y_mesh, Z_proba.reshape(X_mesh.shape), alpha=0.3, color='b')
# ax.plot_surface(X_mesh, y_mesh, Z.reshape(X_mesh.shape), alpha=0.3, color='b') #重回帰の結果を重ねる

ax.set(xlabel='英語', ylabel='数学', zlabel='合否')
plt.grid(), plt.legend()
plt.show()

[OUT]

[IN]
data_1dim = np.array([[70, 1],
                    [75, 1],
                    [80, 0],
                    [60, 0],
                    [45, 0],
                    [55, 1],
                    [90, 1],
                    [85, 1],
                    [40, 0],
                    [50, 0]])

#データにバイアスを追加
_ones_bias = np.ones((data_1dim.shape[0], 1))
datas = np.concatenate([_ones_bias, data_1dim[:, :-1]], axis=1) #バイアス項用の1を追加
target = data_1dim[:, -1] #ラベル

# ロジスティック回帰モデルの定義：回答確認
log = LogisticRegression()
log.fit(data_1dim[:, :-1], data_1dim[:, -1])
print(log.predict(data_1dim[:, :-1]))
print(f'slope\n{log.coef_}', f'intercept\n{log.intercept_}', sep='\n')

[OUT]
[1 1 1 0 0 0 1 1 0 0]
slope
[[0.09370268]]
intercept
[-6.0906743]

[IN]

def linear_reg(X, w):
    return np.dot(X, w) #内積で線形予測値を計算

def sigmoid(x):
    return 1 / (1 + np.exp(-x)) #シグモイド関数

# 尤度関数の定義
def calc_likelihood(prob, label):
    likelihood = prob ** label * (1 - prob) ** (1 - label) #p^y(1-p)^(1-y)
    return likelihood

# 各重みにおける尤度関数の計算・配列化
slopes = np.linspace(-10, 10, 100)
intercepts = np.linspace(-10, 10, 100)
S, I = np.meshgrid(slopes, intercepts) #各重みにおける尤度関数の計算・配列化

outputs = [] #尤度関数を格納する配列

#slopes:100×intercepts:100=10000通りの重みの組み合わせを計算※zip()をそのまま使うと100通りになるため注意
for intercept in intercepts:
    for slope in slopes:
        w = np.array([intercept, slope]) #重みを[intercept, slope]の順に並べる
        linear_pred = linear_reg(datas, w) #線形予測値
        prob = sigmoid(linear_pred) #確率
        likelihoods = calc_likelihood(prob, target) #各データ尤度
        likelihood_func = np.prod(likelihoods)  #尤度関数         outputs.append(likelihood_func) #尤度関数を配列に追加
    
Z = np.array(outputs) #尤度関数の配列
Zs = Z.reshape(S.shape) #尤度関数の配列をグラフ用に整形(100×100)

#グラフの描画
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(10, 10), facecolor='w')
ax = fig.gca(projection='3d')

#尤度関数のプロット
surf = ax.plot_surface(S, I,  Zs , cmap='coolwarm', alpha=0.8)

#x, y軸のラベルとタイトルの設定
ax.set_xlabel('slope')
ax.set_ylabel('intercept')
ax.set_zlabel('Likelihood')
ax.set(title='Logistic Regression Likelihood Surface',
       xticks=np.arange(-10, 10, 1), yticks=np.arange(-10, 10, 1))
plt.title('Logistic Regression Likelihood Surface')

#カラーバーの表示
fig.colorbar(surf)
plt.show()

[OUT]

[IN]
%matplotlib inline
# 等高線図の描画
fig, ax = plt.subplots(figsize=(8, 8), facecolor='w')
# 等高線のプロット
CS = ax.contour(S, I, Zs, levels=10, cmap='coolwarm')

# 等高線のラベルの設定
ax.clabel(CS, inline=True, fontsize=10)

# 軸ラベルとタイトルの設定
ax.set_xlabel('slope', fontsize=12)
ax.set_ylabel('intercept', fontsize=12)
ax.set_title('Logistic Regression Likelihood Surface', fontsize=14)
ax.set(xticks=np.arange(-10, 10, 1), yticks=np.arange(-10, 10, 1))
# グリッド表示
ax.grid(True)
plt.show()

# 等高線図の描画
fig, ax = plt.subplots(figsize=(8, 8), facecolor='w')
# 等高線のプロット
CS = ax.contour(S, I, Zs, levels=10, cmap='coolwarm')

# 等高線のラベルの設定
ax.clabel(CS, inline=True, fontsize=10)

# 軸ラベルとタイトルの設定
ax.set_xlabel('slope', fontsize=12)
ax.set_ylabel('intercept', fontsize=12)
ax.set_title('Logistic Regression Likelihood Surface', fontsize=14)
ax.set(xticks=np.arange(-2, 2, 0.5), yticks=np.arange(-9, -4, 0.5), xlim=(-2, 2), ylim=(-9, -4))
# グリッド表示
ax.grid(True)
plt.show()

[OUT]

[IN]
class HorizontalDisplay:
    def __init__(self, *args):
        self.args = args
        
    def _repr_html_(self):
        template = '<div style="float: left; padding: 10px;">{}</div>'
        return "\n".join(template.format(a._repr_html_()) for a in self.args)

#データ数が少ないデータセットを作成
df_binary = df[df['label']!=2] #virginicaを除外

#乱数値を生成
np.random.seed(0) #乱数の種を設定
nums_index_binary = np.random.permutation(len(df_binary)) #0~データ数までの整数の乱数生成
nums_index_multi = np.random.permutation(len(df)) #0~データ数までの整数の乱数生成

#パターン1：1次元, 2値分類
df_1dim_binary = df_binary.iloc[nums_index_binary[:8], [0,4]] #sepal length (cm)とlabelのみを抽出
#パターン2：1次元, 3値分類
df_1dim_multi = df.iloc[nums_index_multi[:12], [0,4]] #sepal length (cm)とlabelのみを抽出
#パターン3：2次元, 2値分類
df_2dim_binary = df_binary.iloc[nums_index_binary[:8], [0,1,4]] #sepal length (cm)とsepal width (cm)とlabelのみを抽出
#パターン4：2次元, 3値分類
df_2dim_multi = df.iloc[nums_index_multi[:12], [0,1,4]] #sepal length (cm)とsepal width (cm)とlabelのみを抽出
display(HorizontalDisplay(df_1dim_binary, df_1dim_multi))
display(HorizontalDisplay(df_2dim_binary, df_2dim_multi))


#Excelファイルに出力
# dftoexcel = pd.concat([df_1dim_binary.reset_index(drop=True), df_1dim_multi.reset_index(drop=True), df_2dim_binary.reset_index(drop=True), df_2dim_multi.reset_index(drop=True)], axis=1)
# dftoexcel.to_excel('output/iris_dataset.xlsx')

[OUT]

[IN]
#Numpyでロジスティック回帰をスクラッチ実装

class LogisticScratch:
    def __init__(self, lr=0.01, iters=1000, random_state=1):
        self.lr = lr 
        self.iters = iters
        self.random_state = random_state

    def sigmoid(self, x): 
        return 1 / (1 + np.exp(-x)) #シグモイド関数
        
    def fit(self, X, y):
        rgen = np.random.RandomState(self.random_state) #乱数の初期化
        self.w_ = rgen.randn(X.shape[1])/1e5 #重みの初期化
        self.bias_ = rgen.randn(1)/1e5 #バイアスの初期化
        self.cost_ = [] #損失関数の誤差を格納するリスト
        
        for i in range(self.iters):
            input = np.dot(X, self.w_) + self.bias_ #入力値：重回帰モデルと同じ
            output = self.sigmoid(input) #活性化関数：シグモイド関数
            error = (y - output) #誤差：(正解データ-予測値)
            self.w_ = self.w_ + self.lr * np.dot(X.T, error) #重みの更新
            self.bias_ = self.bias_ + self.lr * np.sum(error) #バイアスの更新
            loss = -np.dot(y, np.log(output)) - np.dot((1-y), np.log(1-output)) #損失関数：交差エントロピー誤差
            self.cost_.append(loss) #損失関数の誤差を格納  
        return self
    
    def predict(self, X):
        return np.where(np.dot(X, self.w_) + self.bias_ > 0.5, 1, 0) #予測値：0.5以上なら1, そうでなければ0


#データセット1：1次元, 2値分類
print(f'{"#"*10} データセット1：1次元, 2値分類 {"#"*10}')
x = df_1dim_binary['sepal length (cm)'].values
x = x.reshape(-1, 1) #2次元配列に変換
y = df_1dim_binary['label'].values
print(f'x.shape：{x.shape}, y.shape：{y.shape}\n')

logistic = LogisticScratch(lr=0.05, iters=1000, random_state=1)
logistic.fit(x, y)
y_pred = logistic.predict(x)
df = pd.DataFrame({'y': y, 'y_pred': y_pred})
display(df.T)
print(f'正解率：{np.sum(y_pred==y)/len(y)}, weight：{logistic.w_}, bias：{logistic.bias_}\n')

#データセット2：1次元, 3値分類
print(f'{"#"*10} データセット2：1次元, 3値分類 {"#"*10}')
x = df_1dim_multi['sepal length (cm)'].values
x = x.reshape(-1, 1) #2次元配列に変換
y = df_1dim_multi['label'].values
print(f'x.shape：{x.shape}, y.shape：{y.shape}\n')

logistic = LogisticScratch(lr=0.05, iters=1000, random_state=1)
logistic.fit(x, y)
y_pred = logistic.predict(x)
df = pd.DataFrame({'y': y, 'y_pred': y_pred})
display(df.T)
print(f'正解率：{np.sum(y_pred==y)/len(y)}, weight：{logistic.w_}, bias：{logistic.bias_}\n')

#データセット3：2次元, 2値分類
print(f'{"#"*10} データセット3：2次元, 2値分類 {"#"*10}')
x = df_2dim_binary.iloc[:, :-1].values
y = df_2dim_binary.iloc[:, -1].values
print(f'x.shape：{x.shape}, y.shape：{y.shape}\n')

logistic = LogisticScratch(lr=0.05, iters=1000, random_state=1)
logistic.fit(x, y)
y_pred = logistic.predict(x)
df = pd.DataFrame({'y': y, 'y_pred': y_pred})
display(df.T)
print(f'正解率：{np.sum(y_pred==y)/len(y)}, weight：{logistic.w_}, bias：{logistic.bias_}\n')

#データセット4：2次元, 3値分類
print(f'{"#"*10} データセット4：2次元, 3値分類 {"#"*10}')
x = df_2dim_multi.iloc[:, :-1].values
y = df_2dim_multi.iloc[:, -1].values
print(f'x.shape：{x.shape}, y.shape：{y.shape}\n')

logistic = LogisticScratch(lr=0.05, iters=1000, random_state=1)
logistic.fit(x, y)
y_pred = logistic.predict(x)
df = pd.DataFrame({'y': y, 'y_pred': y_pred})
display(df.T)
print(f'正解率：{np.sum(y_pred==y)/len(y)}, weight：{logistic.w_}, bias：{logistic.bias_}\n')

[OUT]

[IN]
np.set_printoptions(precision=3, suppress=True) #小数点以下3桁, 指数表記を抑制
from sklearn.linear_model import LogisticRegression

#データセット1：1次元, 2値分類
print(f'{"#"*10} データセット1：1次元, 2値分類 {"#"*10}')
X, y = df_1dim_binary.iloc[:, :-1], df_1dim_binary.iloc[:, -1] #データとラベルを取得
logistic = LogisticRegression() #インスタンスを生成
logistic.fit(X, y) #モデルを学習
y_pred = logistic.predict(X) #予測
df = pd.DataFrame({'y': y, 'y_pred': y_pred}) #データフレームに変換
display(df.T)
print(f'正解率：{logistic.score(X, y):.2f}, weight：{logistic.coef_}, bias：{logistic.intercept_}\n')

#データセット2：1次元, 3値分類
print(f'{"#"*10} データセット2：1次元, 3値分類 {"#"*10}')
X, y = df_1dim_multi.iloc[:, :-1], df_1dim_multi.iloc[:, -1] #データとラベルを取得
logistic = LogisticRegression() #インスタンスを生成
logistic.fit(X, y) #モデルを学習
y_pred = logistic.predict(X) #予測
df = pd.DataFrame({'y': y, 'y_pred': y_pred}) #データフレームに変換
display(df.T)
print(f'正解率：{logistic.score(X, y):.2f}, weight：{logistic.coef_.reshape(1,-1)}, bias：{logistic.intercept_}\n')

#データセット3：2次元, 2値分類
print(f'{"#"*10} データセット3：2次元, 2値分類 {"#"*10}')
X, y = df_2dim_binary.iloc[:, :-1], df_2dim_binary.iloc[:, -1] #データとラベルを取得
logistic = LogisticRegression() #インスタンスを生成
logistic.fit(X, y) #モデルを学習
y_pred = logistic.predict(X) #予測
df = pd.DataFrame({'y': y, 'y_pred': y_pred}) #データフレームに変換
display(df.T)
print(f'正解率：{logistic.score(X, y):.2f}, weight：{logistic.coef_.reshape(1,-1)}, bias：{logistic.intercept_}\n')

#データセット4：2次元, 3値分類
print(f'{"#"*10} データセット4：2次元, 3値分類 {"#"*10}')
X, y = df_2dim_multi.iloc[:, :-1], df_2dim_multi.iloc[:, -1] #データとラベルを取得
logistic = LogisticRegression() #インスタンスを生成
logistic.fit(X, y) #モデルを学習
y_pred = logistic.predict(X) #予測
df = pd.DataFrame({'y': y, 'y_pred': y_pred}) #データフレームに変換
display(df.T)
print(f'正解率：{logistic.score(X, y):.2f}, weight：{logistic.coef_.reshape(1,-1)}, bias：{logistic.intercept_}\n')

[OUT]

[IN]
[OUT]

[IN]
import numpy as np
import matplotlib.pyplot as plt

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

def diff_sigmoid(x):
    return (1.0 - sigmoid(x)) * sigmoid(x)

x = np.arange(-5.0, 5.0, 0.1)
fig, ax = plt.subplots(1, 1, figsize=(8, 5))
ax.plot(x, sigmoid(x), label='sigmoid')
ax.plot(x, diff_sigmoid(x), label='\u0394sigmoid')
ax.set(xlabel='x', ylabel='y', title='Sigmoid function',
       xlim=(-5, 5), ylim=(0.0, 1.0), xticks=np.arange(-5, 6, 1), yticks=np.arange(0, 1.1, 0.1))
ax.grid(), ax.legend()
plt.show()

[OUT]

[IN]
#ロジスティック関数
def logisticfunc(x:float, a:float, b:float, c:float):
    return a / (1 + np.exp(-b * (x - c)))

x = np.linspace(-5, 5, 100)
y1 = logisticfunc(x, 1, 1, 0) #a=1, b=1, c=0: シグモイド関数と同じ
y2= logisticfunc(x, 1, 1, 3) #a=1, b=1, c=3: シグモイド関数より右にずらした
y3 = logisticfunc(x, 2, 1, 0) #a=2, b=1, c=0: シグモイド関数より上にずらした
y4 = logisticfunc(x, 1, 3, 0) #a=1, b=3, c=0: シグモイド関数より上にずらした

fig, ax = plt.subplots(figsize=(8, 6), facecolor='white')
ax.plot(x, y1, label='a=1, b=1, c=0: Sigmoid')
ax.plot(x, y2, label='a=1, b=1, c=3')
ax.plot(x, y3, label='a=3, b=1, c=0')
ax.plot(x, y4, label='a=1, b=3, c=0')
ax.plot((0, 0), (0, 10), color='black', lw=1), ax.plot((-5, 5), (0, 0), color='black', lw=1) #x軸, y軸
ax.set(xlabel='x', ylabel='y', title='ロジスティック関数', xticks=np.arange(-5, 6, 1), yticks=np.arange(0, 2, 0.2), ylim=(0, 2))
ax.grid(), ax.legend()
plt.show()

[OUT]

[IN]
def logit_func(x):
    return np.log(x / (1 - x))

x = np.arange(0.01, 1.0, 0.01)
fig, ax = plt.subplots(1, 1, figsize=(8, 5), facecolor='w')
ax.plot(x, logit_func(x), label='logit')
ax.plot((0, 1), (0, 0), color='black', lw=1), ax.plot((0, 0), (-5, 5), color='black', lw=1) #x軸とy軸
ax.set(xlabel='x', ylabel='y', title='Logit function', xticks=np.arange(0, 1.1, 0.1), yticks=np.arange(-5, 6, 1))
ax.grid(), ax.legend()
plt.show()

[OUT]

[IN]
#オッズ
p = np.arange(0, 0.95, 0.01) #確率
p_inv = 1 - p
odds = p / p_inv #オッズ

#グラフ化
fig, ax = plt.subplots(1, 1, figsize=(6, 6), facecolor='w')
ax.plot(p, odds, label='odds')
ax.set(xlabel='確率p', ylabel='odds', title='Odds function', xticks=np.arange(0, 1.1, 0.1), yticks=np.arange(0, 16.1, 1))
ax.grid(), ax.legend()
plt.show()

[OUT]

[IN]
class CoinToss:
    def __init__(self, nums_play: int, nums_toss: int, p: float):
        self.nums_play = nums_play #試行回数(サンプル数)
        self.nums_toss = nums_toss #コイントス回数(サンプルサイズ)
        self.p = p #表の確率
        self.values = np.array([0, 1]) #0: 裏, 1: 表
        
    def __call__(self):
        outputs = [] #全結果
        for _ in range(self.nums_play):
            # n回目のコイントス(試行)
            counts = np.random.choice(self.values, self.nums_toss, p=[1 - self.p, self.p]) #play
            counts_head = np.sum(counts) #表の回数
            outputs.append([counts, counts_head])
        return outputs
    
#コイントス
np.random.seed(12) #乱数の固定
nums_play, nums_toss, p = 100, 10, 0.7 #プレイ数：50, コイントス回数：10, 表の確率：0.7
cointoss = CoinToss(nums_play, nums_toss, p)
outputs = cointoss() #出力：試行結果, 表の回数
heads = [output[1] for output in outputs] #表の回数のみ

#表の回数を集計
counts_head = np.bincount(heads, minlength=11) #表が出た回数
props_head = counts_head / np.sum(counts_head) #表が出た確率

#リンク関数（ロジット関数）による確率変換
def logit(p):
    return np.log(p / (1 - p)) #ロジット関数

logit_probs = logit(props_head[1:]) #ロジット関数による確率変換

#グラフ化
X = np.arange(1, nums_toss+1, 1) #x軸: 表の回数
print(f'形状確認：X={X.shape}, props_head={props_head.shape}, logit_probs={logit_probs.shape}')

fig, axs = plt.subplots(2, 2, figsize=(12,12), facecolor='w')
axs = axs.ravel() #2次元配列を1次元配列に変換
axs[0].bar(X, props_head[1:], align='center', width=0.5, label='counts')
axs[0].set(xlabel=f'コイントス回数{nums_toss}：表の回数', ylabel='表の出現確率:p', title=f'コイントス試行(試行数:{nums_play}, p={p})', xticks=np.arange(1, 11, 1), yticks=np.arange(0, 1.1, 0.1))
[axs[0].text(x-0.5, y+0.02, f'p={y:.2f}', fontsize=8) for x, y in zip(X, props_head[1:])] #確率の表示
axs[0].grid(), axs[0].legend()

axs[1].plot(np.linspace(0, 1, 100), logit(np.linspace(0, 1, 100)), label='logit')
axs[1].set(xlabel='確率p', ylabel='log(p/(1-p))', title='ロジット関数', xticks=np.arange(0, 1.1, 0.1), yticks=np.arange(-5, 6, 1))
axs[1].grid(), axs[1].legend()

axs[2].plot(X, logit_probs, label='logit')
axs[2].set(xlabel=f'コイントス回数{nums_toss}：表の回数', ylabel='log(p/(1-p))', title='リンク関数(ロジット関数)による確率変換:回数 vs logit', xticks=np.arange(1, 11, 1), yticks=np.arange(-5, 6, 1))
axs[2].grid(), axs[2].legend()

axs[3].scatter(props_head[1:], logit_probs, label='logit')
axs[3].set(xlabel='表の出現確率:p', ylabel='log(p/(1-p))', title='リンク関数(ロジット関数)による確率変換:p vs logit', xticks=np.arange(0, 1.1, 0.1), yticks=np.arange(-5, 6, 1))
axs[3].grid(), axs[3].legend()

fig.savefig('output/link関数の理解.png')

[OUT]

[IN]
def softmax(x):
    c = np.max(x)
    exp_x = np.exp(x - c)
    sum_exp_x = np.sum(exp_x)
    y = exp_x / sum_exp_x
    return y

[OUT]