1. はじめに

今回は、こちらの記事の続編として、反実仮想機械学習（CFML; CounterFactual Machine Learning）の技術要素の1つであるオフ方策評価（OPE; Off-Policy Evaluation）の実装例を紹介します。

2. オフ方策評価の実装例

ここでは Open Bandit Pipeline（OBP) を用いたPythonによるオフ方策評価の実装を紹介します。実装にあたっては、特にOBPのリポジトリの synthetic.ipynb および evaluate_off_policy_estimators.py を参考にしています。また、実行環境はPython 3.9.12（Windows 10）です。

2-1. 準備

2-1-1. ライブラリ（パッケージ）の準備

実装の準備として、まず最初に必要なライブラリを一通りimportします。installしていないパッケージに関しては、pipコマンドや各種パッケージで紹介されている方法（obpはこちらを参照）を用いて実行環境にinstallしてください。

import numpy as np
from sklearn.linear_model import LogisticRegression as Logit
from obp.dataset import SyntheticBanditDataset, logistic_reward_function
from obp.policy import IPWLearner, Random
from obp.ope import (
    OffPolicyEvaluation as OPE,
    RegressionModel as Reg_model,
    InverseProbabilityWeighting as IPS,
    DirectMethod as DM,
    DoublyRobust as DR,
    DoublyRobustWithShrinkage as DRos,
    SwitchDoublyRobust as SwitchDR,
    InverseProbabilityWeightingTuning as cIPSTuning,
    DoublyRobustTuning as DRpsTuning,
    DoublyRobustWithShrinkageTuning as DRosTuning,
    SwitchDoublyRobustTuning as SwitchDRTuning
)

2-1-2. データの準備

次に、今回使用する乱数および定数を指定します。今回は、行動の種類数が5、特徴量ベクトルの次元数が6という設定にしています。

RANDOM_STATE = 1017
N_actions = 5  # 行動(actions)の種類数
Dim = 6        # 特徴量ベクトルの次元数

そして、実装例に使用するデータとして、人工的なダミーデータを生成します。

# 報酬は2値(0 or 1)とし、その期待値はロジスティック関数で生成されるものとする
SBdata = SyntheticBanditDataset(
    n_actions=N_actions,
    dim_context=Dim,
    reward_type="binary",
    reward_function=logistic_reward_function,
    random_state=RANDOM_STATE
)
# 方策が決まっていない状態なので、方策で使用するモデルを作成するための学習用データを用意する
data_train = SBdata.obtain_batch_bandit_feedback(n_rounds=20000)
# 方策を評価するためのデータを（学習用データとは別で）用意する
data_valid = SBdata.obtain_batch_bandit_feedback(n_rounds=10000)

生成されたデータ（data_train, data_valid）の内容は以下の通りです。

'n_rounds'       : ログデータの行数
'n_actions'      : 行動(actions)の種類数
'context'        : 特徴量(context)の行列
'action_context' : 各行動を指す（one-hotの）ベクトルを格納したもの
'action'         : 方策によって取られた行動
'position'       : None（データ上に方策が複数存在する場合に、方策を区別するラベルを格納）
'reward'         : 行動によって得られた報酬
'expected_reward': 特徴量と各行動から計算される報酬の期待値
'pi_0'           : 各行動をとる確率
'pscore'         : 実際に行った行動をとる確率

'reward'および'expected_reward'は、特徴量、行動、および特徴量と行動の交互作用項にそれぞれランダムに生成した係数をかけて足し合わせたものをロジット変換して作成しています。つまり、交互作用項なしで単純にロジスティック回帰を用いるだけではうまく方策の価値を推定できないデータとしています。

以下のようなデータが生成されます。

