安定分布からのサンプリングとフィッティングの実践

参考文献, リンクは末尾に述べます.

安定分布という広い分布のクラスがあるのですが, その分布からのサンプリングやフィッティング方法について実装込みでの解説(特に日本語のもの)は多くないように思いますので, 筆者の勉強を兼ねてここに記します. 良い資料などございましたらご教示いただけますと誠に幸いです.

安定分布の定義の方法にはいくつかの種類が知られているため(例えば2023.09.27現在で日本語版と英語版のWikipediaで特性関数の式が微妙に異なります), それらが等価であるかの確認を行うとともにPDF, サンプリング方法, パラメータ推定方法が全て一貫した定義で実装されているか注意が必要と思われます. 本稿では一貫して次の定義を用います. PDF$f(x; \alpha, \beta, \gamma, \delta)$, 特性関数(以下CF)$\phi(t; \alpha, \beta, \gamma, \delta)$は次で与えられるとします.

\begin{align}
f(x; \alpha, \beta, \gamma, \delta) &= \frac{1}{2\pi} \int_{-\infty}^{\infty} \phi(t; \alpha, \beta, \gamma, \delta) e^{-i t x} dt \\
\phi(t; \alpha, \beta, \gamma, \delta) &= \exp \bigg[ it\delta - \gamma |t|^\alpha \big(1 - i \beta \textrm{sign}(t) \omega(t; \alpha, \gamma) \big) \bigg] \\
\omega(t; \alpha, \gamma) &= \begin{cases} 
\tan \frac{\pi \alpha}{2} & \alpha \neq 1\\
-\frac{2}{\pi} \log (\gamma |t|) & \alpha = 1
\end{cases}
\end{align}

他の文献にあたる際は尺度パラメータ$\gamma$の定義, 歪度パラメータ$\beta$の前の符号, $\omega$の定義などに注意します.

安定分布はさまざまな分布を特殊な場合に含みます. 人の名前のついている例として, $(\alpha, \beta) = (2, 0)$はガウス分布, $(\alpha, \beta) = (0.5, 1)$はレヴィ分布, $(\alpha, \beta) = (1, 0)$はコーシー分布に対応しています.

CFの定義式から次が成立します.

\begin{align}
\log (-\log |\phi(t)|) &= \alpha \log |t| + \log \gamma\\
\frac{1}{t} \textrm{arctan} \frac{\Im \phi(t)}{\Re \phi(t)} &= \beta \big( \gamma |t|^{\alpha-1} \omega(t; \alpha, \gamma) \big) + \delta\\
\omega(t; \alpha, \gamma) &= \begin{cases} 
\tan \frac{\pi \alpha}{2} & \alpha \neq 1\\
-\frac{2}{\pi} \log (\gamma|t|) & \alpha = 1
\end{cases}
\end{align}

上の式はCFの絶対値や偏角と$\alpha, \beta, \gamma, \delta$の意味を結びつける式と見ることができ, CFの絶対値は分布の中心付近や裾での振る舞いに関係し, 偏角は分布の歪みに関係していると大雑把に理解できます. CFは一般に有界$|\phi(t)| \leq 1$かつエルミート$\phi(-t) = \overline{\phi(t)}$であるので, 上2式は意味を持ち, かつ$t \geq 0$の領域のみ考えて一般性を失いません. 絶対値は$t$の増加に対して単調減少で, $t=0$の時は$(\Re \phi(0), \Im \phi(0)) = (1, 0)$であり, また$t \rightarrow \infty$で$(\Re \phi(t), \Im \phi(t)) \rightarrow (0, 0)$となることがリーマンルベーグの補題からわかります. 参考までにLevy分布$(\alpha, \beta, \gamma, \delta)=(0.5, 1, 1, 0)$の場合でのCFを次の図に示します. 左と真ん中は$\phi(t)$の実部, 虚部を横軸を波数$t$として描画したもの, 右は複素数平面内の半径1の円板に$\phi(t)$を描画したものです. 波数を等間隔0.1で切って計算しています. ⚪︎は以下に示すサンプリング方法で生成したデータから計算した経験的特性関数(ECF), 点線はCFの解析式です.

