大規模な不均衡データに対するロジスティック回帰（後編）

離散時間 MCMC から連続時間 MCMC へ

Bayesian
Computation
Python
MCMC
Author

司馬博文

Published

7/18/2024

Modified

7/23/2024

概要
ロジットモデルやプロビットモデルの事後分布からのサンプリングには，その混合構造を利用したデータ拡張による Gibbs サンプラーが考案されている．しかし，このような Gibbs サンプラーは不明な理由で極めて収束が遅くなることがよく見られ，そのうちの１つのパターンが 大規模な不均衡データ である．前編ではこの現象がなぜ起こるかに関して考察した．ここでは代替手法として Zig-Zag サンプラーがうまくいくことをみる．

前稿はこちら：

1 Zig-Zag サンプラーによる解決

1.1 はじめに

説明変数 \(X\in\mathcal{L}(\Omega;\mathbb{R}^p)\) と係数 \(\xi\in\mathcal{L}(\Omega;\mathbb{R}^p)\) に関するロジスティックモデル

\[ \operatorname{P}[Y=1\,|\,X,\xi]=g^{-1}(X^\top\xi)=\frac{\exp(X^\top\xi)}{1+\exp(X^\top\xi)}\tag{1} \]

において，\(\{y^i\}\) または \(\{x^i_j\}\) が不均衡であった場合は，前節で見たように Gibbs サンプラーによる方法ではスケールしない．

よく調整された Metropolis-Hastings 法を用いることが解決の１つとしてあり得るが，やはり本当に大きな \(n,p\) に関してスケールしないことが問題である．

1.2 事後分布の特徴

実は，ロジスティック回帰ではポテンシャル \(U\) の勾配が有界になり，簡単なサブサンプリングが可能である．

事後分布の負の対数密度は \[\begin{align*} U(\xi)&:=-\log p_0(\xi)-\sum_{i=1}^n\log\left(\frac{\exp(y^i(x^i)^\top\xi)}{1+\exp((x^i)^\top\xi)}\right)\\ &=:U_0(\xi)+U_1(\xi) \end{align*}\] と表せる．

\[ U_1(\xi)=\frac{1}{n}\sum_{j=1}^nU^j_1(\xi) \] \[ U_1^j(\xi)=-n\log\left(\frac{\exp\left(y^j(x^j)^\top\xi\right)}{1+\exp\left((x^j)^\top\xi\right)}\right) \] であるから， \[ \partial_iU^j_1(\xi)=n\frac{x^j_i\exp\left((x^j)^\top\xi\right)}{1+\exp\left((x^j)^\top\xi\right)}-ny^jx^j_i<nx^j_i(1-y^j) \] を得る．

1.3 Poisson 剪定

特に各 \(\partial_iU_1(\xi)\) のバウンドとして \[ \lvert\partial_iU^j_1(\xi)\rvert\le n\max_{j\in[n]}\lvert x^j_i\rvert\qquad i\in[d] \tag{1}\] を得るから， \[ \lambda_i(\xi,\theta)=\biggr(\theta_i\partial_iU(\xi)\biggl)_+\le\biggr(n\theta_i\max_{j\in[n]}\lvert x^j_i\rvert\biggl)_+ \] を元に剪定を実行できる．

サブサンプリングなしの Zig-Zag 過程のシミュレーションをする関数 ZZ1d() を定義 
using ZigZagBoomerang
using Distributions
using Random
using LinearAlgebra
using Statistics  # just for sure
using StatsFuns

"""
    ∇U(i,j,ξ,x,y)
        i ∈ [d]: 次元を表すインデックス
        j ∈ [n]: サンプル番号を表すインデックス
        ξ: パラメータ空間 R^d 上の位置
        他，観測 (x,y) を引数にとる．
    この関数を実装する際，log の中身をそのまま計算しようとすると大変大きくなり，数値的に不安定になる（除算の後は 1 近くになるはずだが，Inf になってしまう）
"""
∇U(i::Int64, j::Int64, ξ, x::Matrix{Float64}, y::Vector{Float64}) = length(y) * x[i,j] * (logistic(dot(x[:,j],ξ)) - y[j])