print(data_train)
> Output exceeds the size limit. Open the full output data in a text editor
{'n_rounds': 20000, 'n_actions': 5, 'context': array([[ 2.09714853, -0.35579847,  0.06712557,  1.50564762, -0.85839701,
         1.55372896],
       [ 0.48818327,  1.19960695, -0.94140772,  0.32784131,  0.83420094,
         0.56572927],
       [ 0.41007327,  0.52990896, -0.13082159,  1.15790144, -2.04059223,
         1.01678568],
       ...,
       [-1.71399503,  0.17549435, -1.17718535, -0.0945791 ,  2.0135148 ,
        -1.72181698],
       [-0.77021031, -1.71381288,  0.20587305,  1.22108858,  2.05495588,
        -1.02054866],
       [ 1.16307083,  1.07968683,  0.81672033, -0.50489557,  1.83220528,
         0.80984235]]), 'action_context': array([[1, 0, 0, 0, 0],
       [0, 1, 0, 0, 0],
       [0, 0, 1, 0, 0],
       [0, 0, 0, 1, 0],
       [0, 0, 0, 0, 1]]), 'action': array([0, 4, 1, ..., 2, 0, 2], dtype=int64), 'position': None, 'reward': array([1, 1, 0, ..., 1, 1, 0]), 'expected_reward': array([[0.95646199, 0.79277703, 0.66569711, 0.32157778, 0.7671249 ],
       [0.72185158, 0.49735522, 0.75019125, 0.2280894 , 0.74151351],
       [0.73127674, 0.60793383, 0.30475919, 0.27827123, 0.38788611],
       ...,
       [0.33455863, 0.51759728, 0.76876385, 0.20858293, 0.67655739],
       [0.84164778, 0.8783893 , 0.85902601, 0.2545173 , 0.76676251],
       [0.7502211 , 0.42716675, 0.86110977, 0.268607  , 0.72296959]]), 'pi_0': array([[[0.25289532],
        [0.21471052],
        [0.18908767],
        [0.13403372],
        [0.20927278]],

       [[0.22432359],
        [0.17921626],
        [0.23077178],
        [0.1369105 ],
        [0.22877787]],

       [[0.25761148],
        [0.22771837],
        [0.16816321],
        [0.16376739],
        [0.18273955]],

       ...,

       [[0.16572778],
        [0.199016  ],
        [0.25583989],
        [0.14611162],
        [0.2333047 ]],

       [[0.22032619],
        [0.22857186],
        [0.22418853],
        [0.12248378],
        [0.20442964]],

       [[0.22562633],
        [0.16333868],
        [0.25208564],
        [0.13938865],
        [0.21956069]]]), 'pscore': array([0.25289532, 0.22877787, 0.22771837, ..., 0.25583989, 0.22032619,
       0.25208564])}

2-2. オフ方策評価による方策の価値の推定

2-2-1. 新しい方策の学習

次に、オフ方策評価をしたい新たな方策として、ここではロジスティック回帰にもとづくIPWLearner¹によって行動を選択する意思決定モデルを用います。具体的には、学習用のログデータを用いて行動選択に使用するモデルを構築します。このように、ログデータを用いて、方策における行動選択に使用する新たなモデルを構築することを、オフ方策学習（OPL; Off-Policy Learning） といいます。

# ベースラインとしてランダムに選択する方策を作成
rand_pol = Random(n_actions=N_actions, random_state=RANDOM_STATE)
action_dist_random = rand_pol.compute_batch_action_dist(
    n_rounds=data_valid["n_rounds"]
)
# ロジスティック回帰にもとづくIPWLearnerによるオフ方策学習(OPL)
ipw_logit = IPWLearner(n_actions=N_actions, base_classifier=Logit())
ipw_logit.fit(
    context=data_train["context"],
    action=data_train["action"],
    reward=data_train["reward"],
    pscore=data_train["pscore"]
)
action_dist_ipw = ipw_logit.predict(context=data_valid["context"])

2-2-2. 方策の真の価値

今回は人工的にデータを生成しており、報酬の真の期待値も得られるようになっているため、現行の方策および新たな方策（今回の場合はロジスティック回帰にもとづくIPWLearnerで行動を選択する方策）の真の価値が算出できます²。今回は、オフ方策評価の各手法が新たな方策の価値をうまく推定できるかを確認する実装例を紹介しますが、実務でオフ方策評価を行う際は、実データが必ずしも得られるわけではない点にご注意ください。

ground_truth_present_value = SBdata.calc_ground_truth_policy_value(
    expected_reward=data_valid["expected_reward"],
    action_dist=data_valid["pi_0"]
)
print("現行の方策の価値：", ground_truth_present_value)
ground_truth_random_value = SBdata.calc_ground_truth_policy_value(
    expected_reward=data_valid["expected_reward"],
    action_dist=action_dist_random
)
print("ランダム方策の価値：", ground_truth_present_value)
ground_truth_IPWLearner_value = SBdata.calc_ground_truth_policy_value(
    expected_reward=data_valid["expected_reward"],
    action_dist=action_dist_ipw
)
print("新たな方策の価値：", ground_truth_IPWLearner_value)

