机器学习中有哪些涉及统计因果推断的算法？

## 准备工作
import warnings
warnings.simplefilter('ignore')
import os
import urllib.request
import numpy as np
import pandas as pd
from networkx.drawing.nx_pydot import to_pydot
from IPython.display import Image, display
## 导入机器学习包的相关工具
from sklearn.preprocessing import PolynomialFeatures
from sklearn.ensemble import GradientBoostingRegressor
## 导入 EconML 包的相关工具
from econml.dml import LinearDML, CausalForestDML
from econml.cate_interpreter import SingleTreeCateInterpreter, SingleTreePolicyInterpreter
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format = 'retina'
## 数据导入
file_url = "https://msalicedatapublic.blob.core.windows.net/datasets/Pricing/pricing_sample.csv"
train_data = pd.read_csv(file_url)
train_data.head()
train_data.describe()

## 变量转换
## 结果变量 target
train_data["log_demand"] = np.log(train_data




    
["demand"])
train_data["log_price"] = np.log(train_data["price"])
Y = train_data["log_demand"].values ## from Series to np array
## 干预变量 treatment
T = train_data["log_price"].values
## 关键自变量
X = train_data[["income"]].values  
## 混淆变量 - 控制变量
confounder_names = ["account_age", "age", "avg_hours", "days_visited", "friends_count", "has_membership", "is_US", "songs_purchased"]
W = train_data[confounder_names].values
## 构造关键自变量的测试数据
X_test = np.linspace(0, 5, 100).reshape(-1, 1)
X_test_data = pd.DataFrame(X_test, columns=["income"])

## 初始化 EconML 估计器
est = LinearDML(model_y=GradientBoostingRegressor(), model_t=GradientBoostingRegressor(),
              featurizer=PolynomialFeatures(degree=2, include_bias=False))
## 通过 DoWhy 进行拟合
est_dw = est.dowhy.fit(Y, T, X=X, W=W, outcome_names=["log_demand"], treatment_names=["log_price"], feature_names=["income"],
               confounder_names=confounder_names, inference="statsmodels")
