ほくそ笑む

R言語と統計解析について

人生で役立つ Grimshaw’s Trick

概要

人生を生きていると一般化パレート分布のパラメータを最尤推定したいときがあります。しかしこの推定は2変数の非線形最適化問題を解く必要があり数値計算的に不安定なため人生も不安になります。そんなときにGrimshaw’s Trickを使えばこの問題を一元方程式の解を求める問題に帰着でき数値的に安定するので人生も安心です。本記事ではそんな人生に役立つGrimshaw’s Trickを紹介します。

1. はじめに

極値理論は確率変数の極値 (最大値または最小値) に関する理論である。 極値理論では1920年から1940年にかけて中心極限定理に似た美しい結果が得られた。 それは確率変数列の極値はその分布によらず3つのタイプの確率分布 (極値分布) に弱収束するというものである (Fisher–Tippett–Gnedenko定理)。 これは極値理論の大定理となったが、実データに対して極値分布のパラメータを推定するのは困難であった。 原因は期間ごとの極値をサンプルとして扱う必要があるため必然的にサンプルサイズが小さくなることに加えてパラメータが元の分布に依存することにある。 この問題を解決するために Peaks-Over-Threshold (POT) というアプローチが使われるようになった。 このアプローチは極値が従う分布を調べる代わりにピーク (設定した閾値を超える値) が従う確率分布を調べるというものである。 1970年代に POT アプローチにおける重要な定理が証明された。 ピークが従う確率分布は元の分布によらず一般化パレート分布で近似できるというものである (Pickands–Balkema–de Haan定理)。 この定理により極値理論における次の関心事は一般化パレート分布のパラメータをどうやって推定するのかに移った。 この問題に対してモーメント法によるいくつかの推定法が提案されたが、推定量が理論的に良い性質をもつ最尤法による推定手法が望まれた。 しかし、一般化パレート分布の最尤推定のためには2変数の非線形最適化問題を解く必要があり数値計算の安定性の問題をはらむ。 本記事ではより安定した数値計算のために Grimshaw [1993] で提案された変数変換とその解法を紹介する。

2. 一般化パレート分布の最尤推定

一般化パレート分布は1変数の連続型の確率分布であり、次の確率密度関数をもつ。

$$ f(x) = \begin{cases} \frac{1}{\sigma} \left( 1 + \frac{\gamma}{\sigma}x \right) ^ {-\left(1 + \frac{1}{\gamma}\right)} & \gamma \neq 0 \\ \frac{1}{\sigma} \exp\left( -\frac{x}{\sigma} \right) & \gamma = 0 \end{cases} $$

ただし、 \sigma \in (0, \infty) は尺度パラメータ、 \gamma \in (-\infty, \infty) は形状パラメータである。

 \gamma = 0 の場合、 \lambda = \frac{1}{\sigma} の指数分布となるので、 \sigma の最尤推定値はデータの平均値である。

 \gamma \neq 0 の場合、データ  X = { x_1, \dots, x_n } が与えられたとき、対数尤度関数は次になる。

$$ \begin{align} LogLik(\sigma, \gamma) &= \sum_{i=1}^n \log\left( \frac{1}{\sigma} \left( 1 + \frac{\gamma}{\sigma}x_i \right) ^ {-\left(1 + \frac{1}{\gamma}\right)} \right) \\ &= \sum_{i=1}^n -\log \sigma + \sum_{i=1}^n -\left(1 + \frac{1}{\gamma}\right) \log \left( 1 + \frac{\gamma}{\sigma}x_i \right) \\ &= -n \log \sigma - \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \log \left( 1 + \frac{\gamma}{\sigma}x_i \right)\\ \end{align} $$

この対数尤度関数を最大にするような  \sigma \gamma を求めることが本記事の目的である。

ここで、パラメータ  \gamma \sigma のとる範囲を考えてみよう。 対数尤度の式の  \log \left(1 + \frac{\gamma}{\sigma} x_i\right) より、すべての  x_i \in X に対して  1 + \frac{\gamma}{\sigma} x_i > 0 である必要がある。 これは  \gamma > 0 ならば自然に満たす。  \gamma  \lt  0 の場合、データ  X の最大値で成り立てば他も自然に満たすため、

$$ \begin{align} & \sigma > -\gamma \max(X) \\ \end{align} $$

という条件になる。したがって、パラメータのとる領域は次図のように制限される。

f:id:hoxo_m:20181225205825p:plain

右図に  \gamma = 0.25,  \sigma = 2 のときの対数尤度を示した。 最適なパラメータは境界付近にあり、境界付近では勾配が急であることがわかる。 このような場合に勾配を使った数値最適化は不安定になりやすい。 また、後述するように  \gamma = 0 を特別な場合として扱う必要があるという問題もある。

1980年代に一般化パレート分布に対する最尤推定アルゴリズムは多数提案された (DuMouchel [1983], Davison [1984], R. L. Smith [1984, 1987], J. A. Smith [1986], Joe [1984]) が、どれも数値安定性の問題をはらんでいた。 Grimshaw [1993] は、2パラメータの最適化を行うのではなく、変数変換により一元方程式の求根問題に帰着することで数値的に安定した解を得る方法を示した。

3. Grimshaw's Trick

対数尤度関数を最大化するために、勾配が 0 となる  \sigma \gamma を求めたい ( \nabla LogLik(\sigma, \gamma) = 0)。 まずはそれぞれの偏微分を導出しよう。