>現行の方策の価値： 0.5228788589361542
ランダム方策の価値： 0.5228788589361542
新たな方策の価値： 0.6831882078149638

生成したデータにおける現行の方策は、行動をランダムに選択しており、ランダムな方策で使用しているシード（および方法）が同じため、現行の方策の価値はランダムな方策の価値と一致します。また、新たな方策は現行の方策よりも価値が上がっていることから、新たな方策を用いることで、特徴量にもとづきより適切な行動が選択できるといえます。

2-2-3. 方策の価値の推定

新たな方策に使用する意思決定モデルができたので、ここからログデータをもとに新たな方策の価値を評価します。理論編で紹介した各種手法は、OBPライブラリで利用することができます。
ハイパーパラメータ $\lambda$ を持つ推定量に関しては、$\lambda=2.0$ で決め打ちした場合と、Tucker&Lee, 2021のSLOPE++アルゴリズムによりチューニングを行った場合で比較します。SLOPE++アルゴリズムは、推定した価値のBiasと偏差（deviation）を抑えられるようなハイパーパラメータを探索する方法です。

SLOPE++アルゴリズムの詳細を確認したい方はこちらをご参考ください。

ここでは、OBPでの実装に合わせて、SLOPE++アルゴリズムの詳細を説明します。以降では、探索したいハイパーパラメータの $K$ 個の候補を $\{\lambda_k\}_{k=1}^K$ とします。

1．パラメータ $\lambda_k$ をもとに価値を推定した場合のconfidence function（$CNF_k$）を計算します。confidence functionとは、パラメータ $\lambda_k$ をもとに価値を推定した場合の偏差以上の値を（少なくとも $1-\delta (\delta >0)$ の確率で）返す関数です。OBPでは、$i$ 番目のデータに対する報酬の推定値を計算したうえで標本分散を計算し、$1-\delta=0.05$ としたうえで、Thomas et al., 2015を参考に（標本分散の平方根）×（t分布の上位5%点）を $CNF_k$ としています³。

2．$CNF_k$ を値が大きい順に並べ替えます。並べ替えた後のものを $CNF_{(k)}$、対応するハイパーパラメータを $\lambda_{(k)}$、$\lambda_{(k)}$ をもとに推定した価値を $\hat V_{(k)}$ とします。

3．任意の $l<k$ に対し、

|\hat V_{(k)}-\hat V_{(l)}|\leq CNF_{(k)}+(\sqrt{6}-1)CNF_{(l)}

を満たす $k$ の最大値を探索し、求めた $k$ を $\hat k$ とします。

4．最適なパラメータとして $\lambda_{(\hat k)}$ を返します。

SLOPEアルゴリズム（Su et al., 2020b）と同様に、SLOPE++アルゴリズムは、MSEではなくMAEの分解を考えています。つまり、

|\hat V_{k}-V|\leq |E[\hat V_{k}]-V|+|\hat V_{k}-E[\hat V_{k}]|

としたうえで、Bias $|E[\hat V_{k}]-V|$ と偏差 $|\hat V_{k}-E[\hat V_{k}]|$ に分解したうえで、MAEを最小にするような $\lambda_k$ を探索しています。

# 報酬(reward)の推定モデル(ロジスティック回帰)を設定
logit_model = Reg_model(
    n_actions=N_actions,
    action_context=SBdata.action_context,
    base_model=Logit()
)
# Cross-fittingに基づく報酬の推定
estimated_rewards_logit_model = logit_model.fit_predict(
    context=data_valid["context"],
    action=data_valid["action"],
    reward=data_valid["reward"],
    n_folds=5,  # Cross-fittingにおけるデータの分割数
    random_state=RANDOM_STATE
)
# チューニングしたいハイパーパラメータを設定
lambdas = [0.1, 0.5, 1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0,
           200.0, 500.0, 1000.0, 2000.0, 5000.0, 10000.0, np.inf]
