freeeデータに関わる人たち Advent Calendar 2020の13日目、analyticsのxmmmxです。

背景

直前で他のアドベントカレンダーの方と完全にネタが被ってしまいどうしようか悩んでいたところ以下の記事を見つけました。

初めて聞く用語だったのでググってみると、

「原因から結果を導く普通の因果律とは逆に、結果から原因を推定する」ものを言う
引用:逆解析とその応用

結果から原因を推定？これは因果推論に関心のある私にとっては非常に興味深いトピックなので早速試してみることにしました。
が、上記の記事の中で使われているソルバーは有償なのと他の数理最適化で使われるソルバーは使ったことがないためOptunaで実装することにしました、　

目次

逆解析とは
やってみた
課題
参考

逆解析とは

y=x1+x2+x3という式があるとき、yの値(この記事では住宅の価格:price)から任意の数の説明変数を推定すること。
詳しくは以下を参照。
[機械学習モデルを逆解析する] (https://qiita.com/futakuchi0117/items/92d03c61f715d6b2b7c6)
機械学習モデルの逆解析のススメ

やってみた

使うデータはおなじみのkaggleのチュートリアルで使用されるボストン住宅価格のデータセットです。
逆解析を行う上での注意点はこちらの「モデルの逆解析をするときのチェックリスト」を参考にしました。ただ、正規化などの一部手順はすっとばしているので実装する際は注意して下さい。
実際の流れ
①データの読み込み
②目的変数yを予測するモデルの作成
③モデルと実測値の絶対誤差を最小化させるような説明変数の推定。

データの読み込み〜モデル作成

モデリングは単純な線形回帰にしますが、scikit-learnの他のモデルやlightgbmなどの他のライブラリを使うことも可能です。
このくらいならexcelで出来るのでより複雑なモデルが適当
今回はとりあえず回帰係数の一番大きかった1住戸あたりの平均部屋数(RM)を推定します。
*正規化はしてない。
今回に限らず、推定したい変数がモデルに影響を与えてなければ意味がないのでまずは回帰係数や変数重要度を確認する必要があると思います。

import sklearn.datasets
import sklearn.model_selection
import optuna
from optuna.logging import set_verbosity
import pandas as pd
import numpy  as np
from sklearn import linear_model

boston = sklearn.datasets.load_boston()
boston_df=(pd.DataFrame(boston.data, columns=boston.feature_names)
            .assign(price=boston.target)
            )

X=boston_df.drop(["price"],axis=1)
y=boston_df["price"]

model = linear_model.LinearRegression()
model.fit(X, y)

# 訓練データを用いた評価
print(model.score(X, y))
 
# 回帰係数
coeff_df=pd.DataFrame({"Feature":X.columns,"Coefficients":np.transpose(model.coef_)}).sort_values(by=["Coefficients"],ascending=False)
coeff_df

推定用の関数定義

高階関数にしているのはDataFrameを変数として推定の中で使うためです。
目的関数は予測結果y_predと実際の値yの絶対誤差です。この誤差を最小化するような説明変数rmを探索します。
今回は変数1つだけですが、モデルによっては多変量でも出来ます。

def objective_variable(df):

    def objective(trial):
    
        #ハイパーパラメータの定義(元の値のmin-max)
        # age=trial.suggest_int('age', 1, 65)
        rm=trial.suggest_float("rm", 3.56, 8.78,step=0.01)

        # X, yで分割
        X=df.drop(["price"],axis=1)
        y=df["price"]

        # 元のdfに代入
        # X["NOX"]=indus
        X["RM"]=rm

        pred_y=model.predict(X)

        return abs(y-pred_y)

パラメータ推定のテスト実行

基本的にDataFrameを一行ずつ渡して推定することになると思います。まずは1レコードでやってみます。
試しに推定してみるレコードは↓です。推定したいRMは6.575となっています。

# テストログを全部出す
set_verbosity(optuna.logging.DEBUG)
# お試しで1回実行してみる　
study = optuna.create_study(direction='minimize' )
study.optimize(objective_variable(boston_df.head(1)), n_trials=100)

#トライアル結果
trial = study.best_trial

#ベストな結果、パラメータの表示
print(f"Best hyperparameters: {trial.params}")

結果はBest hyperparameters: {'rm': 4.99}になりました。元のパラメータが6.575なので微妙なところ。
単純な線形回帰モデルで予測精度がそこまでよろしくないのでこのような結果になりました。実際に使う場合は変数の前処理や予測モデルと目的関数のデザインをきちんと行う必要があります。

パラメータ推定の本番実行

DataFrameをインデックスでスライシングして一行ずつ推定させています。もっとエレガントなやり方が知りたい・・
実際に使う際は推定したパラメータをどこかに保存するかDataFrameとして持っておくことになると思います。

# 試行中のログはオフにしておかないとレンダリングで落ちる
set_verbosity(optuna.logging.CRITICAL)

for n in range(len(boston_df)) :

    study = optuna.create_study(direction='minimize')
    study.optimize(objective_variable(boston_df[boston_df.index==n]), n_trials=100)

    #トライアル結果
    trial = study.best_trial

    #ベストな結果、パラメータの表示
    print(f"Best hyperparameters: {trial.params}")

課題

~~制約付き最適化などは出来ない。~~
- と思ったらすぐ出来るようになるみたいです。
  
  https://t.co/GPorOThFsy 来週にもコレがOptunaに入って多目的最適化+制約付き最適化ができる様になります。
  — まむ (@mamurai1208) December 10, 2020
当然目的変数yを予測するモデルの精度に依存するので、工数をどこまでかけてどのぐらいの精度のモデルを作るのかの判断が難しい。

参考

逆解析とその応用
[機械学習モデルを逆解析する] (https://qiita.com/futakuchi0117/items/92d03c61f715d6b2b7c6)
機械学習モデルの逆解析のススメ
モデルの逆解析をするときのチェックリスト