$$ \begin{align} \frac{\partial LogLik(\sigma, \gamma)}{\partial \sigma} &= \frac{\partial}{\partial \sigma} \left( -n \log \sigma - \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \log \left( 1 + \frac{\gamma}{\sigma}x_i \right) \right) \\ &= -n \frac{1}{\sigma} - \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \frac{1}{1 + \frac{\gamma}{\sigma}x_i} \left( - \frac{\gamma}{\sigma^2} x_i\right) \\ &= - \frac{n}{\sigma} + \left(1 + \frac{1}{\gamma}\right) \frac{1}{\sigma} \sum_{i=1}^n \frac{\frac{\gamma}{\sigma}x_i}{1 + \frac{\gamma}{\sigma}x_i} \\ &= - \frac{n}{\sigma} + \left(1 + \frac{1}{\gamma}\right) \frac{1}{\sigma} \sum_{i=1}^n \frac{1 + \frac{\gamma}{\sigma}x_i - 1}{1 + \frac{\gamma}{\sigma}x_i} \\ &= - \frac{n}{\sigma} + \left(1 + \frac{1}{\gamma}\right) \frac{1}{\sigma} \sum_{i=1}^n 1 + \left(1 + \frac{1}{\gamma}\right) \frac{1}{\sigma} \sum_{i=1}^n \frac{-1}{1 + \frac{\gamma}{\sigma}x_i} \\ &= - \frac{n}{\sigma} + \left(1 + \frac{1}{\gamma}\right) \frac{n}{\sigma} - \frac{1}{\sigma} \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \frac{1}{1 + \frac{\gamma}{\sigma}x_i} \\ &= \frac{n}{\gamma\sigma} - \frac{1}{\sigma} \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \frac{1}{1 + \frac{\gamma}{\sigma}x_i} \\ \frac{\partial LogLik(\sigma, \gamma)}{\partial \gamma} &= \frac{\partial}{\partial \gamma} \left( -n \log \sigma - \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \log \left( 1 + \frac{\gamma}{\sigma}x_i \right) \right) \\ &= \frac{1}{\gamma^2} \sum_{i=1}^n\log \left( 1 + \frac{\gamma}{\sigma}x_i \right) - \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \frac{\frac{1}{\sigma}x_i}{1 + \frac{\gamma}{\sigma}x_i} \\ &= \frac{1}{\gamma^2} \sum_{i=1}^n\log \left( 1 + \frac{\gamma}{\sigma}x_i \right) - \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \frac{1}{\gamma} \frac{1 + \frac{\gamma}{\sigma}x_i - 1}{1 + \frac{\gamma}{\sigma}x_i} \\ &= \frac{1}{\gamma^2} \sum_{i=1}^n\log \left( 1 + \frac{\gamma}{\sigma}x_i \right) - \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \frac{1}{\gamma} - \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \frac{1}{\gamma} \frac{- 1}{1 + \frac{\gamma}{\sigma}x_i} \\ &= \frac{1}{\gamma^2} \sum_{i=1}^n\log \left( 1 + \frac{\gamma}{\sigma}x_i \right) - \left(1 + \frac{1}{\gamma}\right) \frac{n}{\gamma} + \frac{1}{\gamma} \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \frac{ 1}{1 + \frac{\gamma}{\sigma}x_i} \\ \end{align} $$

したがって、一般化パレート分布の最尤推定を行うには次の連立方程式を解けばよい。

$$ \begin{cases} \begin{align} \frac{n}{\gamma\sigma} - \frac{1}{\sigma} \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \frac{1}{1 + \frac{\gamma}{\sigma}x_i} &= 0 & (1)\\ \frac{1}{\gamma^2} \sum_{i=1}^n\log \left( 1 + \frac{\gamma}{\sigma}x_i \right) - \left(1 + \frac{1}{\gamma}\right) \frac{n}{\gamma} + \frac{1}{\gamma} \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \frac{ 1}{1 + \frac{\gamma}{\sigma}x_i} &= 0 & (2)\\ \end{align} \end{cases} $$

実はこの偏微分の導出過程ですでに Grimshaw のアイデアが使われている。 連立方程式の左辺の後半部分が同じ形になるように式変形を行った。 式 (1) に  \sigma をかけ、式 (2) に  \gamma をかけるとこの部分が完全に同じになる。

$$ \begin{cases} \begin{align} \frac{n}{\gamma} - \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \frac{1}{1 + \frac{\gamma}{\sigma}x_i} &= 0\\ \frac{1}{\gamma} \sum_{i=1}^n\log \left( 1 + \frac{\gamma}{\sigma}x_i \right) - n - \frac{n}{\gamma} + \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \frac{ 1}{1 + \frac{\gamma}{\sigma}x_i} &= 0\\ \end{align} \end{cases} $$

上式と下式を足し合わせることで次の式が得られる。

$$ \begin{align} & \frac{1}{\gamma} \sum_{i=1}^n\log \left( 1 + \frac{\gamma}{\sigma}x_i \right) - n = 0 &\\ & \gamma = \frac{1}{n} \sum_{i=1}^n\log \left( 1 + \frac{\gamma}{\sigma}x_i \right) & (3) \end{align} $$

式 (3) は最適なパラメータが満たすべき条件となる。

例えば、 \gamma=0.25,  \sigma=2 の場合、式 (3) は次図の青線になる。 この線上に最適パラメータ (赤点) が存在する。  \gamma=0 は必ず条件を満たすが、特別に扱う必要があり、これも上記の最適化問題を難しくしている。

f:id:hoxo_m:20181225205918p:plain

そこで、Grimshaw は  \theta = \frac{\gamma}{\sigma} による変数変換を提案した。  (\gamma, \sigma) (\gamma, \theta) に写すことで次の良い性質が得られる。

  • 変換  (\gamma, \sigma) \rightarrow (\gamma, \theta) \gamma = 0 を除いて全単射である。
  •  \gamma=0 ならば  \theta=0 であり、特別な場合が1点に集約される。
  • 解の領域がパラメータに依存せず矩形になる ( \theta > -\frac{1}{\max(X)})。
  • 勾配が急になる境界は  \theta = -\frac{1}{\max(X)} 付近のみである。

f:id:hoxo_m:20181225205942p:plain

この変換により、最適条件を満たす場合、式 (3) から  \gamma \theta によって定まる。

$$ \begin{align} \gamma &= \frac{1}{n} \sum_{i=1}^n\log \left( 1 + \theta x_i \right) \end{align} $$

したがって、解を求めるために連立方程式を解く必要はなくなり、最適な  \theta だけを求めれば良い。

具体的には2変数の尤度関数を  \theta だけの関数に変換する。

$$ \begin{cases} \gamma = \frac{1}{n} \sum_{i=1}^n\log \left( 1 + \theta x_i \right) \\ \sigma = \frac{\gamma}{\theta} \end{cases} $$

