2023年8月11日金曜日

不確実性を減少させる意思決定

今回は環境を取り巻くランダムネスの影響を小さくする意思決定の方法について考えていきます。
ノイズの標準偏差が説明変数を全て含む多変数の一次関数になる場合の不等分散重回帰モデルにおいて、不確実性を減少させるための説明変数の状態について、予算配分問題を例に考えていきます。

予算配分問題

多角的に事業を展開している企業において、事業の成長見込みを基に来年度予算割り当て額を決める問題を考えます。
$y$={事業1の利益, 事業2の利益, ...}、$x$={事業1の割り当て予算, 事業2の割り当て予算,...}、$ \alpha$={事業1の投下資本利益率(ROIC)の予測値,事業2のROICの予測値,...}
とした時、実際のROICは予測値にノイズが加わるため
$$ y_i=x_i( \alpha_i + \varepsilon_i)=\alpha_i x_i + x_i \varepsilon_i, \quad \varepsilon_i \sim N(0, \sigma_i^2), \quad \alpha_i \geq 0 $$
という関係式を考えることが出来ます。
この時の$x_i$を与えた上での$y_i$の条件付き分布は、$y_i|x_i \sim N(\alpha_i x_i, (x_i \sigma_i)^2)$になります。

事業数を$D$として、会社全体の利益は$t=\sum_i^D y_i$となります。 

この節では簡単のため、$D=4$とし、各事業の利益は他事業に依存せず独立に決まるとします(独立性の仮定)。 

$x$を与えた上での$t$の条件付き分散は確率変数の和の公式によって、$V[t|x]=\sum_{i=1}^4 (x_i \sigma_i)^2$となります。
$V[t|x]$が小さくなることはランダムネスの影響を小さくし、予定通りの結果を出しやすくなる確率が上昇し、計画性の向上や信頼性の保持に繋がりますので、これを最小化する$x$を求めます。
$\sum_{i=1}^4 x_i=a$の条件のもとで$V[t|x]$を最小化する$x$はラグランジュ未定乗数法によって
$$ x_i=\frac{\Pi_{j\neq i}^4 \sigma_j^2}{\sum_{j=1}^2\sum_{k=j+1}^3\sum_{l=k+1}^4 (\sigma_j\sigma_k\sigma_l)^2}a $$ となります。 上式は分子の$\sigma_i^2$以外の分散の積が、分母の${}_4 C_{4-1}=4$の組み合わせ全体の中での比率を求める計算に$a$を掛けたものになっているため、$\sigma_1=...=\sigma_4$の時、ランダムネスの影響が最小化される最適点は$x_1=...=x_4=a/4$となります。

従ってこのような単純化された状況ではバランスが取れた予算配分には計画性を向上させるメリットがあることが分かりました。
一方、予算配分をバランスよくさせることは、最もROICが高い事業に全予算を割り当てる場合と比べて利益の期待値である$E[t|x]$が下がるデメリットがあるため、$E[t|x]$と$V[t|x]$の兼ね合いを考える必要があります。
また、事業間シナジーを考慮し、独立性の仮定も外す必要があります。次節でこれらを考慮した最適化問題を考えていきます。

評価指標の導出

事業$i$の予算額$x_i$に応じてもたらされる、観測されない利益$y_i$の$x_i$を与えた上での条件付き分布に正規性の仮定をおきます。 $$y_i|x_i \sim N\left(\alpha_i x_i, (x_i\sigma_i)^2\right), \quad \alpha_i \geq 0 $$ 会社全体の利益は $$ t = \sum_{i=1}^D y_i$$ とします。先ほどとは違い独立性の仮定を置きません。

$t$の条件付き分布は、正規分布の再生性によって、次式になることが分かります。

$$t|x \sim N\left( \sum_{i=1}^D \alpha_i x_i, \sum_{i=1}^D V(y_i|x_i) + 2\sum_{i=1}^D\sum_{j=i+1}^D \Lambda_{ij} \right)\quad (1)$$

