目录

基于约束的方法 (Constraint-based)

基于评分的方法 (Score-based)

基于函数因果模型的方法 (Functional Causal Models)

基于梯度的方法 (Gradient-based)

因果发现是机器学习中一个重要的研究方向，它旨在从观测数据中推断变量之间的因果关系

基于约束的方法 (Constraint-based)

• 核心思想：利用条件独立性检验来推断因果结构
• 代表算法：PC算法、FCI算法、RFCI算法

## 因果发现算法
# 1. 基于约束的方法——PC算法
from pgmpy.estimators import PC
import pandas as pd
import numpy as np

# 生成数据
data = pd.DataFrame(np.random.randn(1000, 3), columns=['X', 'Y', 'Z'])

# PC算法
est = PC(data)
model = est.estimate(variant="orig", alpha=0.05)  
print("因果边:", model.edges())

基于评分的方法 (Score-based)

• 核心思想：定义评分函数评估图结构，搜索得分最高的图
• 代表算法：GES算法、FGES算法

# 2. 基于评分的方法——GES算法
import numpy as np
import pandas as pd
import networkx as nx
import matplotlib.pyplot as plt
from itertools import permutations
from sklearn.linear_model import LinearRegression

np.random.seed(42)
X = np.random.normal(size=1000)
Y = 0.5 * X + np.random.normal(size=1000)
Z = 0.3 * Y + np.random.normal(size=1000)
data = pd.DataFrame({'X': X, 'Y': Y, 'Z': Z})

variables = data.columns.tolist()  
print("变量列表:", variables)  # ['X', 'Y', 'Z']

def compute_bic(data, parents_dict):
    """计算 BIC 评分（线性高斯模型）"""
    score = 0
    for node in data.columns:
        X = data[parents_dict[node]] if parents_dict[node] else np.zeros((len(data), 1))
        y = data[node]
        model = LinearRegression().fit(X, y)
        residuals = y - model.predict(X)
        n = len(data)
        k = len(parents_dict[node])
        sigma2 = np.var(residuals)
        score += -n * np.log(sigma2) / 2 - k * np.log(n) / 2
    return score

# 穷举所有可能的 DAG（检查 edges 生成）
best_score = -np.inf
best_graph = None

all_edges = list(permutations(variables, 2))
print("所有可能的边组合:", all_edges)  # [('X', 'Y'), ('X', 'Z'), ('Y', 'X'), ('Y', 'Z'), ('Z', 'X'), ('Z', 'Y')]

for u, v in all_edges:
    graph = {node: [] for node in variables}
    graph[v].append(u)  # 添加边 u -> v
    score = compute_bic(data, graph)
    if score > best_score:
        best_score = score
        best_graph = graph

# 可视化
G = nx.DiGraph()
for node in best_graph:
    for parent in best_graph[node]:
        G.add_edge(parent, node)

nx.draw(G, with_labels=True, node_color='lightgreen')
plt.title("GES Algorithm (Simplified)")
plt.savefig('ges_graph.png')
plt.show()

基于函数因果模型的方法 (Functional Causal Models)

• 核心思想：假设数据生成过程的函数形式
• 代表算法：LiNGAM、ANM（加性噪声模型）、PNL（后非线性模型）

# 3. 基于函数因果模型的方法——LiNGAM
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import networkx as nx
from itertools import permutations

np.random.seed(42)
X = np.random.normal(size=1000)
Y = 0.5 * X + np.random.normal(size=1000)
Z = 0.3 * Y + np.random.normal(size=1000)
data = pd.DataFrame({'X': X, 'Y': Y, 'Z': Z})

data = (data - data.mean()) / data.std()

try:
    # 尝试使用 causallearn 的 LiNGAM
    from causallearn.search.FCMBased.lingam import DirectLiNGAM

    model = DirectLiNGAM()
    model.fit(data)
    adj_matrix = model.adjacency_matrix_
    print("LiNGAM 邻接矩阵:\n", adj_matrix)
    # [[0.         0.41099056 0.]
    #  [0.         0.         0.]
    # [0.          0.31428138 0.]]

except ImportError:
    print("未找到 causallearn，改用简化版 LiNGAM（基于回归残差）")

    # 简化版 LiNGAM（基于残差独立性检验）
    adj_matrix = np.zeros((3, 3))  # 初始化邻接矩阵
    variables = data.columns.tolist()

    for i, target in enumerate(variables):
        predictors = [var for var in variables if var != target]
        X_pred = data[predictors].values
        y_target = data[target].values

        # 多元线性回归
        coef = np.linalg.lstsq(X_pred, y_target, rcond=None)[0]
        residuals = y_target - X_pred.dot(coef)

        # 检查残差与预测变量的独立性（简化版：相关系数）
        for j, predictor in enumerate(predictors):
            corr = np.corrcoef(residuals, data[predictor])[0, 1]
            if abs(corr) < 0.05:  # 阈值可调整
                adj_matrix[variables.index(predictor), i] = 1  # predictor -> target

# 可视化
G = nx.DiGraph(adj_matrix)
pos = nx.spring_layout(G)
nx.draw(G, pos, with_labels=True, node_size=1000, arrowsize=20,
        node_color='lightgreen', edge_color='gray')
plt.title('LiNGAM Causal Graph (Direct)' if 'model' in locals() else 'Simplified LiNGAM')
plt.show()

基于梯度的方法 (Gradient-based)

• 核心思想：使用神经网络和梯度下降学习因果结构
• 代表算法：DAG-GNN、NOTEARS

# 4. 基于梯度的方法——NOTEARS
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt

# 生成数据
np.random.seed(42)
X = np.random.normal(size=1000)
Y = 0.5 * X + np.random.normal(size=1000)
Z = 0.3 * Y + np.random.normal(size=1000)
data = np.column_stack([X, Y, Z])
data = (data - data.mean(axis=0)) / data.std(axis=0)  # 标准化


# 自定义简化版NOTEARS（梯度下降优化）
def notears_simple(X, lambda1=0.1, max_iter=100):
    n_features = X.shape[1]
    W = np.zeros((n_features, n_features))  # 初始化邻接矩阵

    for _ in range(max_iter):
        # 计算梯度（最小二乘损失 + 无环约束）
        grad = -X.T @ (X - X @ W) / len(X) + lambda1 * np.sign(W)

        # 投影梯度更新（确保无环）
        W -= 0.01 * grad  # 学习率
        W = np.clip(W, -1, 1)  # 限制权重范围

    return (W != 0).astype(int)  # 二值化邻接矩阵


# 运行并可视化
W_est = notears_simple(data)
labels = ['X', 'Y', 'Z']
G = nx.DiGraph(W_est, labels=labels)
nx.draw(G, with_labels=True, node_color='lightblue', arrowsize=20)
plt.title('Simplified NOTEARS Result')
plt.show()