$$ \begin{align} LogLik(\theta) &= -n \log \sigma - \left(1 + \frac{1}{\gamma}\right) \sum_{i=1}^n \log \left( 1 + \frac{\gamma}{\sigma}x_i \right)\\ &= -n \log \left( \frac{\gamma}{\theta} \right) - \left(1 + \frac{1}{\gamma}\right) n\gamma\\ &= -n \log \left( \frac{1}{\theta} \frac{1}{n} \sum_{i=1}^n\log \left( 1 + \theta x_i \right) \right) - n\gamma -n \\ &= -n \log \left( \frac{1}{n} \sum_{i=1}^n \frac{1}{\theta} \log \left( 1 + \theta x_i \right) \right) - \sum_{i=1}^n\log \left( 1 + \theta x_i \right) -n \\ \end{align} $$

このように、尤度関数を注目したいパラメータだけの関数として表したものをプロファイル尤度という。 一般化パレート分布の最尤推定の問題は、この1パラメータのプロファイル対数尤度を最大化する問題となる。

プロファイル対数尤度の微分を求める。

$$ \begin{align} \frac{d\ LogLik(\theta)}{d\theta} &= -n \frac{ \frac{1}{n} \sum_{i=1}^n \left( -\frac{1}{\theta^2} \log(1+\theta x_i) + \frac{1}{\theta}\frac{x_i}{1+\theta x_i} \right) }{ \frac{1}{n} \sum_{i=1}^n \frac{1}{\theta} \log \left( 1 + \theta x_i \right) } - \sum_{i=1}^n \frac{x_i}{1 + \theta x_i} \\ &= -n \frac{ - \sum_{i=1}^n \frac{1}{\theta^2} \log(1+\theta x_i) + \sum_{i=1}^n \frac{1}{\theta}\frac{x_i}{1+\theta x_i} }{ \sum_{i=1}^n \frac{1}{\theta} \log \left( 1 + \theta x_i \right) } - \sum_{i=1}^n \frac{x_i}{1 + \theta x_i} \\ &= \frac{n}{\theta} -n \frac{\sum_{i=1}^n \frac{x_i}{1+\theta x_i}}{\sum_{i=1}^n \log \left( 1 + \theta x_i \right)} - \left( \frac{n}{\theta} - \frac{1}{\theta} \sum_{i=1}^n \frac{1}{1+\theta x_i} \right)\\ &= -n\frac{\left( \frac{n}{\theta} - \frac{1}{\theta} \sum_{i=1}^n \frac{1}{1+\theta x_i} \right)}{\sum_{i=1}^n \log \left( 1 + \theta x_i \right)} + \frac{1}{\theta} \sum_{i=1}^n \frac{1}{1+\theta x_i} \\ &= -\frac{n}{\theta}\left( n - \sum_{i=1}^n \frac{1}{1+\theta x_i} \right) \left( \frac{1}{\sum_{i=1}^n \log \left( 1 + \theta x_i \right)}\right) + \frac{1}{\theta} \sum_{i=1}^n \frac{1}{1+\theta x_i} \end{align} $$

だだし、次を用いた。

$$ \begin{align} \sum_{i=1}^n \frac{x_i}{1+\theta x_i} &= \frac{1}{\theta} \sum_{i=1}^n \frac{1 + \theta x_i + 1}{1+\theta x_i} \\ &= \frac{1}{\theta} \sum_{i=1}^n 1 + \frac{1}{\theta} \sum_{i=1}^n \frac{-1}{1+\theta x_i} \\ &= \frac{n}{\theta} - \frac{1}{\theta} \sum_{i=1}^n \frac{1}{1+\theta x_i} \end{align} $$

プロファイル対数尤度の微分が 0 になるには次を満たせばよい。

$$ \begin{align} -\frac{n}{\theta}\left( n - \sum_{i=1}^n \frac{1}{1+\theta x_i} \right) \left( \frac{1}{\sum_{i=1}^n \log \left( 1 + \theta x_i \right)}\right) + \frac{1}{\theta} \sum_{i=1}^n \frac{1}{1+\theta x_i} &= 0 \\ \left( n - \sum_{i=1}^n \frac{1}{1+\theta x_i} \right) \left( \frac{1}{\sum_{i=1}^n \log \left( 1 + \theta x_i \right)}\right) - \frac{1}{n} \sum_{i=1}^n \frac{1}{1+\theta x_i} &= 0 \\ \left( n - \sum_{i=1}^n \frac{1}{1+\theta x_i} \right) - \frac{1}{n} \sum_{i=1}^n \frac{1}{1+\theta x_i} \left( \sum_{i=1}^n \log \left( 1 + \theta x_i \right)\right) &= 0 \\ n - \left( \sum_{i=1}^n \frac{1}{1+\theta x_i} \right) \left( 1 + \frac{1}{n} \sum_{i=1}^n \log \left( 1 + \theta x_i \right)\right) &= 0 \\ \left( \frac{1}{n}\sum_{i=1}^n \frac{1}{1+\theta x_i} \right) \left( 1 + \frac{1}{n} \sum_{i=1}^n \log \left( 1 + \theta x_i \right)\right) &= 1 \\ \end{align} $$

ここで

$$ \begin{align} u(\theta) &= \frac{1}{n}\sum_{i=1}^n \frac{1}{1+\theta x_i} \\ v(\theta) &= 1 + \frac{1}{n} \sum_{i=1}^n \log \left( 1 + \theta x_i \right) \end{align} $$

とおくと

$$ \begin{align} u(\theta) v(\theta) &= 1 & (4) \end{align} $$

と表せる。

この方程式を満たす  \theta を求めることで、一般化パレート分布の2つのパラメータを次の式で計算できる。

$$ \begin{align} \gamma &= \frac{1}{n} \sum_{i=1}^n \log \left( 1 + \theta x_i \right) \\ \sigma &= \frac{\gamma}{\theta} \end{align} $$

これが Grimshaw's Trick である。

4. 解の存在

Grimshaw's Trick により一般化パレート分布の最尤推定は式 (4) の方程式を解く問題に帰着された。

ここで

$$ h(\theta) = u(\theta) v(\theta) - 1 $$

とおくと、 h(\theta) = 0 の解を求めればよい。

次図に  \gamma=0.25,  \sigma = 2 のときのプロファイル対数尤度と  h(\theta) を示す。

f:id:hoxo_m:20181225210005p:plain

方程式  h(\theta) = 0 に解が存在することを示そう。