ここで$\Lambda$は分散共分散行列で $$ \Lambda_{ij}= \rho_{ij. x}S(y_i|x_i)S(y_j|x_j), \quad 1 \geq \rho_{ij. x} \geq -1 $$と定義されます。$ \rho_{ij . x}$は$x_i$と$x_j$の影響を取り除いた$y_i$と$y_j$の偏相関で、$\varepsilon_i$と$\varepsilon_j$の相関係数とみなすことが出来ます。予算額の類似性による影響を除いても利益に相関関係がみられる場合にそれを代入するものになります。$S(y_i|x_i)$は条件付き標準偏差で、$S(y_i|x_i)=x_i\sigma_i$です。

目的はこの正規分布の期待値の高さと、計画通りの結果を出しやすくする度合いを表す標準偏差の小ささ、この二つの情報が含まれていて最適化の目的関数になり得る指標を導出する事ですが、$t$が期待値$\mu:=\sum_{i=1}^D \alpha_i x_i$以下の値を取ったときの条件付き期待値、すなわち、失敗したときの期待値がそれに該当します。

これから行う条件付き期待値の導出は、「正規分布全体の期待値の導出」と共通する部分が多いため、そちらを見て分かるところは省略します。

\begin{align} E[t| t \leq \mu] &= \frac{E[t, t \leq \mu]}{P(t \leq \mu)} \\ &= 2 \int_{-\infty}^\mu t N(t, \mu, \sigma^2) dt \\ &= \frac{2}{ \sqrt{2 \pi \sigma^2 } } \int_{-\infty}^\mu (t-\mu) e^{-\frac{(t-\mu)^2}{2\sigma^2}} dt + 2\mu \int_{-\infty}^\mu \frac{1}{ \sqrt{2 \pi \sigma^2 } } e^{-\frac{(t-\mu)^2}{2\sigma^2}} dt \\ &= \frac{2}{ \sqrt{2 \pi \sigma^2 } } \int_{-\infty}^\mu (t-\mu) e^{-\frac{(t-\mu)^2}{2\sigma^2}} dt + \mu \end{align} 第1項目について、$z=\frac{t-\mu}{\sigma}$とおくと、$t=\sigma z+\mu$となり、単調関数で表すことが出来るため、定積分の置換積分法を適用すると、 \begin{align} \frac{2}{ \sqrt{2 \pi \sigma^2 } } \int_{-\infty}^\mu (t-\mu) e^{-\frac{(t-\mu)^2}{2\sigma^2}} dt &= \frac{2}{ \sqrt{2 \pi \sigma^2 } } \int_{-\infty}^\frac{\mu-\mu}{\sigma} \sigma z e^{-\frac{z^2}{2}} \sigma dz \\ &= \frac{2 \sigma}{ \sqrt{2 \pi } } \int_{-\infty}^0 z e^{-\frac{z^2}{2}} dz \\ &= -\frac{2 \sigma}{ \sqrt{2 \pi } } \end{align} 最後の等式はガウス積分の公式を使うことで分かります。最後に第1項と第2項を合わせた上で、式(1)の結果を代入すると指標の完成です。 $$  f(x) = \sum_{i=1}^D \alpha_i x_i -\frac{2}{\sqrt{2 \pi }} \sqrt{\sum_{i=1}^D (\sigma_i x_i)^2 + 2\sum_{i=1}^D\sum_{j=i+1}^D \Lambda_{ij}}\quad (2) $$

$\sum_i^D x_i \leq a $の条件の下で$f(x)$が増加するよう予算配分を調整することを考えることが出来ます。しかし、期待値以下の実績を失敗と定義するのはかなり保守的ですので、必要に応じて第2項の重み係数$2/\sqrt{2\pi}$を下げることになると思います。

2023年8月8日火曜日

結合エントロピーと汎化能力

今回は次回の記事で重要になる結合エントロピーについて考えます。

結合エントロピーの展開

結合エントロピーは次式によって求めることが出来ます。