サンプリング方法実装

完成図は以下です. パラメータを指定した際の, ⚪︎印がサンプルデータから作成したPDF, 実線がCFからフーリエ変換して得られたPDFです. 左では$\beta = 0, \gamma = 1, \delta = 0$としたときの様々な$\alpha$に対応するPDF, 中では$\alpha = 0.5, \gamma = 1, \delta = 0$としたときの様々な$\beta$に対応するPDF, 右では$\alpha = 1, \gamma = 1, \delta = 0$とした時の様々な$\beta$に対応するPDFを示しています.

作成時の条件としては, データポイント数は100000, PDFのbin幅は0.1, CFでカットオフ波数$\pm50$までを使っています. 数秒で上3つの図を作れます.

以下のコードのうち, def generate以下でサイズ$n$, パラメータ$\alpha, \beta, \gamma, \delta$の安定分布からのサンプルを生成します. Weronの方法(参考文献[2], [3], URL[4])を実装しています.

生成コード

import numpy as np
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy.fft import fft, fftfreq, ifft, fftshift

def pdfn(ndarray, bins, ax=None, scale='linear', kwargs={}):
  cnts = pd.cut(ndarray, bins).value_counts()
  p = cnts / ndarray.shape[0]
  vals = p / np.diff(bins)
  if scale=='linear':
    mb = np.array([ (bins[i]+bins[i+1])/2. for i in np.arange(bins.shape[0]-1) ])
  elif scale=='log':
    mb = np.array([ np.sqrt(bins[i]*bins[i+1]) for i in np.arange(bins.shape[0]-1) ])
  else:
    raise ValueError("scale should be either 'linear' or 'log'")
  if ax!=None:
    ax.plot(mb, vals, **kwargs)
  return mb, vals, cnts, p

def char(t: float, param: np.ndarray):
  alpha, beta, gamma, delta = param
  if alpha == 1.0:
    omega = - 2.0 / np.pi * np.log(gamma * np.abs(t))
  else:
    omega = np.tan(np.pi * alpha / 2.0)
  return np.exp( 1j * delta * t - gamma * (np.abs(t))**(alpha) * (1 - 1j * beta * np.sign(t) * omega) )


def pdf(param: np.ndarray):
  dt = 0.01
  t = np.arange(-50, 50, dt)
  cf = char(t, param)
  f = 1.0 / (2.0 * np.pi) * dt * np.abs(fftshift(fft(cf)))
  x = 2.0 * np.pi / dt * fftshift(fftfreq(t.shape[0]))
  return x, f

def generate(n: int, param: np.ndarray):
  alpha, beta, gamma, delta = param
  V = np.random.default_rng(seed=42).uniform(low=-np.pi/2.0, high=np.pi/2.0, size=n)
  W = np.random.default_rng(seed=42).exponential(scale = 1.0, size=n)

  if alpha == 1.0:
    X = 2.0 / np.pi * ( (np.pi / 2.0 + beta * V) * np.tan(V) - beta * np.log( (np.pi/2.0 * W * np.cos(V)) / (np.pi/2.0 + beta * V) ) )
  else:
    B = np.arctan(beta * np.tan(np.pi * alpha / 2.0)) / alpha
    S = (1 + beta**2 * (np.tan(np.pi * alpha / 2.0)**2)) ** ( 1.0 / (2.0*alpha) )
    X = S * np.sin(alpha * (V + B)) / (np.cos(V) ** (1.0 / alpha)) * (np.cos(V - alpha * (V + B)) / W) ** ( (1.0-alpha) / alpha )

  return gamma ** (1.0/alpha) * X + delta

def ft(x: np.ndarray):
  dt = 0.1
  t = np.arange(0, 5, dt)
  def dft(t: float):
    return np.exp(1j * t * x).sum()
  np_dft = np.frompyfunc(dft, 1, 1)
  f = (1.0 / x.shape[0] * np_dft(t)).astype('<c8')
  return t, f