まず  \theta = 0 ならば  u(\theta) = 1,  v(\theta)=1 なので  h(\theta) = 0 となる。 したがって  \theta = 0 は常に解である。 ただし上図からわかるように  \theta = 0 はプロファイル尤度を最大にするとは限らない。

 \theta = 0 以外に解が存在することは次で示すことができる。

$$ \begin{align} \theta \rightarrow \frac{-1}{\max(X)} &\ \Rightarrow\ h(\theta) \rightarrow \infty & (a)\\ \theta > 2\frac{\bar{X} - \min(X)}{\min(X)^2} &\ \Rightarrow\ h(\theta) < 0 & (b) \\ \theta \rightarrow 0 &\ \Rightarrow \ h'(\theta) = 0 & (c) \end{align} $$

ただし  \bar{X} はデータ  X の平均値である。

これらの結果より、 h(\theta) \theta = 0 で鞍点でなければ次の範囲内に  h(\theta) = 0 となる  \theta \theta = 0 以外に1つ以上あることがわかる。

$$ \begin{align} &-\frac{1}{\max(X)} < \theta < 2\frac{\bar{X} - \min(X)}{\min(X)^2} & (5) \\ \end{align} $$

式 (a) - (c) の証明については Grimshaw [1993] では次のように書かれている。

Proof. The results are straightforward.
(訳:この証明は非常に難しい)

したがって、本記事では証明しないことにする1

一方、 h(\theta) = 0 の解が1つである保証はない。 そのため次の手順でプロファイル対数尤度が最大となる  \theta を見つける。

  1. 方程式  h(\theta) = 0 のすべての解を求める。
  2. すべての解に対して  LogLik(\theta) を計算し、最大になる  \theta を選ぶ。

ただし、複数の根のある方程式の解法は一般には知られていない2。 Grimshaw [1993] では方程式  h(\theta) = 0 に対して、解の探索区間を細かく区切り、ニュートン・ラフソン法を繰り返し用いることですべての解を求める手順が示されている。 この手順はR言語では mev パッケージの gp.fit(method="Grimshaw") で実行できる。

# データの生成
sigma <- 2
gamma <- 0.25

set.seed(314)
x <- POT::rgpd(200, 0, sigma, gamma)

# Grimshaw の手順でパラメータを最尤推定
library(mev)
fit <- gp.fit(x, 0, method="Grimshaw")
est_sigma <- unname(fit$estimate[1])
est_gamma <- unname(fit$estimate[2])

# 推定結果の表示
c(sigma = est_sigma, gamma = est_gamma)
##     sigma     gamma 
## 1.9039251 0.2957402

5. Sifferらの求根法

方程式  h(\theta) = 0 のすべての解を求める方法として、Siffer et al. [2017] で提案された方法が面白いので紹介したい。

式 (5) より、解の探索範囲として次の2つの区間を考える。ただし  \varepsilon は非常に小さい正の実数とする。

$$ \begin{align} I_1 &= \left[-\frac{1}{\max(X)} + \varepsilon,\ -\varepsilon \right] \\ I_2 &= \left[\varepsilon,\ 2\frac{\bar{X} - \min(X)}{\min(X)^2} - \varepsilon\right] \\ \end{align} $$

この2つの区間において、次で定義する  K 変数関数の最小化を行う。

$$ f(\theta_1, \dots, \theta_K) = h(\theta_1)^2 + \dots h(\theta_K)^2 $$

すなわち、次の最小化問題を解く。

$$ \min_{\theta_1, \dots, \theta_K} \sum_{k=1}^K h(\theta_k)^2 $$

ただし、 \theta_1, \dots, \theta_K の初期値をそれぞれの区間で散らばるように設定する。 これにより  \theta_k h(\theta) = 0 の異なる局所解に収束することが期待される。 したがって、この最小化問題を解くことで、最大で  K 個の根が得られる。

2つの区間それぞれに対して最小化問題を解き、得られた  2K 個の根に 0 (常に解) を加える。 こうして得られた  2K + 1 個の根に対してプロファイル尤度を計算し、尤度が最大となる根を選択するというのが Sifferらの提案した方法である。

Siffer et al. [2017] では  K=10,  \varepsilon = 10^{-8}、最適化手法として L-BFGS-B法が用いられている。

R言語では Sifferらの方法は次のように実装できる。

#' 一般化パレート分布の最尤推定を Sifferらの方法で行う関数
#'
#' @param x データを数値ベクトルとして入力
#' @param K デフォルトは 10
#' @param epsilon 境界を探索範囲から外すための微小値。デフォルトは 10^-8
#' @return パラメータの最尤推定値 c(sigma, gamma) 
gp.fit_siffer <- function(x, K = 10, epsilon = 10^-8) {
  # 式 (4) の方程式
  u <- function(theta) mean(1 / (1 + theta * x))
  v <- function(theta) 1 + mean(log(1 + theta * x))
  h <- function(theta) u(theta) * v(theta) - 1
  
  # 最適化目的関数
  fn <- function(theta_K, logged=FALSE) {
    if(logged) theta_K <- exp(theta_K)
    sum(vapply(theta_K, h, double(1))^2)
  }
  
  # 解の範囲
  lower_bound <- -1 / max(x)
  upper_bound <- 2 * (mean(x) - min(x)) / min(x)^2

  theta_range1 <- c(lower_bound + epsilon, -epsilon)
  theta_range2 <- c(epsilon, upper_bound)

  # 初期値の設定
  theta_init1 <- seq(lower_bound, 0, length.out = K+2)[-c(1, K+2)]
  theta_init2 <- seq(log(epsilon), log(upper_bound), length.out = K+2)[-c(1, K+2)]

  # 最適化の実行
  opt1 <- optim(theta_init1, fn = fn, method = "L-BFGS-B",
                lower = theta_range1[1], upper = theta_range1[2])
  opt2 <- optim(theta_init2, fn = fn, logged = TRUE, method = "L-BFGS-B",
                lower = log(theta_range2[1]), upper = log(theta_range2[2]))

  # 最適解 (推定値の候補)
  opt_theta_K <- unique(c(opt1$par, 0, exp(opt2$par)))

  # プロファイル対数尤度
  n <- length(x)
  profileLogLik <- function(theta) {
    -n * log(mean(log(1 + theta * x))/theta) - n * mean(log(1 + theta * x)) - n
  }

  # 最適解のうちプロファイル対数尤度が最大のものを推定値とする
  pLLs <- vapply(opt_theta_K, profileLogLik, double(1))
  est_theta <- opt_theta_K[which.max(pLLs)]
  
  # パラメータの引き戻し
  est_gamma <- mean(log(1 + est_theta * x))
  est_sigma <- est_gamma / est_theta
  
  # 結果の返却
  c(sigma = est_sigma, gamma = est_gamma)
}