$$ H(x_1,x_2,...,x_D) = \sum_{i=1}^D H(x_i) - \sum_{i=1}^D \sum_{j=i+1}^D I(x_i, x_{j}|x_{j+1},x_{j+2},...,x_D)\quad (1) $$ ここで$I(A,B|C)=H(A|C)+H(B|C)-H(A,B|C)$は条件付き相互情報量です。
証明には次の公式を使います。
・公式1 $H(A,B)=H(A|B)+H(B)\Leftrightarrow H(A|B)=H(A,B)-H(B)$
・公式2 $H(A,B)=H(A)+H(B) - I(A,B)$
チェーンルール $$ H(x_1,x_2,...,x_D)=H(x_1) + \sum_{i=2}^D H(x_i|x_1,...,x_{i-1}) $$ チェーンルールは次のように書き換えることが出来ます。
・公式3 $$ H(x_D,x_{D-1},...,x_1)=H(x_D) + \sum_{i=1}^{D-1} H(x_{D-i}|x_{D-i+1},...,x_D) $$ 先に次の等式を証明します。 $$ H(x_i|x_{i+1},..,x_D) = H(x_i) - \sum_{j=1}^D I(x_i, x_{i+j}|x_{i+j+1},x_{i+j+2},...,x_D)\quad (2) $$ 証明 $$ \begin{align} f_i(j) := H(x_i|x_{i+j},..,x_D) &= H(x_i,x_{i+j},..,x_D) - H(x_{i+j},..,x_D) \\ &= H(x_i,x_{i+j},...,x_{D-1}|x_D) + H(x_D) - H(x_{i+j},..,x_D) \\ &= H(x_i,x_{i+j},...,x_{D-2}|x_{D-1},x_D) + H(x_{D-1}|x_D) + H(x_D) - H(x_{i+j},..,x_D) \end{align} $$ 第1項に立て続けに公式1を適用すると $$ \begin{align} & H(x_i,x_{i+j},...,x_{D-2}|x_{D-1},x_D) + H(x_{D-1}|x_D) + H(x_D) - H(x_{i+j},..,x_D) & \\ & = H(x_i,x_{i+j}|x_{i+j+1},..,x_D) + H(x_D) + \sum_{k=i+j+1}^{D-1} H(x_{D-k}|x_{D-k+1},...,x_D) - H(x_{i+j},..,x_D) & \\ \end{align} $$ 公式3を適用し $$ \begin{align} & = H(x_i,x_{i+j}|x_{i+j+1},..,x_D) + H(x_{i+j+1},...,x_D) - H(x_{i+j},..,x_D)\quad (3) & \end{align} $$ 第1項について、公式2より $$ H(x_i,x_{i+j}|x_{i+j+1},..,x_D) = H(x_i|x_{i+j+1},..,x_D) + H(x_{i+j}|x_{i+j+1},..,x_D) - I(x_i, x_{i+j}|x_{i+j+1},..,x_D) $$ 上式の第2項について、 $$ H(x_{i+j}|x_{i+j+1},..,x_D) = H(x_{i+j},..,x_D) - H(x_{i+j+1},..,x_D) $$ 従って式(3)は、 $$ \begin{align} & H(x_i,x_{i+j}|x_{i+j+1},..,x_D) + H(x_{i+j+1},...,x_D) - H(x_{i+j},..,x_D) & \\ & = H(x_i|x_{i+j+1},..,x_D) + H(x_{i+j},..,x_D) - H(x_{i+j+1},..,x_D) - I(x_i, x_{i+j}|x_{i+j+1},..,x_D) + H(x_{i+j+1},...,x_D) - H(x_{i+j},..,x_D) & \\ & = H(x_i|x_{i+j+1},..,x_D) - I(x_i, x_{i+j}|x_{i+j+1},..,x_D) \end{align} $$ すなわち $$ f_i(j) = f_i(j+1) - I(x_i, x_{i+j}|x_{i+j+1},..,x_D) $$ となり、条件となる変数を一つづつ減らしながら再帰的に展開していくことが出来、$j=D$については、 $$ \begin{align} f_i(D) &= H(x_i|x_D) = H(x_i,x_D)-H(x_D)\\ &= H(x_i)+ H(x_D)-I(x_i,x_D)-H(x_D) \\ &= H(x_i)-I(x_i,x_D) \end{align} $$ となるため式(2)は成り立ちます。 

 式(1)についても公式3と式(2)を適用することで証明することが出来ます。

