用多尺度图相关区分线性与螺旋依赖

散点图中一条明显的曲线,不一定对应很大的线性相关系数。SciPy 的多尺度图相关检验(Multiscale Graph Correlation,MGC)可以用于高维、非线性数据的独立性检验;它还提供局部相关图,帮助观察哪些邻域尺度对当前数据的关系更有信息。

本文依据 SciPy 1.18.0 官方教程完整整理线性和螺旋两个例子,并用对应 API 文档核对参数。原作者为 SciPy 文档贡献者,页面无个人署名,也未给出首次发表日期。代码仅作静态审查,本文没有运行模拟或置换检验。

MGC 工作流程:成对观测生成距离,遍历 X 和 Y 的邻域尺度得到局部相关图,然后计算统计量与置换 p 值;线性和螺旋是两种示例关系。
原创概念示意图。曲线表示关系形状,不是本次随机模拟、局部相关热图或实测检验结果。

先明确三个输出分别回答什么

MGC 为 X 和 Y 分别考虑近邻数,并把一对近邻数 (k, l) 称为一个尺度。不同尺度下的局部相关值组成二维 mgc_map。算法在这些尺度上寻找经过平滑处理的最优相关,因此不是简单选择一个孤立的高值像素。

multiscale_graphcorr(x, y) 的结果包括统计量、置换检验 p 值以及 mgc_dict。统计量描述本样本中检出的关联,p 值用于判断数据与独立性假设是否相容;mgc_dict["opt_scale"] 则是估计的最优 X、Y 邻域尺度。局部图提供解释线索,但不是对真实生成机制的证明。

对于通常的独立性检验,两组输入应为成对观测,样本数相同;多维输入可使用 (n, p) 与 (n, q)。本教程用两个长度为 100 的一维数组。不要把未对齐的两组样本直接解释为成对独立性检验;该 API 还包含两样本检验行为,实际输入形状应按所用模式核对。

准备绘图函数

原文先导入 NumPy、Matplotlib 和 MGC,然后定义一个绘图函数:可以只画散点图,也可以只画局部相关图。下面保留这两个用途,同时做了明确的编辑修订:使用局部图实际维度,不把范围写死为 100;将矩阵转置,使横轴对应 X 的邻域数、纵轴对应 Y;用从 1 开始的尺度坐标显示热图并标记 opt_scale。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multiscale_graphcorr

plt.style.use("classic")

def mgc_plot(x, y, sim_name, mgc_dict=None,
             only_viz=False, only_mgc=False):
    if not only_mgc:
        fig, ax = plt.subplots(figsize=(8, 8))
        ax.set_title(f"{sim_name} simulation")
        ax.scatter(x, y)
        ax.set(xlabel="X", ylabel="Y")
        ax.axis("equal")
        plt.show()

    if not only_viz:
        if mgc_dict is None:
            raise ValueError("绘制局部相关图需要 mgc_dict")
        local_map = np.asarray(mgc_dict["mgc_map"])
        nx, ny = local_map.shape
        fig, ax = plt.subplots(figsize=(8, 8))
        im = ax.imshow(
            local_map.T,
            origin="lower",
            extent=(0.5, nx + 0.5, 0.5, ny + 0.5),
            cmap="YlGnBu",
            aspect="auto",
        )
        fig.colorbar(im, ax=ax, label="Local correlation")
        k, l = mgc_dict["opt_scale"]
        ax.scatter(k, l, marker="X", s=160, color="red")
        ax.set(
            title="Local correlation map",
            xlabel="Neighbors for X",
            ylabel="Neighbors for Y",
        )
        plt.show()

散点图是在原始观测空间中观察 X 与 Y;热图则是在“邻域尺度”空间中观察相关值。两幅图坐标轴的含义不同。红色叉号表示算法返回的最优尺度,而不是数据中的某个观测点。

例一:有噪声的线性关系

原文在 −1 到 1 之间等距取 100 个 X 值,再给 y = x 加上随机扰动。扰动为 0.3 * rng.random(...),因此来自非负区间,并非均值为零的高斯噪声。