推定を実行してみる。

# データの生成
sigma <- 2
gamma <- 0.25

set.seed(314)
x <- POT::rgpd(200, 0, sigma, gamma)

# Sifferらの方法でパラメータを最尤推定
gp.fit_siffer(x)
##     sigma     gamma 
## 1.9040107 0.2956953

Grimshaw の手順と同様の結果が得られた。

6. まとめ

「人生を安定させるにはまず数値計算を安定させよ」という格言があります3。 本記事では人生を安定させるための第一歩として、一般化パレート分布の最尤推定を安定させるための Grimshaw's Trick を紹介しました。 あなたの人生で一般化パレート分布を最尤推定したい、そんなときに役立つ知識が得られたのではないかと思います。

参考文献

  • Scott D. Grimshaw (1993) Computing Maximum Likelihood Estimates for the Generalized Pareto Distribution, Technometrics, 35:2, 185-191, DOI: https://doi.org/10.1080/00401706.1993.10485040
  • Alban Siffer, Pierre-Alain Fouque, Alexandre Termier, and Christine Largouet. 2017. Anomaly Detection in Streams with Extreme Value Theory. In Proceedings of the 23rd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD '17). ACM, New York, NY, USA, 1067-1075. DOI: https://doi.org/10.1145/3097983.3098144

付録

式 (b) を証明する。

$$ \begin{align} \theta > 2\frac{\bar{X} - \min(X)}{\min(X)^2} &\ \Rightarrow\ h(\theta) \lt 0 & (b) \\ \end{align} $$

イエンセンの不等式より

$$ \begin{align} v(\theta) &= 1 + \frac{1}{n} \sum_{i=1}^n \log (1 + \theta x_i)\\ &\leq 1 + \log(1 + \theta \bar{X}) \end{align} $$

また、任意の  x_i に対して  \min(X) \leq x_i より

$$ \begin{align} u(\theta) &= \frac{1}{n} \sum_{i=1}^n \frac{1}{1 + \theta x_i} \\ &\leq \frac{1}{n} \sum_{i=1}^n \frac{1}{1 + \theta \min(X)} \\ &\leq \frac{1}{1 + \theta \min(X)} \end{align} $$

したがって、

$$ \begin{align} h(\theta) &= u(\theta)v(\theta) - 1 \\ &\leq \frac{1 + \log(1 + \theta \bar{X})}{1 + \theta \min(X)} - 1 \\ &= \frac{1 + \log(1 + \theta \bar{X}) - 1 - \theta \min(X)}{1 + \theta \min(X)} \\ &= \frac{\log(1 + \theta \bar{X}) - \theta \min(X)}{1 + \theta \min(X)} \\ \end{align} $$

分子が負ならば  h(\theta)  \lt  0 である。よって

$$ \begin{align} \frac{2(\bar{X} -\min(\theta))}{\min(X)^2} < \theta &\ \Rightarrow\ \bar{X} < \min(X) + \frac{1}{2}\theta\min(X)^2 \\ &\ \Rightarrow\ 1 + \theta \bar{X} < 1 + \theta\min(X) + \frac{1}{2}\theta^2\min(X)^2 \\ &\ \Rightarrow\ 1 + \theta \bar{X} < \exp(\theta\min(X)) \\ &\ \Rightarrow\ \log (1 + \theta \bar{X}) < \theta\min(X) \\ &\ \Rightarrow\ h(\theta) < 0 \\ \end{align} $$

以上で証明された。ただし、指数関数のマクローリン展開

$$ \begin{align} \exp(x) &= 1 + x + \frac{1}{2!} x^2 + \frac{1}{3!} x^3 + \dots \\ &< 1 + x + \frac{1}{2} x^2 \end{align} $$

をもちいた。


  1. これはウィットに富んだジョークである。(b) の証明は付録に示す。

  2. 多項式の場合はスツルム (Sturm) の解法が知られている。

  3. 私が考えた格言です。

非心ベータ分布を Stan で推定する (1)

1. はじめに

ベータ分布は範囲が固定されている連続データの分布を柔軟に表現できるため便利である。 非心ベータ分布 (Noncentral Beta Distribution) はベータ分布を一般化した確率分布であり柔軟性はさらに上がる。

R ではベータ分布に従う乱数を生成する rbeta() 関数に ncp 引数1を指定することで非心ベータ分布に従う乱数を生成できる。

x <- rbeta(2000, 10, 5, ncp = 3)

hist(x, probability = TRUE)
curve(dbeta(x, 10, 5), add = TRUE, col = "red")
curve(dbeta(x, 10, 5, ncp = 3), add = TRUE, col = "blue")

f:id:hoxo_m:20180823223727p:plain

赤線がベータ分布、青線が非心ベータ分布であり、ヒストグラムへの当てはまりは青線の方が良い。

非心ベータ分布はベータ分布の中心をずらした分布になる。 右にずらす場合は Type I、左にずらす場合は Type II と呼ばれる。 rbeta() は Type I にしか対応していないが、Type I と Type II は単に左右反転しただけなので 1 - x を取ることで Type II を表せる。

本記事では Stan を使って非心ベータ分布のパラメータを推定する方法について述べる。

2. 対数尤度関数

Stan には非心ベータ分布を推定するための組み込み関数はない。 したがって、target += 記法を用いてモデルを記述する。 そのために、非心ベータ分布の対数尤度関数を求める必要がある。

定義より、非心ベータ分布の密度関数は次である。

\begin{align} \mathrm{NCBeta}(x \mid \alpha, \beta, \lambda) = \sum_{j = 0}^\infty \mathrm{Pois}(j \mid \lambda) \mathrm{Beta}(x \mid \alpha + j, \beta) \end{align}