$x_i$に正規性の仮定をおくと、$H(A)$は分散に帰着され[1]、$I(A,B|C)$は偏相関に帰着されます[2]。偏相関の計算は回帰分析の残差を使う方法が提案されています[3]。

多変量正規分布のエントロピーの計算には式(1)を使う方法以外にも分散共分散行列の行列式を使って求める方法[4]があります。

結合エントロピーと分散の和の関係

独立な多変量正規分布$P(X_1, X_2,...,X_D)=\Pi_i^D N(\mu_i, \sigma_i^2)$に従う確率変数の結合エントロピーは分散の合計値だけで決まるわけではなく、それぞれの分散のバランス状態によっても変わります。

多変数の独立性から

$$ H(X_1,X_2,...,X_D)=\sum_i^D H(X_i)=\frac{D}{2}\{1+log(2\pi)\}+\frac{1}{2}log(\Pi_i^D \sigma_i^2)$$

ここで$\sum_i^D \sigma_i^2=a$の制約の下で$H(X_1,...,X_D)$が最大化する$\sigma^2$を考えたとき、最適化に影響しない定数及び単調関数を除くと目的関数は$\Pi_i^D \sigma_i^2$となり、相加相乗平均の不等式と累乗根の単調性から、$\sigma_1^2=\sigma_2^2=...=\sigma_D^2$となるときのみ結合エントロピーが最大化することが分ります。

説明変数の結合エントロピー

回帰分析において、平均が0になるよう中心化した説明変数の結合エントロピーはその説明変数の分散共分散行列との関連性が認められるため、そのエントロピーの増大が意味するところは、多重共線性の問題が発生しにくい望ましい状態になる確率が高まることと言えます。

ここで次の線形重回帰モデルを考えます。 $$ y_i = w_1x_{i1} + w_2x_{i2} + \varepsilon_i,\quad \varepsilon_i \sim N(0, \sigma^2),\quad i=1,2,...,n$$ $x$と$y$は中心化されているとします。$w$の最尤推定量は $$ w=(x^Tx)^{-1}x^Ty $$ となります。$x$が中心化されているため、$(x^Tx)$は分散共分散行列を$n$倍したものになります。$(x^Tx)_{ij}=n a_{ij}$とし、上式に基づいて$w_1$を求める式を計算すると、次式になります。 $$ \begin{align} w_1 &= \sum_{i=1}^n y_i \frac{n(x_{i1}a_{22} - x_{i2}a_{12})}{n^2(a_{11}a_{22} - a_{12}a_{21})} \\ &= \frac{1}{n}\sum_{i=1}^n y_i \frac{x_{i1}a_{22} - x_{i2}a_{12}}{a_{11}a_{22} - a_{12}^2} \quad (4) \end{align} $$ $a_{12}$が大きくなれば$(x^T)_1$と$y$との関連性の高さに関係なく$w_1$が大きくなってしまい本来の意味を持ちません。一方$a_{12}=0$とした場合、上式は$x, y$が中心化されていることに留意して $$ \frac{1}{n}\sum_{i=1}^n y_i \frac{x_{i1}a_{22} - x_{i2}a_{12}}{a_{11}a_{22} - a_{12}^2} = \frac{1}{n}\sum_{i=1}^n \frac{y_i x_{i1}}{a_{11}} = \rho_{y,1}\frac{\sigma_y}{\sigma_{x,1}}   $$ となり、本来の意味を持ちます。$\rho_{y,j}$は$y$と$(x^T)_j$との相関係数で、$\sigma_{x,j}$、$\sigma_y$はそれぞれ$(x^T)_j$及び$y$の標準偏差です。

また、$a_{22}$を大きくしても相対的に$a_{12}$が回帰分析に与える影響を小さくすることが出来、本来の意味を持ちます。加えて、式(4)を分散$a_{jj}$の関数として見たとき、分母は2次式、分子は1次式となるため、分散の合計が回帰係数に影響を与えます。さらに言うと$a_{11}+a_{22}=b$の制約下で分母の$a_{11}a_{22}$が最大の値を取るのは前節で解説した通り、結合エントロピーと同じく、$a_{11}=a_{22}$となる時であることが分かります。以上をまとめると、説明変数同士の相関の小ささ、分散の合計、分散のバランス性の3つが回帰係数に影響を与えるため、式(1)及び前節での結論を踏まえると説明変数の結合エントロピーが増加すると、解がスパースになる傾向がうかがえます。このことから、汎化能力の観点から、説明変数のエントロピーの重要性を認識することが出来ます。