"""
    ∇U(i,ξ,x,y)：∇U(i,j,ξ,x,y) を全データ j ∈ [n] について足し合わせたもの
        i ∈ [d]: 次元を表すインデックス
        ξ: パラメータ空間 R^d 上の位置
        他，観測 (x,y) を引数にとる．
"""
function ∇U(i::Int64, ξ, x::Matrix{Float64}, y::Vector{Float64})
    n = length(y)
    U_list = []
    for j in 1:n
        push!(U_list, ∇U(i, j, ξ, x, y))
    end
    return mean(U_list)
end

function  ∇U(ξ, x::Matrix{Float64}, y::Vector{Float64})  # 1次元の場合のショートカット
    return ∇U(1, ξ, x, y)
end

pos(x) = max(zero(x), x)

"""
    λ(i, ξ, θ, ∇U, x, y)：第 i ∈ [d] 次元のレート関数
        i ∈ [d]: 次元を表すインデックス
        (ξ,θ): E 上の座標
        ∇U
        (x,y): 観測
"""
λ(i::Int64, ξ, θ, ∇U, x, y) = pos(θ[i] * ∇U(i, ξ, x, y))
λ(ξ, θ, ∇U, x, y) = pos(θ * ∇U(ξ, x, y))  # 1次元の場合のショートカット

"""
    λ(τ, a, b)：代理レート関数の時刻 τ における値
        τ: 時間
        a,b: 1次関数の係数
"""
λ_bar(τ, a, b) = pos(a + b*τ)

"""
`x`: current location, `θ`: current velocity, `t`: current time,
"""
function move_forward(τ, t, ξ, θ, ::ZigZag1d)
    τ + t, ξ + θ*τ , θ
end

"""
    ZZ1d(∇U, ξ, θ, T, x, y, Flow; rng=Random.GLOBAL_RNG, ab=ab_Global)：ZigZag sampler without subsampling
        `∇U`: gradient of the negative log-density
        `(ξ,θ)`: initial state
        `T`: Time Horizon
        `(x,y)`: observation
        `Flow`: continuous dynamics

        `a+bt`: computational bound for intensity m(t)

        `num`: ポアソン時刻に到着した回数
        `acc`: 受容回数．`acc/num` は acceptance rate
"""
function ZZ1d(∇U, ξ, θ, T::Float64, x::Matrix{Float64}, y::Vector{Float64}, Flow::ZigZagBoomerang.ContinuousDynamics; rng=Random.GLOBAL_RNG, ab=ab_Global)
    t = zero(T)
    Ξ = [(t, ξ, θ)]
    num = acc = 0
    epoch_list = [num]
    a, b = ab(ξ, θ, x, y, Flow)
    t′ =  t + poisson_time(a, b, rand())  # イベントは a,b が定める affine proxy に従って生成する

    while t < T
        τ = t′ - t
        t, ξ, θ = move_forward(τ, t, ξ, θ, Flow)
        l, lb = λ(ξ, θ, ∇U, x, y), λ_bar(τ, a, b)  # λ が真のレート, λ_bar が affine proxy
        num += 1
        if rand()*lb < l
            acc += 1
            if l > lb + 0.01
                println(l-lb)
                println(l)
            end
            θ = -θ
            push!(Ξ, (t, ξ, θ))
            push!(epoch_list, num)
        end
        a, b = ab(ξ, θ, x, y, Flow)
        t′ = t + poisson_time(a, b, rand())
    end

    return Ξ, epoch_list, acc/num
end
a_Global(ξ, θ, x, y) = length(y) * maximum(abs.(vec(x)))
b_Global(ξ, θ, x, y) = 0

ab_Global(ξ, θ, x, y, ::ZigZag1d) = (a_Global(ξ, θ, x, y), b_Global(ξ, θ, x, y))