ただし、Pois はポアソン分布の密度関数、Beta はベータ分布の密度関数である。

ここから尤度関数を求めたいが、 j についての無限和の存在が問題になる。 ここで、ベータ分布の密度関数を  j に対するオーダーで書き直すと

\begin{align} \mathrm{Beta}(x \mid \alpha + j, \beta) &= \frac{x^{\alpha+j-1}(1-x)^{\beta-1}}{\mathrm{B}(\alpha + j, \beta)} \\ &= \frac{O(x^j)}{O\left(\frac{(\alpha+j-1)!(\beta - 1)!}{(\alpha + j + \beta - 1)!}\right)} \\ &= \frac{O(x^j)}{O\left(\frac{1}{j^\beta}\right)} \\ &= O(x^j j^\beta) \end{align}

となる。ただし B はベータ関数であり、そのオーダーとして  \alpha \beta が整数の場合を考えた。

一方、ポアソン分布の密度関数の  j に対するオーダーは次のようになる。

\begin{align} \mathrm{Pois}(j \mid \lambda) &= \frac{\left(\lambda/2\right)^j}{j!} e^{-\lambda} \\ &= \frac{O(\lambda^j)}{O(j!)} \end{align}

すなわち、ベータ分布は  x = 1 のときでも高々  O(j^\beta) でしか増加しないのに対し、ポアソン分布は  O(j!) の速度で減少する。

したがって、非心ベータ分布の密度関数の  j に関して次のことが言える。

  •  j が十分大きいとき項は 0 とみなせる。
  • ポアソン分布の重みが支配的である。

すなわち、0 とみなせる  j についてはポアソン分布の重みのみを考えれば良い。

ポアソン分布の相対的な重みを比較するには次のようにする。

lambda <- 7
j <- rpois(10000, lambda / 2)
round(table(j) / max(table(j)), 2)
j
   0    1    2    3    4    5    6    7    8    9   10   11   12   13 
0.14 0.49 0.81 1.00 0.82 0.58 0.37 0.18 0.08 0.03 0.01 0.00 0.00 0.00 

これより、 \lambda \leq 7 に限れば、 j は 0 から 10 までを考えれば十分である。  \lambda の上限を上げると  j の範囲も広がり、そのぶん計算にかかる時間は増加することに注意。

以上によりベータ分布の対数尤度関数は次のように表される。

\begin{align} \log\mathrm{NCBeta}(x \mid \alpha, \beta, \lambda) &= \log\left( \sum_{j = 0}^{10} \mathrm{Pois}(j \mid \lambda) \mathrm{Beta}(x \mid \alpha + j, \beta) \right) \\ &= \log \sum_{j=0}^{10} \exp\left( \log \mathrm{Pois}(j \mid \lambda) + \log \mathrm{Beta}(x \mid \alpha + j, \beta) \right) \end{align}

これは log_sum_exp の形になっていることに注意しよう。

3. Stan による実装

非心ベータ分布の対数尤度関数がわかったので、これを Stan で実装すると次のようになる。

data {
  int<lower=1> N;
  real<lower=0, upper=1> values[N];
}
parameters {
  real<lower=1> alpha;
  real<lower=1> beta;
  real<lower=0, upper=7> lambda;
}
model {
  real logliks[11];
  for (i in 1:N) {
    for (j in 0:10) {
      logliks[j + 1] = poisson_lpmf(j | lambda / 2);
      logliks[j + 1] += beta_lpdf(values[i] | alpha + j, beta);
    }
    target += log_sum_exp(logliks);
  }
}

これを ncbeta01.stan という名前でファイルに保存する。

推定を実行する R コードは以下である。

set.seed(314)

alpha <- 10
beta <- 5
ncp <- 3
N <- 200

x <- rbeta(N, alpha, beta, ncp)

standata <- list(N = N, values = x) # Type I 
# standata <- list(N = N, values = 1 - x) # Type II

library(rstan)
# options(mc.cores = 3)
model <- stan_model("ncbeta01.stan")
fit <- sampling(model, data = standata, chains = 3)

print(fit)
         mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
alpha    8.93    0.03 1.15   6.70   8.18   8.91   9.70  11.29  1160    1
beta     4.57    0.01 0.42   3.80   4.28   4.56   4.86   5.43  1491    1
lambda   3.17    0.06 1.97   0.17   1.45   3.04   4.81   6.73  1255    1
lp__   150.20    0.04 1.25 147.12 149.60 150.52 151.13 151.65   928    1

真の値  \alpha = 10,  \beta = 5,  \lambda = 3 がうまく推定できたようである。

まとめ

非心ベータ分布を Stan で推定する方法を述べた。 上記のコードだと私の PC では 1 chain 80秒程度かかる。

次回はこれを高速化して 1 chain 10秒程度にする方法を書こうと思う。

StanとRでベイズ統計モデリング (Wonderful R)

StanとRでベイズ統計モデリング (Wonderful R)


  1. npc は Non-Centrality Parameter の略

施策効果測定におけるメタアナリシスの応用

はじめに

マーケティング施策を行うときに、その施策効果を測定するために、コントロールグループ(施策を適用しないユーザ)を作る場合がある。 例えば、販促メールを送るという施策を行うときに、一部のユーザには送らないようにする。 仮にメールを送らなかったユーザの平均売上が1000円であり、メールを送ったユーザは平均1100円だとすると、その差である100円が1人あたりの施策効果となる。

しかし、施策を適用しないユーザの数を増やすと、全体の売上効果はそれだけ減少してしまう。 つまり、全体のユーザ数を1万人とすると、全員にメールを送れば100万円の売上効果があるが、半分の5千人をコントロールとすると売上効果も半分の50万円となる。 したがって、コントロールグループに割り当てる人数はなるべく小さくしたいという要求がある。

ただし、コントロールグループの人数を少なくすると、効果測定の精度が落ちるという問題がある。 今、1万人のユーザがいるECサイトで、メール施策を行うことを考えよう。 このとき、コントロールグループは全体の 1% である 100人しか取れないとする。 ユーザの作り出す売上がポアソン分布に従うとして、

  • メールを送らないユーザは平均 1000 円の売上
  • メールを送ると平均 1100 の売上