## 画出因果图
try:
    # Try pretty printing the graph. Requires pydot and pygraphviz
    display(
        Image(to_pydot(est_dw._graph._graph).create_png())
except:
    # Fall back on default graph view
    est_dw.view_model() 

## 输出估计结果
identified_estimand = est_dw.identified_estimand_
print(identified_estimand)
[Out: ]
Estimand type: nonparametric-ate
### Estimand : 1
Estimand name: backdoor
Estimand expression:
────────────(Expectation(log_demand|friends_count,income,avg_hours,account_age
d[log_price]                                                                  
,is_US,days_visited,age,has_membership,songs_purchased))
Estimand assumption 1, Unconfoundedness: If U→{log_price} and U→log_demand then P(log_demand|log_price,friends_count,income,avg_hours,account_age,is_US,days_visited,age,has_membership,songs_purchased,U) = P(log_demand|log_price,friends_count,income,avg_hours,account_age,is_US,days_visited,age,has_membership,songs_purchased)
### Estimand : 2
Estimand name: iv
No such variable found!
### Estimand : 3
Estimand name: frontdoor
No such variable found!

## 基于真实的数据生成过程，估计因果效应
def gamma_fn(X):
    return -3 - 14 * (X["income"] < 1)
def beta_fn(X):
    return 20 + 0.5 * (X["avg_hours"]) + 5 * (X["days_visited"] > 4)
def demand_fn(data, T):
    Y = gamma_fn(data) * T + beta_fn(data)
    return Y
def true_te(x, n, stats):
    if x < 1:
        subdata = train_data[train_data["income"] < 1].sample(n=n, replace=True)
    else:
        subdata = train_data[train_data["income"] >= 1].sample(n=n, replace=True)
    te_array = subdata["price"] * gamma_fn(subdata) / (subdata["demand"])
    if stats == "mean":
        return np.mean(te_array)
    elif stats == "median":
        return np.median(te_array)
    elif isinstance(stats, int):
        return np.percentile(te_array, stats)
## 估计结果和置信区间
truth_te_estimate = np.apply_along_axis(true_te, 1, X_test, 1000, "mean")  
truth_te_upper = np.apply_along_axis(




    
true_te, 1, X_test, 1000, 95)  
truth_te_lower = np.apply_along_axis(true_te, 1, X_test, 1000, 5)  

plt.figure(figsize=(16, 10))
plt.plot(X_test.flatten(), te_pred, label="Sales Elasticity Prediction")
plt.plot(X_test.flatten(), truth_te_estimate, "--", label="True Elasticity")
plt.fill_between(
    X_test.flatten(),
    te_pred_interval[0].flatten(),
    te_pred_interval[1].flatten(),
    alpha=0.2,
    label="95% Confidence Interval",
plt.fill_between(
    X_test.flatten(),
    truth_te_lower,
    truth_te_upper,
    alpha=0.2,
    label="True Elasticity Range",
plt.xlabel("Income")
plt.ylabel("Songs Sales Elasticity")
plt.title("Songs Sales Elasticity vs Income")
plt.legend(loc="lower right")

## 显示X和截距项的估计结果信息
est_dw.summary(feature_names=["income"])

## 初始化
est_nonparam = CausalForestDML(model_y=GradientBoostingRegressor(), model_t=GradientBoostingRegressor())
## 通过 DoWhy 进行拟合
est_nonparam_dw = est_nonparam.dowhy.fit(Y, T, X=X, W=W, outcome_names=["log_demand"], treatment_names=["log_price"],
                                         feature_names=["income"], confounder_names=confounder_names, inference="blb")
## 计算干预效应及其置信区间
te_pred = est_nonparam_dw.effect(X_test).flatten()
te_pred_interval = est_nonparam_dw.effect_interval(X_test)
## 比较非参估计的结果和真实结果
plt.figure(figsize=(16, 10))
plt.plot(X_test.flatten(), te_pred, label="Sales Elasticity Prediction")
plt.plot(X_test.flatten(), truth_te_estimate, "--", label="True Elasticity")
plt.fill_between(
    X_test.flatten(),
    te_pred_interval[0].flatten(),
    te_pred_interval[1].flatten(),
    alpha=0.2,
    label="95% Confidence Interval",
plt.fill_between(
    X_test.flatten(),
    truth_te_lower,
    truth_te_upper,
    alpha=0.2,
    label="True Elasticity Range",
plt.xlabel("Income")
plt.ylabel("Songs Sales Elasticity")
plt.title("Songs Sales Elasticity vs Income")
plt.legend(loc="lower right")

res_random = est_nonparam_dw.refute_estimate(method_name="random_common_cause")
print(res_random)
[Out: ]
Refute: Add a Random Common Cause
Estimated effect:-0.9561140448923923
New effect:-0.963688374616349

res_unobserved = est_nonparam_dw.refute_estimate(
    method_name="add_unobserved_common_cause",
    confounders_effect_on_treatment="linear",
    confounders_effect_on_outcome="linear",
    effect_strength_on_treatment=0.1,
    effect_strength_on_outcome=0.1,
print(res_unobserved)
[Out: ]
Refute: Add an Unobserved Common Cause
Estimated effect:-0.9561140448923923
New effect:0.1932125423770007

res_placebo = est_nonparam_dw.refute_estimate(
    method_name="placebo_treatment_refuter", placebo_type="permute", 
    num_simulations=3
print(res_placebo)
[Out: ]
Refute: Use a Placebo Treatment
Estimated effect:-0.9561140448923923
New effect:-0.00769355893605958
p value:0.21924855198612314

res_subset = est_nonparam_dw.refute_estimate(
    method_name="data_subset_refuter", subset_fraction=0.8, 
    num_simulations=3)
print(res_subset)
Refute: Use a subset of data
Estimated effect:-0.9561140448923923
New effect:-0.957908388178223
p value:0.022684332959365117

intrp = SingleTreeCateInterpreter(include_model_uncertainty=True, max_depth=2, min_samples_leaf=10)
intrp.interpret(est_nonparam_dw, X_test)
plt.figure(figsize=(25, 5))
intrp.plot(feature_names=["income"], fontsize=12)

intrp = SingleTreePolicyInterpreter(risk_level=0.05, max_depth=2, min_samples_leaf=1, min_impurity_decrease=0.001)
intrp.interpret(est_nonparam_dw, X_test, sample_treatment_costs=-1)
plt.figure(figsize=(25, 5))
intrp.plot(feature_names=["income"], treatment_names=["Discount", "No-Discount"], fontsize=12)

## 定义折扣策略收益函数
def revenue_fn(data, discount_level1, discount_level2, baseline_T, policy):
    policy_price = baseline_T * (1 - discount_level1) * policy + baseline_T * (1 - discount_level2) * (1 - policy)
    demand = demand_fn(data, policy_price)
    rev = demand * policy_price
    return rev
policy_dic = {}
## 本文的策略
policy = intrp.treat(X)
policy_dic["Our Policy"] = np.mean(revenue_fn(train_data, 0, 0.1, 1, policy))
## 原始策略
policy_dic["Previous Strategy"] = np.mean(train_data["price"] * train_data["demand"])
## 对所有人都打折
policy_dic["Give Everyone Discount"] = np.mean(revenue_fn(train_data, 0.1, 0, 1, np.ones(len(X))))
## 对所有人都不打折
policy_dic["Give No One Discount"] = np.mean(revenue_fn(train_data, 0, 0.1, 1, np.ones(len(X))))
## 对穷人打折，对富人涨价