そこで，さらにタイトに，次の affine 関数による評価を考えることになる： \[ \overline{m}_i(t):=a_i+b_it,\qquad i\in[d], \] \[ a_i:=(\theta_i\partial_iU(\xi_*))_++C_i\lvert\xi-\xi_*\rvert \] \[ b_i:=C_i\sqrt{d},\qquad C_i:=\frac{n}{4}\max_{j\in[n]}\lvert x^j_i\rvert\lvert x^j\rvert. \]

観測を生成 
using StatsFuns
using Distributions

ξ0 = [1.0] # True value
n_list = [10, 100, 1000]  # 実験で用いるサンプルサイズの列

Σ = [2]
x = rand(MvNormal(ξ0, Σ), n_list[end])
y = rand.(Bernoulli.(logistic.(ξ0*x)))  # BitVector になってしまう
y = Float64.(vec(y))  # Vector{Float64} に変換
using Statistics
using LinearAlgebra

"""
    U(ξ, x, y)：ポテンシャル関数
        ξ: パラメータ空間上の点
        (x,y): 観測
"""
function U(ξ, x, y)
    n = length(y)
    U_list = []
    for j in 1:n
        push!(U_list, U(j, ξ, x, y))
    end
    return mean(U_list)
end
function U(j, ξ, x, y)
    n = length(y)
    product = dot(x[:,j],ξ)
    return -n * log(exp(y[j] * product) / (1 + exp(product)))
end

using Optim

result = optimize(ξ -> U(ξ, x, y), [0.0], LBFGS())
ξ_star = Optim.minimizer(result)

function C(ξ, θ, x, y)
    n = length(y)
    max_value = maximum(x.^2)
    return n * max_value / 4
end

a_Affine(ξ, θ, x, y) = pos(θ * ∇U(ξ_star,x,y)) + C(ξ, θ, x, y) * abs(ξ - ξ_star[1])
b_Affine(ξ, θ, x, y) = C(ξ, θ, x, y)

# computational bounds for intensity m(t)
ab_Affine(ξ, θ, x, y, ::ZigZag1d) = (a_Affine(ξ, θ, x, y), b_Affine(ξ, θ, x, y))

1.4 サブサンプリング

Poisson 過程の強度関数を確率化し，\(K\sim\mathrm{U}([n])\) に対して \[ m_i^K(t):=\biggr(\theta_i E^K_i(x+\theta t)\biggl)_+ \] \[ E^K_i(x):=\partial_iU(\xi_*)+\partial_iU^K(\xi)-\partial_iU^K(\xi_*) \] としても，引き続き同様の上界を持つ．

ここで，\(m_i^K\) の評価は \(m_i\) より \(n\) 倍軽量になっていることに注意．

サブサンプリングありの Zig-Zag 過程のシミュレーションをする関数 ZZ1d_SS() と ZZ1d_CV() を定義 
using ZigZagBoomerang
using Distributions
using Random
using LinearAlgebra
using Statistics  # just for sure
using StatsFuns

function λj_Global(j::Int64, ξ, θ, ∇U, x, y)
    Eʲ = ∇U(1, j, ξ, x, y)
    return pos(θ * Eʲ)
end

function ZZ1d_SS(∇U, ξ, θ, T::Float64, x::Matrix{Float64}, y::Vector{Float64}, Flow::ZigZagBoomerang.ContinuousDynamics; rng=Random.GLOBAL_RNG, ab=ab_Global)
    t = zero(T)
    Ξ = [(t, ξ, θ)]
    num = acc = 0
    epoch_list = [num]
    a, b = ab(ξ, θ, x, y, Flow)
    t′ =  t + poisson_time(a, b, rand())  # イベントは a,b が定める affine proxy に従って生成する

    while t < T
        τ = t′ - t
        t, ξ, θ = move_forward(τ, t, ξ, θ, Flow)
        j = rand(1:length(y))
        l, lb = λj_Global(j, ξ, θ, ∇U, x, y), λ_bar(τ, a, b)  # λ が真のレート, λ_bar が affine proxy
        num += 1
        if rand()*lb < l
            acc += 1
            if l > lb + 0.01
                println(l-lb)
            end
            θ = -θ
            push!(Ξ, (t, ξ, θ))
            push!(epoch_list, num)
        end
        a, b = ab(ξ, θ, x, y, Flow)
        t′ = t + poisson_time(a, b, rand())
    end

    return Ξ, epoch_list, acc/num
