順序尺度の相関係数

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib_fontja
from scipy.stats import multivariate_normal, pearsonr, spearmanr, kendalltau
from ordinalcorr import polychoric, polyserial

def gen_bivariate_normal(rho=0.5, size=1000, seed=0):
    # 2変量標準正規分布によるデータ生成
    mean = [0, 0]
    std = [1, 1]
    cov = rho * std[0] * std[1]
    Cov = np.array([[std[0] ** 2, cov], [cov, std[1] ** 2]])
    X = multivariate_normal.rvs(mean=mean, cov=Cov, size=size, random_state=seed)
    df = pd.DataFrame(X, columns=["x", "y"])
    return df

def plot_data(ax, continuous_data, rho, n):
    sns.scatterplot(data=continuous_data, x="x", y="y", ax=ax)
    ax.set(title=f"2変量標準正規分布に従うデータ（ρ={rho}, n={n}）")

def plot_polyserial(ax, continuous_data, rho):
    # 離散化のカテゴリ数を変えつつ相関係数を算出
    results = []
    for k in range(2, 11):
        x = continuous_data["x"]
        y, _ = pd.cut(continuous_data["y"], bins=k).factorize(sort=True)
        results += [
            dict(method="Pearson", value=pearsonr(x, y).statistic, k=k),
            dict(method="Polyserial", value=polyserial(x, y), k=k),
        ]
    results = pd.DataFrame(results)

    # データをグラフにplot
    sns.lineplot(x="k", y="value", data=results, hue="method", marker="o", ax=ax)
    ax.axhline(rho, label=f"真値 (ρ={rho})", color="black")
    ax.set(
        xlabel="カテゴリ数",
        ylabel=r"推定された相関係数 $\hat{\rho}$",
        title="連続変数xと離散化されたyの相関係数"
    )
    ax.legend()

def plot_polychoric(ax, continuous_data, rho):
    # 離散化のカテゴリ数を変えつつ相関係数を算出
    results = []
    for k in range(2, 11):
        x, _ = pd.cut(continuous_data["x"], bins=k).factorize(sort=True)
        y, _ = pd.cut(continuous_data["y"], bins=k).factorize(sort=True)
        results += [
            dict(method="Pearson", value=pearsonr(x, y).statistic, k=k),
            dict(method="Polychoric", value=polychoric(x, y), k=k),
            dict(method="Spearman", value=spearmanr(x, y).statistic, k=k),
            dict(method="Kendall", value=kendalltau(x, y).statistic, k=k),
        ]
    results = pd.DataFrame(results)

    # データをグラフにplot
    sns.lineplot(x="k", y="value", data=results, hue="method", marker="o", ax=ax)
    ax.axhline(rho, label=f"真値 (ρ={rho})", color="black")
    ax.set(
        xlabel="カテゴリ数",
        ylabel=r"推定された相関係数 $\hat{\rho}$",
        title="離散化された x, y の相関係数"
    )
    ax.legend()

# 実験
n = 1000
fig, axes = plt.subplots(figsize=[12, 9], nrows=3, ncols=3, tight_layout=True)
rhos = [-0.5, 0, 0.5]
for i, rho in enumerate(rhos):
    continuous_data = gen_bivariate_normal(rho=rho, size=n)
    plot_data(axes[i, 0], continuous_data, rho, n)
    plot_polyserial(axes[i, 1], continuous_data, rho)
    plot_polychoric(axes[i, 2], continuous_data, rho)