def plot_samples():

  color = ['k', 'b', 'g', 'r', 'm']
  marker = ['o', 'x', '^', 'v', 'd']

  alpha = [0.5, 1.0, 1.5, 2.0]
  beta  = 0.5
  gamma = 1.0
  delta = 0.0

  bins = np.linspace(-5, 5)
  fig, ax = plt.subplots()
  ax.set_xlim([-5, 5])
  ax.set_ylim([0, 0.7])
  for i in range(4):
    param = [alpha[i], beta, gamma, delta]
    x = generate(n = 100000, param = param)
    mb, val, cnt, p = pdfn(x, bins = bins, scale='linear')
    ax.scatter(mb, val, fc='none', ec=color[i], marker=marker[i], label=fr'$\alpha$={alpha[i]:.1f}')
    print(pdf(param))
    x, f = pdf(param)
    ax.plot(x, f, c=color[i])
  ax.legend(fontsize=12)
  fig.savefig(f'generate-stable-alpha.png')

  alpha = 2.0
  beta  = [0.0, 0.25, 0.5, 0.75, 1.0]
  gamma = 1.0
  delta = 0.0

  bins = np.linspace(-5, 5)
  fig, ax = plt.subplots()
  ax.set_xlim([-5, 5])
  ax.set_ylim([0, 0.7])
  for i in range(5):
    param = [alpha, beta[i], gamma, delta]
    x = generate(n = 100000, param = param)
    mb, val, cnt, p = pdfn(x, bins = bins, scale='linear')
    ax.scatter(mb, val, fc='none', ec=color[i], marker=marker[i], label=fr'$\beta$={beta[i]:.2f}')
    x, f = pdf(param)
    ax.plot(x, f, c=color[i])
  ax.legend(fontsize=12)
  fig.savefig(f'generate-stable-beta-1.png')

  alpha = 1.0
  beta  = [0.0, 0.25, 0.5, 0.75, 1.0]
  gamma = 1.0
  delta = 0.0

  bins = np.linspace(-5, 5)
  fig, ax = plt.subplots()
  ax.set_xlim([-5, 5])
  ax.set_ylim([0, 0.7])
  for i in range(5):
    param = [alpha, beta[i], gamma, delta]
    x = generate(n = 100000, param = param)
    mb, val, cnt, p = pdfn(x, bins = bins, scale='linear')
    ax.scatter(mb, val, fc='none', ec=color[i], marker=marker[i], label=fr'$\beta$={beta[i]:.2f}')
    x, f = pdf(param)
    ax.plot(x, f, c=color[i])
  ax.legend(fontsize=12)
  fig.savefig(f'generate-stable-beta-1.png')

if __name__ == '__main__':
  plot_samples()

フィッティング方法実装

フィッティングの際に気になることがあったので方法の詳細とフィッティングの過程を示します.

Koutrouvelisの方法の第1近似

データ$x_j, j=1, 2, \cdots, N$から経験的特性関数(ECF)$\hat{\phi}_{N}(t)$を次式で生成します.

\hat{\phi}_{N}(t) = \frac{1}{N} \sum_{j=1}^{N} e^{i t x_j}

データが安定分布から得られているのであれば, ECFを本稿冒頭に示したCFの絶対値と偏角が満たす2本の式を上から順にフィットすることにより係数を決めていくことができます.

Koutrouvelisの方法はこれらの式のフィッティングとデータの標準化, さらに波数の分点の取り直しを組み合わせて再帰的にパラメータ推定を繰り返していくものです.

例

以下の2例では, 正解値を用いてWeronの方法でサンプルを生成し, 上に説明した第一近似の方法でそれに対するフィッティングを行って推定値と正解が一致するか見てみます.

例 レヴィ分布: $\alpha = 0.5, \beta = 1.0, \gamma = 1.0, \delta = 0.0$の場合