# 価値(Value)の推定量を指定
estimators = [
    DM(),
    IPS(estimator_name="ips"),
    DR(),
    IPS(lambda_=2.0, estimator_name="clipped-ips"),
    DR(lambda_=2.0, estimator_name="DRps"),
    DRos(lambda_=2.0, estimator_name="DRos"),
    SwitchDR(lambda_=2.0),
    IPSTuning(lambdas=lambdas, estimator_name="clipped-ips(tuned)"),
    DRTuning(lambdas=lambdas, estimator_name="DRps(tuned)"),
    DRosTuning(lambdas=lambdas, estimator_name="DRos(tuned)"),
    SwitchDRTuning(lambdas=lambdas, estimator_name="switch-dr(tuned)")
]
# 評価用データに対してオフ方策評価を実施
ope = OPE(
    bandit_feedback=data_valid,
    ope_estimators=estimators
)

ブートストラップ法により価値の推定を1000回繰り返した場合の平均と信頼区間を計算し、これを推定結果とします。

# 結果の要約（価値の推定値の平均と信頼区間を表示）
policy_value, interval = ope.summarize_off_policy_estimates(
    action_dist=action_dist_ipw,
    estimated_rewards_by_reg_model=estimated_rewards_logit_model,
    n_bootstrap_samples = 1000
)
print(interval, '\n')

>                        mean  95.0% CI (lower)  95.0% CI (upper)
dm                  0.562628          0.559643          0.565619
ips                 0.679760          0.649460          0.710363
dr                  0.684839          0.666926          0.703331
clipped-ips         0.326688          0.312600          0.341010
DRps                0.623109          0.614027          0.632170
DRos                0.575970          0.572426          0.579649
switch-dr           0.562737          0.559617          0.565576
clipped-ips(tuned)  0.679402          0.648790          0.710057
DRps(tuned)         0.685315          0.667023          0.704211
DRos(tuned)         0.667792          0.651954          0.684107
switch-dr(tuned)    0.685747          0.667748          0.702109

今回のデータでは、真の報酬は（交互作用項のない）単純なロジスティック回帰では表せないため、報酬の推定モデルのみで計算するDMでは、新しい方策の価値をうまく推定できていないことがわかります。また、$\lambda$ を決め打ちした推定量だと、新たな方策よりもランダムの方策の方が価値が高いなどの誤った結果が導かれる場合が見受けられるため、$\lambda$ のチューニングが必要だといえます。加えて、IPS推定量やcIPS推定量よりも、DR推定量やその発展型の方が信頼区間は狭くなっており、推定結果が安定しやすいことがわかります。

2-3. 推定精度の評価

最後に、推定精度を評価するため、Squared Error（SE）とrelative estimation error（relative_ee）を確認します。

SE= (\hat V_{\pi} - V_{\pi})^2

relative\_ee = |(\hat V_{\pi} - V_{\pi})/V_{\pi}|

どちらも0に近い値をとるほど、利用した評価手法による推定の精度が高いと解釈できます。

# Squared Error(SE)による精度評価
se_for_ipw = ope.summarize_estimators_comparison(
    ground_truth_policy_value=ground_truth_IPWLearner_value,
    action_dist=action_dist_ipw,
    estimated_rewards_by_reg_model=estimated_rewards_logit_model,
    metric="se"
)

print(se_for_ipw)
>                          se
dm                  0.014532
ips                 0.000012
dr                  0.000004
clipped-ips         0.127155
DRps                0.003629
DRos                0.011509
switch-dr           0.014532
clipped-ips(tuned)  0.000012
DRps(tuned)         0.000004
DRos(tuned)         0.000241
switch-dr(tuned)    0.000007

# relative estimation error(relative_ee)による精度評価
relative_ee_for_ipw = ope.summarize_estimators_comparison(
    ground_truth_policy_value=ground_truth_IPWLearner_value,
    action_dist=action_dist_ipw,
    estimated_rewards_by_reg_model=estimated_rewards_logit_model,
    metric="relative-ee"
)

print(relative_ee_for_ipw)
>                    relative-ee
dm                     0.176449
ips                    0.004994
dr                     0.002908
clipped-ips            0.521947
DRps                   0.088177
DRos                   0.157032
switch-dr              0.176449
clipped-ips(tuned)     0.004994
DRps(tuned)            0.002929
DRos(tuned)            0.022739
switch-dr(tuned)       0.003838

