ポリシリアル相関係数

import numpy as np
from scipy.stats import norm
from scipy.optimize import minimize_scalar
import matplotlib.pyplot as plt
import matplotlib_fontja


def normalize_ordinal(x: np.ndarray[int]) -> np.ndarray[int]:
    """Normalize ordinal variable to be integer-coded starting from 0."""
    unique_values = np.unique(x)
    value_to_code = {value: code for code, value in enumerate(unique_values)}
    return np.vectorize(value_to_code.get)(x)

def polyserial_correlation(x: np.ndarray, y: np.ndarray, plot: bool = False) -> float:
    """
    Estimate the polyserial correlation coefficient between a continuous variable x
    and an ordinal variable y using maximum likelihood estimation.
    """
    x = np.asarray(x)
    y = np.asarray(y)

    if not np.issubdtype(y.dtype, np.integer):
        raise ValueError("y must be an integer-coded ordinal variable.")

    # Standardize x
    z = (x - np.mean(x)) / np.std(x)

    # Estimate thresholds τ from y
    def estimate_thresholds(y):
        levels = np.sort(np.unique(y))
        thresholds = []
        for level in levels[:-1]:  # exclude top category
            p = np.mean(y <= level)
            thresholds.append(norm.ppf(p))
        return np.concatenate(([-np.inf], thresholds, [np.inf]))

    y = normalize_ordinal(y)
    tau = estimate_thresholds(y)

    def neg_log_likelihood(rho):
        log_likelihood = 0.0
        for i in range(len(z)):
            j = y[i]
            tau_lower = (tau[j] - rho * z[i]) / np.sqrt(1 - rho**2)
            tau_upper = (tau[j + 1] - rho * z[i]) / np.sqrt(1 - rho**2)
            p_i = norm.cdf(tau_upper) - norm.cdf(tau_lower)
            p_i = max(p_i, 1e-6)  # soft clipping
            if np.isnan(p_i):
                continue
            log_likelihood += np.log(p_i)
        return -log_likelihood
    
    result = minimize_scalar(neg_log_likelihood, bounds=(-0.999, 0.999), method='bounded')
    
    if plot:
        # plot for debug
        rho_range = np.linspace(-0.999, 0.999, 200)
        likelihoods = np.array([neg_log_likelihood(rho) for rho in rho_range])
        rho_hat = rho_range[np.argmin(likelihoods)]

        fig, ax = plt.subplots(figsize=[4,3])
        ax.plot(rho_range, likelihoods, color="dimgray")
        ax.set(xlabel=r"$\rho$", ylabel="log likelihood", title="Maximum Likelihood Estimate")
        l = likelihoods[~np.isinf(likelihoods)]
        # y = -(l.max() - l.min()) / 2
        y = np.min(likelihoods)
        ax.text(rho_hat * 1.1, y * 1.5, r"$\hat{\rho}$"+f"={rho_hat:.3f}", color="steelblue")
        ax.axvline(rho_hat, color="steelblue")
        fig.show()

    return result.x

import pandas as pd
from scipy.stats import multivariate_normal
n = 100
mean = [50, 0]
std = [10, 3]
rho = 0.5
cov = rho * std[0] * std[1]
Cov = np.array([
    [std[0]**2, cov],
    [cov, std[1]**2]
])

rho = 0.75
Cov = np.array([[1, rho], [rho, 1]])
X = multivariate_normal.rvs(mean=mean, cov=Cov, size=n, random_state=0)
df = pd.DataFrame(X, columns=["x", "y"])
df["y"], _ = pd.cut(df["y"], bins=3).factorize(sort=True)

polyserial_correlation(df["x"], df["y"], plot=True)

import numpy as np
from scipy.stats import norm
from scipy.special import ndtr  # 標準正規CDFの計算本体(引数検証なし・C実装)

rng = np.random.default_rng(0)
v = rng.normal(size=10_000)

# 方法1: Python ループ + norm.cdf(改善前の polyserial と同じパターン)
%timeit -r 3 -n 1 [norm.cdf(vi) for vi in v]

# 方法2: norm.cdf に配列を一括で渡す(ベクトル化)
%timeit -r 3 norm.cdf(v)

# 方法3: ndtr に配列を一括で渡す(ベクトル化 + 低レベル関数)
%timeit -r 3 ndtr(v)

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(0)


def gen(n, k=5, rho=0.5):
    """真の相関 rho の潜在2変量正規から、順序変数 x, y と連続変数 z を作る"""
    latent = rng.multivariate_normal([0, 0], [[1, rho], [rho, 1]], size=n)
    cuts0 = np.quantile(latent[:, 0], np.linspace(0, 1, k + 1)[1:-1])
    cuts1 = np.quantile(latent[:, 1], np.linspace(0, 1, k + 1)[1:-1])
    x = np.digitize(latent[:, 0], cuts0)  # 順序変数 (0..k-1)
    y = np.digitize(latent[:, 1], cuts1)  # 順序変数 (0..k-1)
    return x, y, latent[:, 0]  # z は x の元になった連続変数