end

function λj(j::Int64, ξ, θ, ∇U, x, y)
    Eʲ = ∇U(ξ_star, x, y) + ∇U(1, j, ξ, x, y) - ∇U(1, j, ξ_star, x, y)
    return pos(θ * Eʲ)
end

function ZZ1d_CV(∇U, ξ, θ, T::Float64, x::Matrix{Float64}, y::Vector{Float64}, Flow::ZigZagBoomerang.ContinuousDynamics; rng=Random.GLOBAL_RNG, ab=ab_Affine)
    t = zero(T)
    Ξ = [(t, ξ, θ)]
    num = acc = 0
    epoch_list = [num]
    a, b = ab(ξ, θ, x, y, Flow)
    t′ =  t + poisson_time(a, b, rand())  # イベントは a,b が定める affine proxy に従って生成する

    while t < T
        τ = t′ - t
        t, ξ, θ = move_forward(τ, t, ξ, θ, Flow)
        j = rand(1:length(y))
        l, lb = λj(j, ξ, θ, ∇U, x, y), λ_bar(τ, a, b)  # λ が真のレート, λ_bar が affine proxy
        num += 1
        if rand()*lb < l
            acc += 1
            if l > lb + 0.01
                println(l-lb)
            end
            θ = -θ
            push!(Ξ, (t, ξ, θ))
            push!(epoch_list, num)
        end
        a, b = ab(ξ, θ, x, y, Flow)
        t′ = t + poisson_time(a, b, rand())
    end

    return Ξ, epoch_list, acc/num
end

1.5 数値実験による性能比較

サブサンプリングなしの実験を実行する関数 experiment_ZZ() を定義 
using Statistics

function ESS(samples::Vector{Float64}, T, dt)
    B = T / dt
    V = (dt / T) * sum(samples.^2) - ((dt / T) * sum(samples))^2
    Y = samples .* sqrt(T / B)
    ESS = T * V / var(Y)
    return ESS
end

function getESSperEpoch(ab, T ,dt, x, y; ξ0=0.0, θ0=1.0)
    trace, epochs, acc = ZZ1d(∇U, ξ0, θ0, T, x, y, ZigZag1d(); ab=ab)
    traj = discretize(trace, ZigZag1d(), dt)
    return ESS(traj.x, T, dt) / epochs[end]
end

N = 10
T = 500.0
dt = 0.1

function experiment_ZZ(N, T, dt; ξ0=0.0, θ0=1.0, n_list=[10, 100, 1000])  # サブサンプリングなしの ZZ() に関して N 回実験
    ESSs_sum_Affine = zero(n_list)
    ESSs_sum_Global = zero(n_list)

    for _ in 1:N
        ESSs_Affine = []
        ESSs_Global = []
        for n in n_list
            push!(ESSs_Affine, getESSperEpoch(ab_Affine, T, dt, x[:,1:n], y[1:n]; ξ0=ξ0, θ0=θ0))
            push!(ESSs_Global, getESSperEpoch(ab_Global, T, dt, x[:,1:n], y[1:n]; ξ0=ξ0, θ0=θ0))
        end
        ESSs_sum_Affine = [ESSs_sum_Affine ESSs_Affine]
        ESSs_sum_Global = [ESSs_sum_Global ESSs_Global]
    end
    return mean(ESSs_sum_Affine, dims=2), var(ESSs_sum_Affine, dims=2), mean(ESSs_sum_Global, dims=2), var(ESSs_sum_Global, dims=2)
end

# ESS_Affine, var_ESS_Affine, ESS_Global, var_ESS_Global = experiment_ZZ(2, T, dt; ξ0=0.0, θ0=1.0, n_list=n_list)
実験には 10分 かかるので，保持した実行結果を用いる 
using JLD2