ECFの絶対値, 偏角に対応するフィッティングの図を示します. 上の1本目の式に対応する散布図とフィッティングの様子を左図に示します. 横軸が$\log |t|$, 縦軸が左辺です. 青⚪︎がデータ, 赤線がフィッティング結果です. また, 左図での赤線に対応する$\alpha, \gamma$を用いて2本目の式に対応する散布図とフィッティングの様子を右図に示します. 横軸が$\gamma^{\alpha} \omega(t; \alpha) |t|^{\alpha-1}$, 縦軸が左辺です. 青⚪︎がデータ, 赤線が$\arg \phi(t)>0$の部分のみ用いたフィッティング結果です. 解析的な式と非常に近い振る舞いをしています.

正解は$(\alpha, \beta, \gamma, \delta) = (0.5, 1, 1, 0)$, 推定結果は$(\alpha, \beta, \gamma, \delta) = (0.499(1), 1.007(2), 0.999(1), -0.001(1))$でした.

例 ガウス分布： $\alpha = 2, \beta = 0, \gamma = 1, \delta = 0$の場合

上の例で示した方法と同様の流れで今回のパラメータで結果を示します. 青⚪︎はデータから求めたECFを表します. 左図について, 黒線が全データを利用したフィッティング結果, 赤線が縦軸最大値をとる点より左だけでフィッティングした結果です. 右図については見にくいですが黒線が(0, 0.05)付近にあり全データを用いたフィッティング結果を表し, 赤線が$\arg \phi(t)>0$の部分のみ用いたフィッティング結果です. 複素数の偏角についてはデータが複素数平面上でブランチを跨いだかどうか考慮する必要があり, 特定のブランチのみに属するデータだけでフィットする工夫を行いました.

上の左図において, 左側の領域はCFの波数0付近に対応しており, PDFでいう裾にあたります. 右側のデータまで用いた黒線のフィッティングがずれてしまっています. しかしながら黒線の$R^2$は0.924であり, さもよくフィットできているかのように思えてしまうので注意が必要です. 裾だけフィットした赤線に対応する$R^2$は1.000でした. また, 上の右図において, $\alpha \approx 2$において横軸の量がほぼ0になることと, $\arg \phi(t)$が完全にゼロにならなかったことが影響してフィットが正しくできなくなっています. 安定分布において$\alpha = 2$の時は$\beta$の値は何であっても分布は同じになるため理論上は問題ないのですが, $\delta$の推定に影響が出る可能性があります. 正解が$\beta = \delta = 0$であるのもあり, この状況でパラメータ推定する場合には上右図の代替手段で$\delta$の推定値を求める方法を持っておきたいです.

また, 考察を進めるために今回のパラメータの正解値におけるCFを上に示しました. ECFの絶対値のフィッティングの際の$\log |t|$が大きい部分は複素円板におけるCFの原点付近にあたり, その情報はPDFの滑らかさに関係しているので, 波数が大きい時のECFの原点への向かい方がサンプルサイズの有限性やECF計算時の波数の分点の取り方などによって影響を受けやすい領域と見ることができます. ECFの絶対値のフィッティングの図において縦軸の最大値は最も原点に近くなったECFの値に対応しているので, 原点付近の影響を除くのに縦軸最大値を利用するのは理に適っているように思われます. 適切な閾値を定めてPDFの裾側のみを$\alpha$のフィッティングに用いるのが有効に見えます.

上の赤線のフィッティングで, 正解$(\alpha, \beta, \gamma, \delta) = (2, 0, 1, 0)$, 推定結果$(\alpha, \beta, \gamma, \delta) = (2.014(4), 1.59(33), 1.014(3), -0.005(18))$でした.

上の知見を踏まえて, ECFを求めたのちにその絶対値と偏角を求め,$t$増加とともに絶対値が単調に減少している部分の絶対値と偏角だけ抜き出してフィッティングに使用した場合, 次の図を得ました.

推定結果は$(\alpha, \beta, \gamma, \delta) = (2.014(4), 2.177(308), 1.014(3), -0.030(11))$でした.

$\alpha \approx 2$の時や$\beta \approx 0$の時に$|t| \gg 1$領域で波数の分点の取り方やサンプルの影響を受けてECFの実部や虚部が正になったり負になったりすると推定精度に悪影響が及びやすく, そのような場合の対処方法を用意しておくことが重要となりそうです.

