傾向スコアを用いた逆確率重み付け法(IPTW)

# 乱数のシードを設定
import random
import numpy as np

np.random.seed(1234)
random.seed(1234)

# 使用するパッケージ（ライブラリと関数）を定義
# 標準正規分布の生成用
from numpy.random import *

# グラフの描画用
import matplotlib.pyplot as plt

# SciPy 平均0、分散1に正規化（標準化）関数
import scipy.stats

# シグモイド関数をimport
from scipy.special import expit

# その他
import pandas as pd

# データ数
num_data = 200

# 年齢
x_1 = randint(15, 76, num_data)  # 15から75歳の一様乱数

# 性別（0を女性、1を男性とします）
x_2 = randint(0, 2, num_data)  # 0か1の乱数


# x_1, x_2に依存する変数Zを作成する
# ノイズの生成
e_z = randn(num_data)

# シグモイド関数に入れる部分
z_base = x_1 + (1-x_2)*15 - 40 + 5*e_z

# シグモイド関数を計算
z_prob = expit(0.1*z_base)

# テレビCMを見たかどうかの変数（0は見ていない、1は見た）
Z = np.array([])

for i in range(num_data):
    Z_i = np.random.choice(2, size=1, p=[1-z_prob[i], z_prob[i]])[0]
    Z = np.append(Z, Z_i)


# 購入量Yを作成
# ノイズの生成
e_y = randn(num_data)

Y = -x_1 + 20*x_2 + 15*Z + 80 + 10*e_y


# データフレームにまとめる
df = pd.DataFrame({'年齢': x_1,
                   '性別': x_2,
                   'CMを見た': Z,
                   '購入量': Y,
                   })


# 平均値を確認
print('CMを見たデータ')
print(df[df["CMを見た"] == 1.0].mean())
print("--------")
print('CMを見てないデータ')
print(df[df["CMを見た"] == 0.0].mean())

# scikit-learnからロジスティク回帰をimport
from sklearn.linear_model import LogisticRegression

# 説明変数
X = df[["年齢", "性別"]]

# 被説明変数（目的変数）
Z = df["CMを見た"]

# 回帰の実施
reg = LogisticRegression().fit(X,Z)

# 回帰した結果の係数を出力
print("係数beta：", reg.coef_)
print("係数alpha：", reg.intercept_)

Z_pre = reg.predict_proba(X)
print('傾向スコア')
print(Z_pre[0:5])  # 5人ほどの結果を見てみる
print("----")
print('CMを見たかの正解')
print(Z[0:5])  # 5人ほどの正解

ATE_i = Y/Z_pre[:, 1]*Z - Y/Z_pre[:, 0]*(1-Z)
ATE = 1/len(Y)*ATE_i.sum()
print("推定したATE", ATE)