@load "Logistic2_Experiment1.jld2" ESS_Affine var_ESS_Affine ESS_Global var_ESS_Global
サブサンプリング付きの実験を実行する関数 experiment_ZZ() を定義 
using Statistics

# function ESS(samples::Vector{Float64}, T, dt)
#     B = T / dt
#     V = (dt / T) * sum(samples.^2) - ((dt / T) * sum(samples))^2
#     Y = samples .* sqrt(T / B)
#     ESS = T * V / var(Y)
#     return ESS
# end

function getESSperEpoch_SS(ab, ZZ, T ,dt, x, y; ξ0=0.0, θ0=1.0)
    trace, epochs, acc = ZZ(∇U, ξ0, θ0, T, x, y, ZigZag1d(); ab=ab)
    traj = discretize(trace, ZigZag1d(), dt)
    return ESS(traj.x, T, dt) * length(y) / epochs[end]  # サブサンプリングをしているので length(y) で補正する必要あり
end

N = 10
T = 500.0
dt = 0.1

function experiment_ZZ_SS(N, T, dt; ξ0=0.0, θ0=1.0, n_list=[10, 100, 1000])  # サブサンプリングなしの ZZ() に関して N 回実験
    ESSs_sum_CV = zero(n_list)
    ESSs_sum_SS = zero(n_list)

    for _ in 1:N
        ESSs_CV = []
        ESSs_SS = []
        for n in n_list
            push!(ESSs_CV, getESSperEpoch_SS(ab_Affine, ZZ1d_CV, T, dt, x[:,1:n], y[1:n]; ξ0=ξ0, θ0=θ0))
            push!(ESSs_SS, getESSperEpoch_SS(ab_Global, ZZ1d_SS, T, dt, x[:,1:n], y[1:n]; ξ0=ξ0, θ0=θ0))
        end
        ESSs_sum_CV = [ESSs_sum_CV ESSs_CV]
        ESSs_sum_SS = [ESSs_sum_SS ESSs_SS]
    end
    return mean(ESSs_sum_CV, dims=2), var(ESSs_sum_CV, dims=2), mean(ESSs_sum_SS, dims=2), var(ESSs_sum_SS, dims=2)
end

# ESS_CV, var_ESS_CV, ESS_SS, var_ESS_SS = experiment_ZZ_SS(2, T, dt; ξ0=0.0, θ0=1.0, n_list=n_list)
結果をプロット 
q = addPlot(q, n_list, ESS_CV, sqrt.(var_ESS_CV); label="ZZ-CV", color="darkorange")
q = addPlot(q, n_list, ESS_SS, sqrt.(var_ESS_SS); label="ZZ-SS", color="blue")
q = plot!(q, legend=:bottomleft)
display(q)

実際に用いたコードはこちら．

1.6 比較対象：MALA

MALA によるサンプリングを実行 
using AdvancedHMC, AdvancedMH, ForwardDiff
using LogDensityProblems
using LogDensityProblemsAD
using StructArrays
using LinearAlgebra

struct LogTargetDensity
    x::Matrix{Float64}
    y::Vector{Float64}
end

LogDensityProblems.logdensity(p::LogTargetDensity, ξ) = -U(ξ, p.x, p.y)
LogDensityProblems.dimension(p::LogTargetDensity) = 1
LogDensityProblems.capabilities(::Type{LogTargetDensity}) = LogDensityProblems.LogDensityOrder{0}()

model_with_ad = LogDensityProblemsAD.ADgradient(Val(:ForwardDiff), LogTargetDensity(x, y))

σ² = 0.0001
spl = MALA(x -> MvNormal((σ² / 2) .* x, σ² * I))

chain = sample(model_with_ad, spl, 2000; initial_params=ξ0, chain_type=StructArray, param_names=["ξ"], stats=true)

traj_MALA = chain.ξ

1.7 有効サンプル数について