この第一近似は理解しやすく実装も容易ですが, データを変えるたびに個別対応が必要そうなので, 実データに適用することを見据えた自動化がこのままでは容易でないように思えます.

生成コード

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy.fft import fft, fftfreq, ifft, fftshift

def char(t: float, param: np.ndarray):
  alpha, beta, gamma, delta = param
  if alpha == 1.0:
    omega = - 2.0 / np.pi * np.log(gamma * np.abs(t))
  else:
    omega = np.tan(np.pi * alpha / 2.0)
  return np.exp( 1j * delta * t - gamma * (np.abs(t))**(alpha) * (1 - 1j * beta * np.sign(t) * omega) )

def pdf(param: np.ndarray):
  dt = 0.01
  t = np.arange(-50, 50, dt)
  cf = char(t, param)
  f = 1.0 / (2.0 * np.pi) * dt * np.abs(fftshift(fft(cf)))
  x = 2.0 * np.pi / dt * fftshift(fftfreq(t.shape[0]))
  return x, f

def generate(n: int, param: np.ndarray):
  alpha, beta, gamma, delta = param
  V = np.random.default_rng(seed=42).uniform(low=-np.pi/2.0, high=np.pi/2.0, size=n)
  W = np.random.default_rng(seed=42).exponential(scale = 1.0, size=n)

  if alpha == 1.0:
    X = 2.0 / np.pi * ( (np.pi / 2.0 + beta * V) * np.tan(V) - beta * np.log( (np.pi/2.0 * W * np.cos(V)) / (np.pi/2.0 + beta * V) ) )
  else:
    B = np.arctan(beta * np.tan(np.pi * alpha / 2.0)) / alpha
    S = (1 + beta**2 * (np.tan(np.pi * alpha / 2.0)**2)) ** ( 1.0 / (2.0*alpha) )
    X = S * np.sin(alpha * (V + B)) / (np.cos(V) ** (1.0 / alpha)) * (np.cos(V - alpha * (V + B)) / W) ** ( (1.0-alpha) / alpha )

  return gamma ** (1.0/alpha) * X + delta

def ft(x: np.ndarray):
  dt = 0.1
  t = np.arange(0, 5, dt)
  def dft(t: float):
    return np.exp(1j * t * x).sum()
  np_dft = np.frompyfunc(dft, 1, 1)
  f = (1.0 / x.shape[0] * np_dft(t)).astype('<c8')
  return t, f