# 编辑补充固定随机种子,便于读者在自己的环境复核。
rng = np.random.default_rng(20261005)
x = np.linspace(-1, 1, num=100)
y = x + 0.3 * rng.random(x.size)

mgc_plot(x, y, "Linear", only_viz=True)

stat, pvalue, mgc_dict = multiscale_graphcorr(
    x, y, reps=1000, workers=1, random_state=20261005
)
print("MGC test statistic:", stat)
print("P-value:", repr(pvalue))
print("Optimal scale:", mgc_dict["opt_scale"])
mgc_plot(x, y, "Linear", mgc_dict, only_mgc=True)

原教程展示的统计量四舍五入到一位小数后为 1.0,p 值显示为 0.0。这里必须区分显示结果与实际数值:四舍五入后的 0.0 不是 p 值恰好为零。上面的代码保留浮点输出,并显式给出默认数量级的 1000 次置换。

原文的线性局部相关图在全局尺度取得最优值。直观上,对这组近似直线的数据,考虑更多邻居有助于捕捉整体关系。这个解释针对原教程所展示的样本,不应改写为“所有线性数据都必定得到同一最优尺度”。本文没有复制原结果作为本地运行输出;固定种子后的数值也可能与原图不同。

例二:同一个随机参数生成螺旋

第二组数据从 0 到 5 的均匀分布抽取 100 个参数值,再用余弦和正弦生成螺旋坐标,并为 Y 加上非负噪声。X 与 Y 由同一个参数生成,因此即使关系不是直线,也具有共同结构。

unif = np.array(rng.uniform(0, 5, size=100))
x = unif * np.cos(np.pi * unif)
y = unif * np.sin(np.pi * unif) + 0.4 * rng.random(x.size)

mgc_plot(x, y, "Spiral", only_viz=True)

stat, pvalue, mgc_dict = multiscale_graphcorr(
    x, y, reps=1000, workers=1, random_state=20261006
)
print("MGC test statistic:", stat)
print("P-value:", repr(pvalue))
print("Optimal scale:", mgc_dict["opt_scale"])
mgc_plot(x, y, "Spiral", mgc_dict, only_mgc=True)

原文对这一随机例子展示的统计量约为 0.2,p 值同样被保留一位小数显示为 0.0。它的局部相关图在局部尺度取得最优值,说明在这一例子中,小范围邻域比把全部点视为邻居更有帮助。不要机械比较 1.0 与 0.2 就得出“螺旋没有重要关系”;统计量、置换 p 值和数据几何形状需要一起阅读。

怎样解释,怎样保留不确定性

MGC 的 p 值通过置换估计:打乱样本配对形成独立性假设下的参考分布,再与观测统计量比较。有限次重采样会带来精度限制与随机波动。分析记录应保留完整 p 值、样本数、置换次数、随机种子和软件版本;如果小 p 值对决策关键,可在明确计算成本后增加置换次数,而不应仅为了显示好看而把它写成零。

低 p 值支持“数据存在统计依赖”的判断,不能推出因果关系。最优全局或局部尺度,也不是自动给数据贴上“线性”“螺旋”的分类标签。噪声、样本量、尺度选择和输入预处理都会影响结果。这个教程的价值,是把可见的两类模拟关系与 MGC 输出连接起来,帮助理解局部相关图。

安全审查方面,这些示例仅在内存中生成数组并绘图,没有外部下载、系统命令、硬编码凭据或文件删除操作。置换检验可能消耗较多计算资源,示例保留 workers=1;本文只做静态核对,未运行并行进程或宣称测试通过。

原始版权与许可

SciPy 示例的 BSD 3-Clause 许可及原始版权声明保留如下;中文整理与代码修订已在正文标明。 官方许可来源

Copyright (c) 2001-2002 Enthought, Inc. 2003, SciPy Developers.
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:

1. Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.

2. Redistributions in binary form must reproduce the above
   copyright notice, this list of conditions and the following
   disclaimer in the documentation and/or other materials provided
   with the distribution.

3. Neither the name of the copyright holder nor the names of its
   contributors may be used to endorse or promote products derived
   from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
© 版权声明
THE END
喜欢就支持一下吧
点赞0 分享
评论 抢沙发

请登录后发表评论

    暂无评论内容