時区間 \([0,T]\) における ZigZag 過程の，関数 \(h\in\mathcal{L}^2(\mathbb{R}^d)\) に関する 有効サンプル数 (ESS) とは \[ \widehat{\operatorname{ESS}}:=T\frac{\widehat{\mathrm{V}_\pi[h]}}{\widehat{\sigma^2_h}} \] \[ \widehat{\mathrm{V}_\pi[h]}:=\frac{1}{T}\int^T_0h(X_s)^2\,ds-\left(\frac{1}{T}\int^T_0h(X_s)\,ds\right)^2, \] \[ \widehat{\sigma^2_h}:=\frac{1}{B-1}\sum_{i=1}^B(Y_i-\overline{Y})^2,\quad Y_i:=\sqrt{\frac{B}{T}}\int^{\frac{iT}{B}}_{\frac{(i-1)T}{B}}h(X_s)\,ds \] で定まる値である．

例えば次のようにして計算できる：

function ESS(samples::Vector{Float64}, T, dt)
    V = (dt / T) * sum(samples.^2) - ((dt / T) * sum(samples))^2
    Y = samples .* sqrt(T / B)
    ESS = T * V / var(Y)
    return ESS
end

1.8 重点サブサンプリング (Sen et al., 2020)

一様でないサブサンプリングを導入することで，Zig-Zag サンプラーを不均衡データにも強くすると同時に，サブサンプリングの効率を上げることもできる．

強度関数 \[ m_i^K(t)=\biggr(\theta_iE^K_i(x+\theta t)\biggl)_+ \] は， \[ \operatorname{E}\biggl[E^K_i(\xi)\biggr]=\partial_iU(\xi) \] を満たす限り，\(K\sim\mathrm{U}([n])\) に限る必要はなかったのである．

すなわち，\((p_x)\) をある \([n]\) 上の分布 \(\nu\in\mathcal{P}([n])\) の質量関数として \[ \partial_iV_1^J(\xi):=\frac{1}{p_J}\partial_iU^J(\xi)\qquad J\sim\nu \] と定めると， \[ \lvert\partial_iV_i^j(\xi)\rvert\le\max_{j\in[n]}\frac{\lvert x_i^j\rvert}{p_j} \] が成り立つ．式 (1) は \(p_j\equiv1/n\) の場合であったのである．

換言すれば， \[ p_j\,\propto\,\lvert x^j_i\rvert \] と定めることで，Poisson 強度関数 \(m^j_i\) の上界をタイトにすることができ，その結果剪定の効率が上がる．

a_IS(ξ, θ, x, y) = sum(abs.(x))
b_IS(ξ, θ, x, y) = 0

ab_IS(ξ, θ, x, y, ::ZigZag1d) = (a_IS(ξ, θ, x, y), b_IS(ξ, θ, x, y))
重点サブサンプリングによる ZigZag サンプラー ZZ1d_IS() を定義 
using StatsBase

function λj_IS(j::Int64, ξ, θ, ∇U, x, y)
    pj = abs(x[1,j]) / sum(abs.(x))  # x がスパースだと 0 になりやすいことに注意
    Eʲ = ∇U(1, j, ξ, x, y) / (length(y) * pj)
    return pos(θ * Eʲ)
end

function ZZ1d_IS(∇U, ξ, θ, T::Float64, x::Matrix{Float64}, y::Vector{Float64}, Flow::ZigZagBoomerang.ContinuousDynamics; rng=Random.GLOBAL_RNG, ab=ab_IS)
    t = zero(T)
    Ξ = [(t, ξ, θ)]
    num = acc = 0
    epoch_list = [num]
    a, b = ab(ξ, θ, x, y, Flow)
    t′ =  t + poisson_time(a, b, rand())  # イベントは a,b が定める affine proxy に従って生成する
    n = length(y)

    while t < T
        τ = t′ - t
        t, ξ, θ = move_forward(τ, t, ξ, θ, Flow)
        j = sample(1:n, Weights(abs.(vec(x))))
        l, lb = λj_IS(j, ξ, θ, ∇U, x, y), λ_bar(τ, a, b)  # λ が真のレート, λ_bar が affine proxy
        num += 1
        if rand()*lb < l
            acc += 1
            if l > lb + 0.01
                println(l-lb)
            end
            θ = -θ
            push!(Ξ, (t, ξ, θ))
            push!(epoch_list, num)
        end
        a, b = ab(ξ, θ, x, y, Flow)
        t′ = t + poisson_time(a, b, rand())
    end

    return Ξ, epoch_list, acc/num
