「自然科学の統計学」第8章演習問題5-市場調査においてプロビットモデル・ロジットモデルをPythonで実装してみた

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import py4macro
import wooldridge

from py4etrics.hetero_test import *
from scipy.stats import logistic, norm, chi2
from statsmodels.formula.api import ols, logit, probit

# 警告メッセージを非表示
import warnings
warnings.filterwarnings("ignore")

df = pd.read_csv('market_research.csv', index_col = 0)
display(df)

### (i)
# 回帰式の立式
formula = '車購入の有無 ~ 収入 + 性別 + 居住地域'

# プロビットモデルの導出 
res_probit = probit(formula, data =df).fit()
print(res_probit.summary())

print()
print()

# ロジットモデルの導出
res_logit = logit(formula, data = df).fit()
print(res_logit.summary())

### (iii)
# プロビットモデルの係数を使って、回帰式を作成する
def probit_formula(income, sex, region):
    y = 0.004885  * income -1.3460 * sex -0.4042 * region -1.4173
    return y

# ロジットモデルの係数を使って、回帰式を作成する
def logit_formula(income, sex, region):
    y = 0.009059  * income -2.3155 * sex -0.5922 * region -2.8062
    return y

# 具体的な値を代入
probit = probit_formula(500, 0, 0)
# 標準正規分布の累積分布関数に代入する
pi = norm.cdf(probit)

print(f'プロビットモデルでのPi: {pi}')

### 分散を求める式が間違いの可能性あり
sigma = pi * (1 -  pi)/2

print(f'プロビットモデルでの分散: {sigma}')

# 信頼区間の導出
lower = pi - 1.96 * np.sqrt(sigma)
upper = pi + 1.96 * np.sqrt(sigma)
if lower <0:
    lower = 0
if upper > 1:
    upper = 1
print(f'プロビットモデルによる信頼区間: [{lower}, {upper}]')

print()

# 具体的な値を代入
logit = logit_formula(500, 0, 0)
# ロジスティック分布の累積分布関数に代入する
pi = logistic.cdf(logit)

print(f'ロジットモデルでのPi: {pi}')

### 分散を求める式が間違いの可能性あり
sigma = pi * (1 -  pi)/2

print(f'ロジットモデルでの分散: {sigma}')

# 信頼区間の導出
lower = pi - 1.96 * np.sqrt(sigma)
upper = pi + 1.96 * np.sqrt(sigma)
if lower <0:
    lower = 0
if upper > 1:
    upper = 1
print(f'ロジットモデルによる信頼区間: [{lower}, {upper}]')

df = pd.read_csv('market_research.csv', index_col = 0)
display(df)

### (i)
# 回帰式の立式
formula = '車購入の有無 ~ 収入 + 性別 + 居住地域'

# プロビットモデルの導出 
res_probit = probit(formula, data =df).fit()
print(res_probit.summary())

print()
print()

# ロジットモデルの導出
res_logit = logit(formula, data = df).fit()
print(res_logit.summary())

### (iii)
# プロビットモデルの係数を使って、回帰式を作成する
def probit_formula(income, sex, region):
    y = 0.004885  * income -1.3460 * sex -0.4042 * region -1.4173
    return y

# ロジットモデルの係数を使って、回帰式を作成する
def logit_formula(income, sex, region):
    y = 0.009059  * income -2.3155 * sex -0.5922 * region -2.8062
    return y

# 具体的な値を代入
probit = probit_formula(500, 0, 0)
# 標準正規分布の累積分布関数に代入する
pi = norm.cdf(probit)

print(f'プロビットモデルでのPi: {pi}')

### 分散を求める式が間違いの可能性あり
sigma = pi * (1 -  pi)/2

print(f'プロビットモデルでの分散: {sigma}')

# 信頼区間の導出
lower = pi - 1.96 * np.sqrt(sigma)
upper = pi + 1.96 * np.sqrt(sigma)
if lower <0:
    lower = 0
if upper > 1:
    upper = 1
print(f'プロビットモデルによる信頼区間: [{lower}, {upper}]')

print()

# 具体的な値を代入
logit = logit_formula(500, 0, 0)
# ロジスティック分布の累積分布関数に代入する
pi = logistic.cdf(logit)

print(f'ロジットモデルでのPi: {pi}')

### 分散を求める式が間違いの可能性あり
sigma = pi * (1 -  pi)/2

print(f'ロジットモデルでの分散: {sigma}')

# 信頼区間の導出
lower = pi - 1.96 * np.sqrt(sigma)
upper = pi + 1.96 * np.sqrt(sigma)
if lower <0:
    lower = 0
if upper > 1:
    upper = 1
print(f'ロジットモデルによる信頼区間: [{lower}, {upper}]')

### (iii)
# プロビットモデルの係数を使って、回帰式を作成する
def probit_formula(income, sex, region):
    y = 0.004885  * income -1.3460 * sex -0.4042 * region -1.4173
    return y

# ロジットモデルの係数を使って、回帰式を作成する
def logit_formula(income, sex, region):
    y = 0.009059  * income -2.3155 * sex -0.5922 * region -2.8062
    return y

# 具体的な値を代入
probit = probit_formula(500, 0, 0)

...

# 具体的な値を代入
logit = logit_formula(500, 0, 0)

# 標準正規分布の累積分布関数に代入する
pi = norm.cdf(probit)

print(f'プロビットモデルでのPi: {pi}')

...

# ロジスティック分布の累積分布関数に代入する
pi = logistic.cdf(logit)

print(f'ロジットモデルでのPi: {pi}')