すなわちメール施策の真の効果は 100円としてシミュレーションしてみる。

set.seed(17744)
# 施策未実施ユーザ 100人の個々人の売上
control_group   <- rpois( 100, 1.0) * 1000
# 施策実施ユーザの 9900人の個々人の売上
treatment_group <- rpois(9900, 1.1) * 1000

# それぞれのグループに対して10人の売上を表示
head(control_group, 10)
#> [1] 1000 2000 2000 2000 1000 1000    0    0    0 4000
head(treatment_group, 10)
#> [1] 4000 1000 3000 1000    0 1000 2000 1000 1000 1000

こうして得られた各ユーザの売上を使って、施策効果を推定する。 各グループの平均売上の差が施策効果である。

effect <- mean(treatment_group) - mean(control_group)
effect
#> [1] 47.3

施策効果は 47 円と算出され、真の効果である 100円から大きくずれてしまった。

また、施策に本当に効果があったかは統計的仮説検定を用いて判断することが多い。 ここでは Welch t-test を使ってみよう。

t.test(treatment_group, control_group)
 Welch Two Sample t-test

data:  treatment_group and control_group
t = 0.41425, df = 100.68, p-value = 0.6796
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -179.1100  273.6555
sample estimates:
mean of x mean of y 
 1087.273  1040.000 

検定の結果より、p-value = 0.68 なのでこの施策には有意な効果はないと判断されてしまう。

つまり、コントロールグループの人数が少ないことが原因で

  • 施策効果の推定値が真の値から大きくずれてしまう
  • 検定を行なっても有意な効果がないと判断されてしまう

という問題が起こる。

ただし、このようなマーケティング施策は繰り返し行われることが多い。 例えば、販促メールを毎週水曜日に送るなどである。 この場合、メタアナリシスを用いて上記の問題を解決できる。

今、上記のメール施策を毎週水曜日に5回繰り返したとする。 このときのメタアナリシスの結果を下図に示す。

f:id:hoxo_m:20180809192359p:plain

各週のメール施策において算出された効果は Effect 列に記載されている。 これはばらつきが大きく、マイナスになることもあることがわかる(2018-08-08)。 右側のグラフは、四角が Effect であり、そこから左右に伸びる線が信頼区間を表している。縦線は 0円であり、2018-08-22を除いて各施策は 0 をまたいでいるため、検定結果は有意でない。

メタアナリシスによってこれらの結果を統合することで、推定精度を向上させることができる。 メタアナリシスによる推定結果は上図の Summary に記されている。 これによると、効果の推定値は 114 円であり、真の効果 100円に近い。 また、信頼区間は 26円から 202円となり、0円をまたがないので検定は有意になる。

本記事ではRでこれを行う方法を説明する。

メタアナリシスの実践

次のように、5回の施策に対して、その効果と信頼区間を求める。

# 8/1(水) から 8/29(水) までの水曜日の日付を取得
begin <- as.Date("2018-08-01")
end   <- as.Date("2018-08-29")
dates <- seq(begin, end, by = "weeks")

result <- data.frame()
for (i in seq_along(dates)) {
  date <- dates[i]
  
  # 売上のシミュレーション
  set.seed(as.integer(date))
  control_group   <- rpois( 100, 1.0) * 1000
  treatment_group <- rpois(9900, 1.1) * 1000
  
  # 売上効果の算出
  effect <- mean(treatment_group) - mean(control_group)
  
  # Welch t-test による信頼区間の推定
  t <- t.test(treatment_group, control_group, conf.level = 0.95)
  lower <- t$conf.int[1]
  upper <- t$conf.int[2]

  # # ブートストラップ法による信頼区間の推定
  # library(simpleboot)
  # boot <- two.boot(treatment_group, control_group, mean, R = 1000)
  # ci <- boot.ci(boot, conf = 0.95, type = "perc")
  # lower <- ci$perc[4]
  # upper <- ci$perc[5]

  res <- data.frame(date, effect, lower, upper)
  
  result <- result %>% bind_rows(res)
}

print(result)
        date     effect     lower    upper
1 2018-08-01  47.272727 -171.8910 268.2418
2 2018-08-08  -4.545455 -224.9838 181.0983
3 2018-08-15  18.181818 -171.9945 194.4213
4 2018-08-22 306.969697  150.0565 465.5196
5 2018-08-29  92.323232 -129.1142 297.0604

それぞれの施策実施日について、効果と信頼区間が算出された。 ここでは信頼区間の推定にt検定を用いたが、売上が正規分布に従わないことが不安ならブートストラップ法を使うと良い。 ただし以下では信頼区間は平均について対称であることを仮定している。

この結果に対して『Rで楽しむ統計』(p.108) を参考にメタアナリシスを行う。

Rで楽しむ統計 (Wonderful R 1)

Rで楽しむ統計 (Wonderful R 1)

library(dplyr)

summary <- result %>%
  mutate(weight = 1/(upper - lower)^2, 
         effect = effect * weight / sum(weight)) %>%
  summarise(effect = sum(effect), sd = sqrt(1/sum(weight))/2) %>%
  mutate(lower = effect - sd, upper = effect + sd) %>%
  select(-sd) %>%
  mutate(date = "Summary") %>%
  select(date, everything())

print(summary)
     date   effect    lower    upper
1 Summary 114.2419 26.44446 202.0394

それぞれの施策に対して信頼区間の幅を2乗して逆数を取ったものを重みとする。 それぞれの施策効果に対して重み付けを行い、全てを足し合わせたものが統合された施策効果となる。 また、信頼区間は重みの総和の逆数の平方根によって算出される。 メタアナリシスによって算出された施策効果は 114円であり、真の効果100円と近い値が算出された。 また、信頼区間は 26円から202円であり、0円をまたがないため、統計的に有意な結果であることが言える。 ここではメタアナリシスのモデルとして固定効果モデルを用いたが、一般的には施策効果にはぶれが生じるため、それを考慮したい場合はランダム効果モデルを用いる。

メタアナリシスの可視化にはフォレストプロットというものが使われる。 これは次のように作成できる。

library(forestplot)

nrow <- nrow(result) # 施策の回数

# 施策結果とメタアナリシスの結果を統合
df <- result %>% 
  mutate(date = as.character(date)) %>% 
  bind_rows(summary)