end
実験を実行する関数 experiment_ZZ_IS() を定義 
function experiment_ZZ_IS(N, T, dt; ξ0=0.0, θ0=1.0, n_list=[10, 100, 1000])  # 重点サブサンプリング ZZ1d_IS() に関して N 回実験
    ESSs_sum_IS = zero(n_list)

    for _ in 1:N
        ESSs_IS = []
        for n in n_list
            push!(ESSs_IS, getESSperEpoch_SS(ab_IS, ZZ1d_IS, T, dt, x[:,1:n], y[1:n]; ξ0=ξ0, θ0=θ0))
        end
        ESSs_sum_IS = [ESSs_sum_IS ESSs_IS]
    end
    return mean(ESSs_sum_IS, dims=2), var(ESSs_sum_IS, dims=2)
end

ESS_IS, var_ESS_IS = experiment_ZZ_IS(2, 500.0, 1.0; ξ0=0.0, θ0=1.0, n_list=n_list)
結果をプロット 
r = addPlot(q, n_list, ESS_IS, sqrt.(var_ESS_IS); label="ZZ-IS", color="green")
display(r)

制御変数による方法は \(n\to\infty\) の漸近論に基づいているので，観測数が増えるほど効率は上がっていく．

1.9 大規模不均衡データ

大規模不均衡データでは，事後分布が十分な集中性を持たないために制御変数による方法 ZZ-CV が十分な効率改善を示さないが，重点サブサンプリングによれば Poisson 強度関数のタイトな上界を引き続き構成できる．

ここでは，\(\xi_0=1\) を真値とし，次のような１次元データを考える： \[ X^j\overset{\text{iid}}{\sim}(1-\alpha)\delta_0+\alpha\mathrm{N}(1,2) \] \[ \operatorname{P}[Y^j=1]=\frac{1}{1+e^{-X^j}} \]

ξ0 = [1.0] # True value
Σ = [2]
n = 1000

function sample_SparseData(n::Int64, α::Float64; ρ=MvNormal(ξ0, Σ))
    x = []
    while length(x) < n
        rand() < α ? push!(x, rand(ρ)[1]) : push!(x, 0.0)
    end
    x = Float64.(reshape(x,1,:))
    y = rand.(Bernoulli.(logistic.(ξ0*x)))
    y = Float64.(vec(y))
    return x, y
end

α_list = [1, 0.1, 0.01]

x_Sparse, y_Sparse = [], []

for α in α_list
    x, y = sample_SparseData(n, α)
    push!(x_Sparse, x)
    push!(y_Sparse, y)
end
using Optim

a_Sparse(ξ_star, ξ, θ, x, y) = pos(θ * ∇U(ξ_star,x,y)) + C(ξ, θ, x, y) * abs(ξ - ξ_star[1])
ξ_star_list = []
α_list = [1, 0.1, 0.01]

for α in 1:length(α_list)
    result = optimize(ξ -> U(ξ, x_Sparse[α], y_Sparse[α]), [0.0], LBFGS())
    ξ = Optim.minimizer(result)
    push!(ξ_star_list, ξ)
end