参考文献

  1. 正規分布のエントロピーを丁寧に計算する
  2. 正規分布に従う連続確率変数の相互情報量の推定は相関係数の推定に帰着される
  3. Partial correlation - Wikipedia
  4. Entropy of the multivariate Gaussian

2014年2月24日月曜日

【C言語】ベイジアンネットワークの学習アルゴリズム

今回はベイジアンネットワークの学習アルゴリズムを中心に解説を行います。
■確率的因果
事象A、Bがあり、BがAの原因であるとは、Bが発生したとき、Aが発生する確率を高めることをいいます。

例)
1000人の体型を調査したところ肥満と判断される確率は0.3となりました。ただし体重が80Kg以上の人に限っては肥満であると判断される確率は0.6に高まったため、
体重80Kg以上は肥満の原因であるといいます。

事象A、B、C、...があり、いかなる事象を同時にとってもBが発生したとき、Aが発生する確率を高める場合、BはAの真の原因であるといいます。

例)
体重が80Kg以上の人が肥満と判断される確率は0.6ですが、体重が80Kg以上かつ身長が183cm以上の人が肥満であると判断される確率は元の0.3に戻ったため
体重80Kg以上は肥満の真の原因ではありません。

■ベイジアンネットワーク
ベイジアンネットワークは確率的因果関係をグラフ理論を用いてモデリングする技術です。
ベイジアンネットワークはノードとアーク(矢印)、条件付き確率表から構成されます。ノードは確率変数、各ノード間をつなぐアークは条件付き確率に対応します。
また、ベイジアンネットワークはノードをアークの方向にたどっていって同じノードに戻らないよう構築します。
■ベイジアンネットワークの学習
データからベイジアンネットワークを構築する方法をベイジアンネットワークの構造学習のアルゴリズムといいます。
構造学習のアルゴリズムには大きく分けて、厳密解と、精度を落として計算速度を高めた近似解がありますが、最近の研究のトレンドは厳密解とのことなので、厳密解について説明します。
もっとも基本的な学習アルゴリズムの手順は、考えられるすべての構造パターンのベイジアンネットワークを構築し、個々のベイジアンネットワークについてBICなどの情報量基準を学習スコアとして求め、学習スコアが最小または最大となる構造を持つベイジアンネットワークを真の構造と仮定し採択します。

[最適化する前の厳密解構造学習]
ただし、ベイジアンネットワークの学習アルゴリズムの計算量はノードの数に応じて指数倍的にふえていくため一般的なデスクトップPCではノード数が6個か7個以上の計算は難しくなってくると思います。
この問題に対してSilander and Myllymaki(2006)は動的計画法を使った最適化を提案しています。アルゴリズムの詳細は記事末で紹介している参考書に詳しく書かれているのでここではソースコードの提供のみを行います。動的計画法で最適化されたBIC基準の厳密解学習アルゴリズムと列挙法による確率推論のソースコードを下記のgithubに上げています。
https://github.com/y-mitsui/DpBayesianNet

■検証
Kaggleのタイタニックコンペで提供されているデータを使ってベイジアンネットワークを検証しました。
設定した確率変数は次のとおりです。
性別。年齢。グレード。同一のチケット番号の乗客が全員死亡:2、全員生存:0、それ以外:1。部屋の位置が左か右か。家族と同乗か否か。
予測精度は77%となりました。
全ソースはgithubに上げています。
https://github.com/y-mitsui/kaggleTitanic

参考書籍
ベイジアンネットワーク 植野真臣

2014年1月15日水曜日

【C言語】高速ガウス変換を使ってカーネル密度推定を高速化