def fit(x: np.ndarray, ax: None):
  t, f = ft(x)
  preprocessdata = np.vstack((t, f)).T
  preprocessdata = preprocessdata[np.log(np.abs(preprocessdata[:, 0])) > -5]
  t, f = preprocessdata[:, 0], preprocessdata[:, 1]
  fitdata = np.vstack((np.log(np.abs(t)), np.log(-np.log(np.abs(f))))).T
  fitdata = np.nan_to_num(fitdata, nan = np.nan, posinf = np.nan, neginf = np.nan)
  fitdata = fitdata[~np.isnan(fitdata).any(axis=1)]
  X = sm.add_constant(fitdata[:, 0])
  y = fitdata[:, 1]
  res = sm.OLS(y, X).fit()
  print(' --- alpha and gamma fit using all data --- ')
  print(res.summary())
  alpha = res.params[1]
  gamma = np.exp(res.params[0])

  ax[0].set_xlabel(r'$\log |t|$')
  ax[0].set_ylabel(r'$\log (-\log |\phi(t)|)$')
  ax[0].scatter(fitdata[:, 0], fitdata[:, 1], fc='none', ec='b')
  ax[0].plot(X[:, 1], res.predict(X), c='k', lw=0.5)

  xmax = np.argmax(fitdata[:, 1])
  # ax[0].axvline(fitdata[xmax, 0], c='r', ls='--')
  fitdata = fitdata[fitdata[:, 0] < fitdata[xmax, 0]]
  X = sm.add_constant(fitdata[:, 0])
  y = fitdata[:, 1]
  res = sm.OLS(y, X).fit()
  print(' --- alpha and gamma fit using tail info --- ')
  print(res.summary())
  alpha = res.params[1]
  gamma = np.exp(res.params[0])
  ax[0].plot(X[:, 1], res.predict(X), c='r', ls='--')

  if alpha == 1.0:
    omega = - 2.0/np.pi * np.log(np.abs(t))
  else:
    omega = np.tan(np.pi * alpha / 2.0)

  # fitdata = np.vstack((gamma * omega * np.abs(t)**(alpha - 1), 1.0 / t * np.arctan(f.imag / f.real))).T.astype('<f8')
  fitdata = np.vstack((gamma * omega * np.abs(t)**(alpha - 1), 1.0 / t * np.angle(f))).T.astype('<f8')
  fitdata = np.nan_to_num(fitdata, nan = np.nan, posinf = np.nan, neginf = np.nan)
  fitdata = fitdata[~np.isnan(fitdata).any(axis=1)]
  X = sm.add_constant(fitdata[:, 0])
  y = fitdata[:, 1]
  res = sm.OLS(y, X).fit()
  print(' --- beta and delta fit using all data --- ')
  print(res.summary())
  beta = res.params[1]
  delta = res.params[0]

  ax[1].set_xlabel(r'$\gamma t^{\alpha-1} \omega(t; \alpha)$')
  ax[1].set_ylabel(r'$\frac{1}{t}$ arg $\phi(t)$')
  ax[1].set_xlim([-2, 2])
  ax[1].scatter(fitdata[:, 0], fitdata[:, 1], fc='none', ec='b')
  ax[1].plot(X[:, 1], res.predict(X), c='k')

  fitdata = fitdata[fitdata[:, 1]>0]
  X = sm.add_constant(fitdata[:, 0])
  y = fitdata[:, 1]
  res = sm.OLS(y, X).fit()
  print(' --- beta and delta fit using tail info --- ')
  print(res.summary())
  beta = res.params[1]
  delta = res.params[0]
  ax[1].plot(X[:, 1], res.predict(X), c='r', ls='--')

  return alpha, beta, gamma, delta


def plot_cf(t: np.ndarray, cf: np.ndarray, x: np.ndarray, ax = None):

  ax[0].set_xlabel('t')
  ax[0].set_xlim([0, 5])
  ax[0].set_ylim([-1.1, 1.1])
  ax[1].set_xlabel('t')
  ax[1].set_xlim([0, 5])
  ax[1].set_ylim([-1.1, 1.1])
  ax[0].axhline(0, c='gray', ls='--', alpha=0.5)
  ax[1].axhline(0, c='gray', ls='--', alpha=0.5)
  ax[0].plot(t, cf.real, c='b', ls='--', label='Real')
  ax[1].plot(t, cf.imag, c='r', ls='--', label='Imag')
  t, f = ft(x)
  ax[0].scatter(t, f.real, fc='none', ec='b')
  ax[1].scatter(t, f.imag, fc='none', ec='r')
  ax[0].legend(fontsize=12)
  ax[1].legend(fontsize=12)

  circle = plt.Circle((0, 0), 1, color='k', fill=False)
  ax[2].axhline(0, c='gray', ls='--', alpha=0.5)
  ax[2].axvline(0, c='gray', ls='--', alpha=0.5)
  ax[2].set_xlabel('Re')
  ax[2].set_ylabel('Im', labelpad=-6)
  ax[2].set_xlim([-1.1, 1.1])
  ax[2].set_ylim([-1.1, 1.1])
  ax[2].add_patch(circle)
  ax[2].scatter(f.real, f.imag, fc='none', ec='m')
  ax[2].plot(cf.real, cf.imag, c='m', ls='--')