function experiment_Sparse(N, T, dt; ξ0=0.0, θ0=1.0, α_list=[1, 0.1, 0.01, 0.001])
    # ESSs_sum_CV = zero(α_list)
    ESSs_sum_SS = zero(α_list)
    ESSs_sum_IS = zero(α_list)

    for _ in 1:N
        # ESSs_CV = []
        ESSs_SS = []
        ESSs_IS = []
        for i in 1:length(α_list)
            # ab_Sparse(ξ, θ, x, y, ::ZigZag1d) = (a_Sparse(ξ_star_list[i], ξ, θ, x, y), b_Affine(ξ, θ, x, y))
            # push!(ESSs_CV, getESSperEpoch_SS(ab_Sparse, ZZ1d_CV, T, dt, x_Sparse[i], y_Sparse[i]; ξ0=ξ0, θ0=θ0))
            push!(ESSs_SS, getESSperEpoch_SS(ab_Global, ZZ1d_SS, T, dt, x_Sparse[i], y_Sparse[i]; ξ0=ξ0, θ0=θ0))
            push!(ESSs_IS, getESSperEpoch_SS(ab_IS, ZZ1d_IS, T, dt, x_Sparse[i], y_Sparse[i]; ξ0=ξ0, θ0=θ0))
        end
        # ESSs_sum_CV = [ESSs_sum_CV ESSs_CV]
        ESSs_sum_SS = [ESSs_sum_SS ESSs_SS]
        ESSs_sum_IS = [ESSs_sum_IS ESSs_IS]
    end
    # return mean(ESSs_sum_CV, dims=2), var(ESSs_sum_CV, dims=2), mean(ESSs_sum_SS, dims=2), var(ESSs_sum_SS, dims=2), mean(ESSs_sum_IS, dims=2), var(ESSs_sum_IS, dims=2)
    return mean(ESSs_sum_SS, dims=2), var(ESSs_sum_SS, dims=2), mean(ESSs_sum_IS, dims=2), var(ESSs_sum_IS, dims=2)
end

N = 2
T = 500.0
dt = 0.1

ESS_SS, var_ESS_SS, ESS_IS, var_ESS_IS = experiment_Sparse(N, T, dt; ξ0=0.0, θ0=1.0, α_list=α_list)
プロット 
using LaTeXStrings
p = startPlot(α_list, ESS_SS, sqrt.(var_ESS_SS); label="ZZ-SS", xlabel=L"Sparsity $\alpha$", color="blue"#, background_color=true
)
p = addPlot(p, α_list, ESS_IS, sqrt.(var_ESS_IS); label="ZZ-IS", color="green")
p = addPlot(p, α_list, ESS_CV, sqrt.(var_ESS_CV); label="ZZ-CV", color="darkorange")
display(p)

ここでは問題にしないが，圧倒的に実行時間が重点サブサンプリングの方が短い．

\(\alpha<10^{-3}\) の領域では動作が不安定になる．論文 (Sen et al., 2020) でもこの領域は触れられていない．しかし， \[ \#\left\{i\in[n]\mid y^i=1\right\}\approx500 \] であるため，特に理由は見つからない．

2 高次元へのスケーリング

さらに，(Bierkens et al., 2022) による方法で，スパースデータに対する効率化が可能である．

using StatsFuns
using Distributions

"""
    U(ξ, x, y)：ポテンシャル関数
        ξ: パラメータ空間上の点
        (x,y): 観測
"""
function U(ξ, x, y)
    n = length(y)
    U_list = []
    for j in 1:n
        push!(U_list, U(j, ξ, x, y))
    end
    return mean(U_list)
end
function U(j, ξ, x, y)
    n = length(y)
    product = dot(x[:,j],ξ)
    return -n * exp(y[j] * product) / (1 + exp(product))
end

ξ0 = [1,2]  # True value
n_list = [10, 100, 1000, 10000]  # 実験で用いるサンプルサイズの列

Σ = [2 0; 0 2]
x = rand(MvNormal(ξ0, Σ), n_list[end])
y = rand.(Bernoulli.(logistic.(ξ0'*x)))  # BitVector になってしまう
y = Float64.(vec(y))  #  Vector{Float64} に変換
using Optim

result = optimize(ξ -> U(ξ,x,y), [0.0,0.0], LBFGS())
ξ_star = Optim.minimizer(result)

References

Bierkens, J., Grazzi, S., Meulen, F. van der, and Schauer, M. (2022). Sticky PDMP Samplers for Sparse and Local Inference Problems. Statistics and Computing, 33(1), 8.
Sen, D., Sachs, M., Lu, J., and Dunson, D. B. (2020). Efficient posterior sampling for high-dimensional imbalanced logistic regression. Biometrika, 107(4), 1005–1012.