x, y, z = gen(1000)
print("x:", x[:10], "...")
print("y:", y[:10], "...")
print("z:", np.round(z[:10], 2), "...")

def estimate_thresholds(values):
    """周辺比率から閾値 τ を推定(±100 を擬似的な ±∞ として使う)"""
    inf = 100
    _, counts = np.unique(values, return_counts=True)
    cum_p = np.cumsum(counts)[:-1] / values.size  # P(Y <= i)、最上位カテゴリは除く
    return np.concatenate(([-inf], norm.ppf(cum_p), [inf]))


# --- 改善前: 観測ごとの Python ループ ---
def nll_loop(rho, z, y, tau):
    ll = 0.0
    for i in range(len(z)):
        j = y[i]
        lo = (tau[j] - rho * z[i]) / np.sqrt(1 - rho**2)
        hi = (tau[j + 1] - rho * z[i]) / np.sqrt(1 - rho**2)
        p = norm.cdf(hi) - norm.cdf(lo)
        ll += np.log(max(p, 1e-6))
    return -ll


x, y, z_raw = gen(5000)
z = (z_raw - z_raw.mean()) / z_raw.std()
tau = estimate_thresholds(y)

%timeit -r 3 -n 1 nll_loop(0.5, z, y, tau)

# fancy indexing のミニ例
tau_demo = np.array([-100.0, -0.8, 0.3, 100.0])  # 3カテゴリの閾値
y_demo = np.array([0, 2, 1, 1, 0])
print("tau[y]   =", tau_demo[y_demo])      # 各観測の下側閾値 τ_{y_i}
print("tau[y+1] =", tau_demo[y_demo + 1])  # 各観測の上側閾値 τ_{y_i+1}

# --- 改善後: 全観測を1回のベクトル演算で ---
def nll_vec(rho, z, y, tau):
    t_lo = tau[y]      # τ_{y_i}      (長さ n の配列)
    t_hi = tau[y + 1]  # τ_{y_i + 1}  (長さ n の配列)
    s = np.sqrt(1 - rho**2)
    p = ndtr((t_hi - rho * z) / s) - ndtr((t_lo - rho * z) / s)
    return -np.sum(np.log(np.maximum(p, 1e-6)))


%timeit -r 3 nll_vec(0.5, z, y, tau)

# 結果は一致する
print("loop:", nll_loop(0.5, z, y, tau))
print("vec: ", nll_vec(0.5, z, y, tau))

import warnings
from typing import Any, Sequence

import numpy as np
import numpy.typing as npt
from scipy.optimize import minimize_scalar
from scipy.special import ndtr, owens_t
from scipy.stats import norm


def polyserial_fast(x: npt.ArrayLike, y: npt.ArrayLike) -> float:
    x = np.asarray(x)
    y = np.asarray(y)

    z = (x - np.mean(x)) / np.std(x, ddof=0)
    y = normalize_ordinal(y)
    tau = estimate_thresholds(y)
    tau_lower = tau[y]  # τ_{y_i}
    tau_upper = tau[y + 1]  # τ_{y_i + 1}

    def neg_log_likelihood(rho: float) -> float:
        scale = np.sqrt(1 - rho**2)
        p = univariate_cdf((tau_lower - rho * z) / scale, (tau_upper - rho * z) / scale)
        p = np.maximum(p, 1e-6)  # soft clipping
        return -np.sum(np.log(p, where=~np.isnan(p), out=np.zeros_like(p)))

    result = minimize_scalar(neg_log_likelihood, bounds=(-1, 1), method="bounded")
    return float(result.x)


def estimate_thresholds(values: npt.NDArray[Any]) -> npt.NDArray[np.float64]:
    r"""Estimate thresholds from empirical marginal proportions"""
    inf = 100  # to make log-likelihood smooth, use large value instead of np.inf
    _, counts = np.unique(values, return_counts=True)
    cum_p = np.cumsum(counts)[:-1] / values.size  # P(X ≤ i), exclude top category
    thresholds = norm.ppf(cum_p)  # τ_i = Φ⁻¹(P(X ≤ i))
    return np.concatenate(([-inf], thresholds, [inf]))


def normalize_ordinal(x: npt.NDArray[Any]) -> npt.NDArray[np.int_]:
    r"""Normalize ordinal variable to be integer-coded starting from 0."""
    _, codes = np.unique(x, return_inverse=True)
    return codes

def univariate_cdf(
    lower: npt.ArrayLike, upper: npt.ArrayLike
) -> npt.NDArray[np.float64]:
    """Compute the univariate cumulative distribution function (CDF) for a standard normal distribution.
    P(lower < X <= upper) = Φ(upper) - Φ(lower)
    where Φ is the CDF of the standard normal distribution.
    """
    return ndtr(upper) - ndtr(lower)

# 改善前
%timeit -r 3 polyserial_correlation(x, y)

# 改善後
%timeit -r 3 polyserial_fast(x, y)