SEの結果をみると、DR推定量とチューニングしたDRps推定量にはほぼ差がなさそうに見えますが、relative_eeの結果を見ると、チューニングしたDRps推定量よりもDR推定量の方がわずかに高い精度を示すことがわかります。また、DRos推定量よりもDRps推定量やDR推定量の方が精度が高いことから、重み全体を調整するDRos推定量の方が、極端な重みにのみ対処するDRps推定量よりもバイアスがかかりやすい可能性があることが示唆されます。

3. まとめ

この記事では、反実仮想機械学習におけるオフ方策評価の実装について紹介しました。理論に加えて実際に実装を試すことで、理論上の数式から想定される現象を確認でき、より推定量に対する理解を深めるきっかけになるかと思います。例えば、2-3.の最後に考察した「DRps推定量よりもDRos推定量の方が重み全体を調整する分バイアスがかかりやすい可能性がある」といった内容は、数式だけでも想像できますが、実際に数値例を見るとわかりやすいです。また、今回は簡易のためデータは数万行にし、報酬の推定においてはロジスティック回帰を用いましたが、より大量のデータを用いたり、ランダムフォレストなど機械学習ベースのモデルに変えてみると、また違った結果になるかと思います。

理論編と今回の実装編が、オフ方策評価における理論と実装の両面から入門される方の参考になれば幸いです。

4. 参考文献

（英語（アルファベット順）→日本語（あいうえお順）で記載）

Cross-fittingに関する論文：Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W. and Robins, J. (2018), Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21: C1-C68.

Open Bandit Pipelineに関する文献：Saito, Y., Aihara, S., Matsutani, M., & Narita, Y. (2021). Open bandit dataset and pipeline: Towards realistic and reproducible off-policy evaluation. Proceedings of the Neural Information Processing Systems Track on Datasets and Benchmarks, 1.

clipped-IPSに関する文献：Strehl, A., Langford, J., Li, L., & Kakade, S. M. (2010). Learning from logged implicit exploration data. Advances in neural information processing systems, 23.

DRos、DRpsに関する文献(Su et al.(2020a))：Su, Y., Dimakopoulou, M., Krishnamurthy, A., & Dudik, M. (2020). Doubly robust off-policy evaluation with shrinkage. In International Conference on Machine Learning (pp. 9167-9176). PMLR.

SLOPEに関する文献(Su et al.(2020b))：Su, Y., Srinath, P., & Krishnamurthy, A. (2020). Adaptive estimator selection for off-policy evaluation. In International Conference on Machine Learning (pp. 9196-9205). PMLR.

Confidenceに関する文献：Thomas, P., Theocharous, G., & Ghavamzadeh, M. (2015). High confidence policy improvement. In International Conference on Machine Learning (pp. 2380-2388). PMLR.

SLOPE++に関する文献：Tucker, G., & Lee, J. (2021). Improved Estimator Selection for Off-Policy Evaluation. In Workshop on Reinforcement Learning Theory at the 38th International Conference on Machine Learning.

SwitchDRに関する論文：Wang, Y. X., Agarwal, A., & Dudık, M. (2017). Optimal and adaptive off-policy evaluation in contextual bandits. In International Conference on Machine Learning (pp. 3589-3597). PMLR.

反実仮想機械学習に関する書籍：齋藤優太, 安井翔太 (2021). 施策デザインのための機械学習入門データ分析技術のビジネス活用における正しい考え方, 技術評論社

今回は報酬が二値であることから、ロジスティック回帰にもとづくIPWLearnerは、通常のロジスティック回帰において（行動によって得られた報酬）/（実際に行った行動をとる確率）を重みとして使用することと同じです。（参考：齋藤, 安井 (2021)　p.77） ↩
厳密には、ログデータを所与とした場合の報酬の期待値を方策の価値と見做しています。そのため、報酬が決まる過程がオフ方策評価後に変容するとき（例えば、感染症の流行などによって消費者の行動傾向が大きく変容する出来事が起きたときなど）、実際に方策を適用した場合の価値と相違が出る可能性が高いです。 ↩
実装では $CNF_k$ は「報酬の価値の推定量-(t分布のlower bound)」とされていますが、t分布のlower boundが「報酬の価値の推定量-(標本分散の平方根)×(t分布の上位5%点)」となっているため、同じものを指しています。ただし、Thomas et al.(2015)ではt分布のlower boundの計算に関して、$i$番目のデータに対する報酬の推定値にバイアスがない(unbiasedである)ことを想定しているため、厳密にはバイアスを許容しているcIPSやDRosなどには不適である可能性があります。 ↩