Alpha安定分布（Alpha-Stable Distribution）は、確率密度関数（PDF）が解析的に表現できない分布である。ただし、ガウス分布やコーシー分布などの特殊ケースを除き、数値的手法や近似アルゴリズムによって計算可能である。

分布の基本パラメータ

Alpha安定分布は以下の4つのパラメータで特徴づけられる：

• 安定性指数（α）：分布の裾の厚さを調整する（0 < α ≤ 2）。α = 2 のとき、ガウス分布に一致する。
• 歪度係数（β）：分布の非対称性を表す（-1 ≤ β ≤ 1）。
• 尺度パラメータ（γ）：分散のような役割を持つが、厳密には異なる（γ > 0）。
• 位置パラメータ（μ）：分布の中心を示す。

この分布の確率密度関数は特性関数を通じて定義される。

MATLABによる実装手法

1. 外部ツールボックスの活用

Stable Distribution Toolbox を使用することで、PDFの評価が容易になる。

% Alpha=1.5, β=0, γ=1, μ=0におけるPDFの例
alpha = 1.5;
beta = 0;
gamma = 1;
mu = 0;
x = linspace(-5, 5, 1000);
pdf_vals = stable_pdf(x, alpha, beta, gamma, mu);  % ツールボックスが必要
plot(x, pdf_vals);

2. 自作数値積分関数

任意のパラメータに対応したPDFを求めるために、数値積分を用いた独自関数を構築する。

function pdf_result = stable_pdf_custom(x_point, alpha, beta, gamma, mu)
    if nargin < 5, mu = 0; end
    if nargin < 4, gamma = 1; end
    if nargin < 3, beta = 0; end
    if nargin < 2, alpha = 1.5; end

    % 積分区間設定
    t_vals = linspace(eps, 10, 1000);
    dt = t_vals(2) - t_vals(1);

    % 特性関数の被積分関数
    integrand_func = @(t) sin((x_point - mu) * t / gamma) ./ t .* ...
        exp(-abs(t).^alpha .* (1 - 1i * beta * sign(t) * tan(pi * alpha / 2)));

    % 台形則による積分
    integral_value = trapz(t_vals, integrand_func(t_vals));
    pdf_result = (1 / (pi * gamma)) * integral_value;
end

3. 近似アルゴリズムの適用

Nolanによる近似法を元にした簡易的なPDF計算関数。

function approx_pdf = fast_stable_pdf(x_input, alpha, beta, gamma, mu)
    assert(alpha > 0 && alpha <= 2, 'Alpha must be within (0, 2]');
    assert(abs(beta) <= 1, 'Beta must be within [-1, 1]');

    beta_mod = beta * tan(pi * alpha / 2) / (alpha ~= 1);
    t_domain = linspace(eps, 10, 1000);

    integrand_part = zeros(size(t_domain));
    for idx = 1:length(t_domain)
        t_current = t_domain(idx);
        integrand_part(idx) = exp(-1i * t_current * (x_input - mu) + ...
            abs(t_current)^alpha * (1 - 1i * beta_mod * sign(t_current)));
    end

    integrated_sum = trapz(t_domain, integrand_part);
    approx_pdf = (1 / (pi * gamma)) * real(integrated_sum);
end

全体コード例

% パラメータ指定
alpha_val = 1.3;
beta_val = 0.4;
gamma_val = 1.0;
mu_val = 0;
range_x = linspace(-5, 5, 1000);

% 各手法によるPDF計算
pdf_numerical = arrayfun(@(x) stable_pdf_custom(x, alpha_val, beta_val, gamma_val, mu_val), range_x);
pdf_fast = arrayfun(@(x) fast_stable_pdf(x, alpha_val, beta_val, gamma_val, mu_val), range_x);

% グラフ表示
figure;
plot(range_x, pdf_numerical, 'b-', 'LineWidth', 1.5);
hold on;
plot(range_x, pdf_fast, 'r--', 'LineWidth', 1.5);
xlabel('x'); ylabel('PDF');
legend('Numerical Integration', 'Fast Approximation');
title(sprintf('Alpha-Stable PDF (α=%.1f, β=%.1f)', alpha_val, beta_val));
grid on;

性能向上テクニック

• 積分範囲の調整：αが小さい場合、より広い範囲での積分が必要。
• ベクトル化処理：arrayfun や parfor を使用して高速化。
• 数値安定性の確保：複素数の指数計算は対数領域で行うことで精度を向上。

応用例

金融データへの適合

株式収益率のような厚尾分布のモデリングに適している。

% 株価データ読み込み
data_table = readtable('stock_returns.csv');
returns_data = data_table.Returns;

% パラメータ推定（最尤推定などが必要）
[est_alpha, est_beta] = fit_alpha_stable_params(returns_data);

% 分布フィッティング結果のプロット
histogram(returns_data, 'Normalization', 'pdf');
hold on;
fit_range = linspace(min(returns_data), max(returns_data), 1000);
plot(fit_range, stable_pdf_custom(fit_range, est_alpha, est_beta, 1, mean(returns_data)), 'r');

信号処理におけるノイズモデル

レーダーのクラッターなど、重尾ノイズを考慮したシミュレーション。

% 信号生成
time_axis = 0:0.01:1;
signal_clean = sin(2 * pi * 5 * time_axis);
noise_alpha = generate_alpha_noise(1.5, 0, 0.5, 0, length(signal_clean));
signal_noisy = signal_clean + noise_alpha';

% ノイズ除去（ウェーブレット変換）
signal_denoised = wdenoise(signal_noisy, 4);