高速ガウス変換の紹介とそれを用いてカーネル密度推定を高速化します。
 x1, x2, ..., xN を独立かつ同一な分布に従う標本としたときの予測モデル、

このモデルの類似モデルとして、ガウシアンカーネルを使った、「1次元カーネル
密度推定」、「ミーンシフト」、「カーネル回帰」などが挙げられますが、このモデルは多項式のNが標本数と同じになるため、予測回数をMとするとき、その計算量は非常に高価で、O(MN)となります。
この問題に対して高速多重極法を起源とする高速ガウス変換では、エルミート展開を用いてガウシアンの近似式を作成することで、上式の多項式のNを誤差率に応じて決定される小さな非負整数Kに書き換えられ、その計算量はO(MK)となります。

この高速ガウス変換を改良したIFGTがありますが、
今回はオリジナルの高速ガウス変換を使ったカーネル密度推定のCライブラリを開発し、githubに上げました。
http://github.com/y-mitsui/fastKDE
 また、32bit Linux環境ではSSE2で最適化されたアセンブラコードが利用できます。
詳細はREADME.mdをご覧ください。

また、前述のとおり高速ガウス変換には Kに応じた誤差が伴いますので、厳密な計算が必要な場面で利用する際には注意が必要です。

■実行結果
1000000 個の標準正規乱数からなる標本、予測回数を5000回 とし、Kを10と設定した場合の
1次元カーネル密度推定のCPU時間は最適化されていない場合の444秒に対して、
高速ガウス変換を使用した場合は0.6秒となり、大幅に速度が向上していることがわかります。
一方、Kの値が小さすぎた為、誤差が大きくなっています。

direct time:444.13400sec optimized time:0.60800sec
directValue[00000]:0.221897     optimizedValue[00000]:0.242221
directValue[00250]:0.233130     optimizedValue[00250]:0.252617
directValue[00500]:0.243662     optimizedValue[00500]:0.262466
directValue[00750]:0.253349     optimizedValue[00750]:0.271657
directValue[01000]:0.262056     optimizedValue[01000]:0.280076
directValue[01250]:0.269658     optimizedValue[01250]:0.287611
directValue[01500]:0.276041     optimizedValue[01500]:0.294152
directValue[01750]:0.281110     optimizedValue[01750]:0.299593
directValue[02000]:0.284789     optimizedValue[02000]:0.303839
directValue[02250]:0.287019     optimizedValue[02250]:0.306803
directValue[02500]:0.287767     optimizedValue[02500]:0.308415
directValue[02750]:0.287020     optimizedValue[02750]:0.308619
directValue[03000]:0.284790     optimizedValue[03000]:0.307384
directValue[03250]:0.281112     optimizedValue[03250]:0.304698
directValue[03500]:0.276040     optimizedValue[03500]:0.300573
directValue[03750]:0.269653     optimizedValue[03750]:0.295048
directValue[04000]:0.262046     optimizedValue[04000]:0.288185
directValue[04250]:0.253331     optimizedValue[04250]:0.280067
directValue[04500]:0.243634     optimizedValue[04500]:0.270803
directValue[04750]:0.233090     optimizedValue[04750]:0.260516

CPU: Core i5(Nehalem)  MEM:4GB OS:Windows 7 COMPILER:mingw 4.62

2013年11月20日水曜日

【GAS】アセンブラで名前付きローカル変数を利用できるようにするm4マクロ

アセンブラ開発でしばしば使われるm4マクロを使って
高級言語のように名前付き変数を実現する方法を紹介します。以下のコードはLinux GCC専用です。