# フォレストプロットの左側に表示される表データの作成(日付と効果)
labeltext <- list(c("Date", as.character(df$date)),
                  c("Effect", round(df$effect)))

# フォレストプロットの描画
forestplot(labeltext, mean = c(NA, df$effect), 
           lower = c(NA, df$lower), upper = c(NA, df$upper),
           is.summary = c(TRUE, rep(FALSE, nrow), TRUE),
           hrzl_lines = gpar(col="#444444"))

f:id:hoxo_m:20180809192359p:plain

それぞれの施策における効果と信頼区間、そして統合された結果が可視化された。

まとめ

マーケティング施策効果の測定において、

  • 施策を実施しないコントロールグループが作られる
  • コントロールグループの人数が少なく効果測定の精度が低い
  • 同じ施策が繰り返される

という場合において、メタアナリシスを用いて推定精度を向上させる方法を説明した。

参考文献

Rで楽しむ統計 (Wonderful R 1)

Rで楽しむ統計 (Wonderful R 1)

Stochastic Online Anomaly Analysis for Streaming Time Series 読んだ

Stochastic Online Anomaly Analysis for Streaming Time Series (IJCAI 2017) https://doi.org/10.24963/ijcai.2017/445

Student-t process を使ったストリーミング時系列に対する新しい異常検知手法を提案する。 時系列関数が Student-t process に従うと仮定し予測分布を推定する。 得られた予測分布に対して実際の観測値が確率 p の予測区間から外れた場合を異常とみなす。 モデルの学習はフルベイズで行われ、カーネル関数の持つハイパーパラメータについても学習する。 (カーネル関数の選択はある程度重要になってくると思われるが論文中では選択法については触れていない) 通常はバッチ学習だがストリーミングデータに対応するために SGD を用いたオンライン学習手法に改変する。 実験では、既存手法より提案手法の方が、異常検知精度と予測精度の両方において良いという結果が得られた。 [Smith+ 2012] ではガウス過程を使った同様の手法が提案されているが、学習データに外れ値が含まれている場合に予測精度が落ちる。 これに対し、Student-t process では外れ値に対してロバストであることを実験により示した。

Anomaly Detection in Streams with Extreme Value Theory 読んだ

Anomaly Detection in Streams with Extreme Value Theory (KDD 2017) http://www.kdd.org/kdd2017/papers/view/anomaly-detection-in-streams-with-extreme-value-theory

ほとんどの異常検知手法は異常スコアを算出するが、異常かどうかを判断するためのしきい値を設定するのは難しい課題である。 本論文ではより直感的なパラメータとして偽陽性率を設定するだけでしきい値が自動決定される新しい異常検知手法を提案する。 さらに本手法は他の手法で見られるようなデータがなんらかの確率分布に従うという仮定を置かずに利用できる。 これには極値理論を利用する。極値理論は確率分布の最大値(最小値)のふるまいに関する理論である。 中心極限定理と同様に、極値理論では任意の分布の最大値がある分布に分布収束することが知られている。 ただし、一般の分布に対してこのパラメータを推定するのは難しいため Peaks-Over-Threshold (POT) アプローチを用いる。 POT アプローチはあるしきい値を超える値が一般化パレート分布に従うというもので、パラメータ推定はある方程式の根を求めることに帰着する。 複数の根を求める必要があるが、これを最小化問題に帰着し、L-BFGS-B などの一般的で高速な最適化手法で解く。 POT をストリーミングデータに適用するためのアルゴリズム SPOT と基準値が変化(drift)する場合に拡張したアルゴリズム DSPOT を提案する。 これらの手法をいくつかのデータに適用し有用性を示した。 特に汎用PCで1秒間に1000サンプル以上をさばくことができるため、高頻度なストリーミングデータに適用可能である。

参考

www.yasuhisay.info

追記

上の式から下の式の導出が間違っているのでは疑惑がある。

f:id:hoxo_m:20181221141918p:plain

Sparse Gaussian Markov Random Field Mixtures for Anomaly Detection 読んだ

Sparse Gaussian Markov Random Field Mixtures for Anomaly Detection (ICDM 2016) http://ide-research.net/papers/2016_ICDM_Ide.pdf

複数の動作モードに対応可能かつ変数ごとの異常スコアを算出できる新しい異常検知手法を提案する。 通常、異常検知を行いたいシステムは複数の動作モードを持つことが多いが従来の手法ではこれに対応できない。 また、ホテリング T2 などの古典的な手法では多変量であっても異常スコアは観測ごとにしか算出されない。 変数ごとに異常スコアが算出できればシステムのどこに異常が生じたかを突き止めやすくなる。 これらの問題に対応するため、ガウスマルコフ確率場の混合モデルを考え、その変分ベイズ推定アルゴリズムを導出する。 ガウスマルコフ確率場により、ある変数の異常スコアを同じ観測の他の変数の値から求めることができる。 また、混合モデルにより複数の動作モードを表現できる。 動作モードの数(混合数)は不明なため、大きめの数を設定しておけば重みをスパースに推定し、混合数を自動決定する仕組みを取り入れる。 実験として、合成データに対して混合数、混合比率、動作モードについてうまく推定できることを示した。 また、オイル生産コンプレッサーの実データを用いて、他の手法と比較して異常検知性能が良いことを示した。

参考

www.yasuhisay.info

Putting MRFs on a Tensor Train 読んだ

Putting MRFs on a Tensor Train (ICML 2014) http://proceedings.mlr.press/v32/novikov14.pdf

マルコフ確率場(MRF)のパラメータ推定に使われる最尤訓練法では、分配関数の良い近似が必要とされる。 離散変数の正規化されていない同時分布をテンソル(多次元配列)とみなすと、分配関数はその要素の総和として求めることができる。 ただし、テンソルはそのままではメモリに載らないため、テンソル分解の1つであるテンソルトレイン(TT)形式を用いることで容量を削減する。 TT形式の要素の総和を直接的に求めると、TTランクが指数的に増加し、計算時間も増加する。 そこで、TTランクの低い因子(factor)の状態のままで分配関数の近似を計算するアルゴリズムを提案する。 このアルゴリズムにより求められた分配関数について、誤差の上限を理論的に証明し、既存手法より誤差が小さいことを実験により示した。