def fit_test(figname: str):

  alpha = 2.0
  beta  = 0.0
  gamma = 1.0
  delta = 0.0
  dt = 0.1
  t = np.arange(0, 100, dt)
  cf = char(t, [alpha, beta, gamma, delta])
  x = generate(n = 100000, param = [alpha, beta, gamma, delta])
  fig, ax = plt.subplots(1, 3, figsize=(12, 4))
  plot_cf(t, cf, x, ax)
  fig.tight_layout()
  fig.savefig(f'fit-char{figname}.png')

  fig, ax = plt.subplots(1, 2, figsize=(8, 4))
  alpha_est, beta_est, gamma_est, delta_est = fit(x, ax)
  fig.tight_layout()
  fig.savefig(f'fit-stable{figname}.png')

  print(alpha, beta, gamma, delta)
  print(alpha_est, beta_est, gamma_est, delta_est)

if __name__ == '__main__':
  fit_test(figname='-levy')
  fit_test(figname='-gaussian')

Koutrouvelisの方法

第一近似では精度不十分な場合や個別対応が必要な場合がありそうです. また, 上の例で見たように, サンプル依存性や波数の分点依存性などフィッティング精度に直結しうる点がいくつかあります. また, $\gamma = 1, \delta = 0$という標準化がされていない場合には, データに乗ったノイズと相まってフィッティングの式の$y$切片を正しく推定できず線形回帰の精度が悪くなる可能性があります. そこで, 以下の流れの反復法を説明します.

1. $\gamma, \delta$の初期推定値$\gamma_0, \delta_0$を与える
2. $\gamma_0, \delta_0$を用いてデータを標準化する
3. 新たに波数の分点を取り直す
4. $\delta_c$を$s^{j} = \frac{s^{j-1}}{\gamma} - \delta_c$のECFの偏角が連続になるようにする
5. 1. で求めた$\delta_c$に対する$s^{j}$からECFを計算する
6. 絶対値や偏角をフィットして$\alpha, \beta, \gamma, \delta$の推定値を更新する
7. 2に戻る; $j=j+1$

現在実装中です, 余裕がでたら実装完成させます。

参考文献, リンク

今回採用した安定分布の定義の原典
[1] Janicki, Weron. Statist. Sci. 9(1): 109-126 (1994)

サンプリング方法: Weronの方法

原典: Theorem 3.1 on
[2] Weron, Rafał. Stat prob lett 28.2 (1996): 165-171.
[3] Weron, Ralf. Rapport technique HSC/95/01, Hugo Steinhaus Center, Wroclaw University of Technology (1995). pdfへの直リンクです

Weron方法の紹介(日本語), 端的にアルゴリズムを述べており, お勧めしたいです. 今回の安定分布の定義と記法が一貫しています.
[4] https://yonesuke.github.io/posts/stable_simulations/

Weronの方法については, 90年代が安定分布が盛んに研究されていた時期だったというのもあってか, 同時期に複数の論文が書かれています.

フィッティング方法: Koutrouvelisの方法

データから経験的な特性関数を求め, 解析式を用いて線型回帰する回帰型と呼ばれる方法. 反復法である.

原典
[5] I. A. Koutrouvelis, J. Am. Stat. Assoc. 75, 918 (1980).

この原典ではECFの原点(0, 0)付近の振る舞いや揺らぎについて議論しており, 波数の分点の取り方の重要性の理由づけを行っている.

フィッティング方法のレビュー
[6] Akgiray, Vedat, and Lamoureux. J Bus Econ Stat 7.1 (1989): 85-93.

上の文献[6]の後にもフィッティング方法がいくつか提案されているようです.

FFT周り

scipyのFFTの利用方法について, 以下の記事で勉強しました.
https://qiita.com/hete_mura/items/4a1eb63b93ff4b879e30
https://qiita.com/kaityo256/items/64a54bb2e2c477cc6fa1

FFTの利用について慣れていなかったので苦戦しました. 特に, フーリエ変換の定義の違いを修正するためにどういう補正が必要かや, $\int f(x) e^{itx} dx$を求めるために$\int f(x) e^{-itx} dx$を計算するルーチンをどう利用するかなど適切な変数変換を理解するために時間を使って考える必要があると思いました.