define(`tmpcount',0)
define(`setVars',`ifelse($#,1,,eval($#%2),0,`dnl
        define(`count$1',eval(tmpcount+$2)) dnl
        define(`tmpcount',eval(tmpcount+$2)) dnl
        define($1,`-count$1 (%ebp)') dnl
        setVars(shift($@))',`setVars(shift($@))') dnl
')

上記のマクロを定義して、以下のように、setvarsマクロによって変数の宣言を行います。

hello: .string "hello!"
.global main
main:
        pushl %ebp
        movl  %esp,%ebp
        subl  $8,%esp
        setVars(var1,4,var2,4)
        movl    $1,var1
        movl    $2,var2
        pushl $hello
        call puts
        leave
        ret

setVarsマクロはsetVars(変数1の名前,変数1のバイト数,変数2の名前,変数2のバイト数,...)
といった具合に宣言を行います。

2013年10月19日土曜日

【C言語】カーネル回帰分析の実装と検証

非線形な関係を持つデータ間をカーネル法をつかって回帰分析してみます。 また、後半で本手法の問題点も示します。
このようなデータを当てはめてみることにします。
このようになります。高い表現力があることがわかります。

コードはこちら。行列計算にはmeschというライブラリを使いました。
// kernelRegression
// x:data
// dimention: x's dimention
// y: label of x 
// num: number of data
// result: coeffience
double *kernelRegression(double **x,int dimention,double *y,int num,double gamma,double regilarization){
 VEC *label;
 MAT *gram,*ident,*inv;
 int i,j;
 double *result;

 ident = m_get(num,num);
 label = v_get(num);
 memcpy(label->ve,y,sizeof(double)*num);
 gram = m_get(num,num);
 ident = m_get(num,num);
 for(i=0;i<num;i++){
  for(j=0;j<num;j++){
   gram->me[i][j]=kernel(x[i],x[j],dimention,gamma);
  }
 }
 inv=m_add(gram,sm_mlt(regilarization,m_ident(ident),NULL),NULL);
 inv=m_inverse(inv,NULL); 
 VEC *coefficient=v_get(num);
 mv_mlt(inv,label,coefficient);
 result=malloc(sizeof(double)*num);
 memcpy(result,coefficient->ve,sizeof(double)*num);

 m_free(ident);
 m_free(gram);
 m_free(inv);
 v_free(label);
 v_free(coefficient);
 return result;
}}
カーネルがこちら。
double kernel(double *x1,double *x2,int num,double gamma){
 VEC *xVector,*yVector,*diffVector;
 int i;

 xVector=v_get(num);
 memcpy(xVector->ve,x1,sizeof(double)*num);
 yVector=v_get(num);
 memcpy(yVector->ve,x2,sizeof(double)*num);

 diffVector=v_sub(xVector,yVector,NULL);

 double sum=0.0;
 for(i=0;i<num;i++){
  sum+=diffVector->ve[i]*diffVector->ve[i];
 }

 v_free(diffVector);
 v_free(xVector);
 v_free(yVector);
 
 return exp(-gamma*sum);
}
全ソースはgithubに上げています。
https://github.com/y-mitsui/kernelRegression
このように高い表現力を持つカーネル回帰ですが、多項式の次数がサンプル数と同じなのでサンプル数が増えると計算量が増えます。
また、線形の回帰分析と同様、誤差関数に二乗誤差を用いているので外れ値の影響を受けやすく、外れ値に対して積極的に当てはめようと するので、他のデータの当てはまりが悪ることがあります。
これらの問題を抑制させる回帰分析がサポートベクトル回帰で、そちらのほうがよく使われるように思います。

2013年9月10日火曜日

【C言語】ミーンシフトやKDEのバンド幅を自動調整

カーネル密度推定(KDE)はデータからノンパラメトリック(特定の確率分布を想定しない)に確率密度を推定する方法です。
また、ミーンシフトはKDEの理論を利用して関数の極値を求める方法で、カーネル数を事前に設定しなくてもよいクラスタリングとしてよく利用されていると思います。
しかし、どちらもバンド幅を設定する必要があります。
そのバンド幅をデータから自動調整する方法の一つにPlugin法があります。
Plugin法はデータの確率分布が正規分布に従うことを想定してます。
double variance(double *data,int num){
 int i;
 double mean=0,sum=0;
 
 for(i=0;i<num;i++){
  mean+=data[i];
 }
 mean/=num;

 for(i=0;i<num;i++){
  sum+=pow(mean-data[i],2);
 }

 return sum/num;
}
//プラグイン法によるバンド幅
//data:データ配列
//num:配列の要素数
//戻り値:バンド幅
double calcBandWidth(double *data,int num){
 return 1.06*sqrt(variance(data,num))*pow(num,-0.2);
}