package info.tsujimotter.math;
/**
* 計算機統計計算のためのライブラリクラスです。
*
* @version 1.0, 27 July, 2012
* @author Junpei Tsuji
*/
public class Statistics {
/**
* log sum exponentialを計算する関数です。
*
* W(x, y) = log (exp(x) + exp(y))
*
* のような式を計算する際に、x, yがある程度大きい（あるいは小さい）場合、exp(x), exp(y)の値がオーバーフロー（アンダーフロー）してしまいます。
* この場合に、次のように変換すると、オーバーフロー（またはアンダーフロー）を防ぐことが出来ます。
*
* W(x, y) = log exp(x) (1.0 + exp(y-x))
* ∴W(x, y) = x + log (1.0 + exp(y-x))
*
* （ただし、x > y を仮定します。）
*
* これにより、指数部の計算がexp(y-x)に収まります。
*
* logsumexp(x, y)関数は、上記の要領で、W(x, y)を計算する関数です。
*
* 下式のような三項間の和を求める際も、二項間の和 W(x, y)を用いて次のように計算可能です。
*
* W(x, y, z) = log (exp(x) + exp(y) + exp(z))
*
* W(x, y, z) = W(W(x, y), z)
* = W(log(exp(x) + exp(y)), z)
* = log(exp(log(exp(x) + exp(y))) + exp(z))
* = log(exp(x) + exp(y) + exp(z))
*
* n 個の指数和を for 文などで求める場合も同様ですが、最初の和だけ注意する必要があります。
*
* double[] x = {x1, x2, x3};
* double w;
* for(int i=0; i< n; i++)
* w = logsumexp(w, x[i], i==0);
*
* i==0 のときの和は通常、 w が初期化されていないため不定になってしまいます。
* logsumexp関数では、この問題を防ぐため、flg引数を用意し、初期状態のときのみ y をそのまま返すようにしています。
*
* @param x 引数1
* @param y 引数2
* @param flg 初期状態を示すフラグ変数。trueのとき、logsumexp関数は 第二引数 y を返す。
* @return logsumexp(x, y)
*/
public static double logsumexp (double x, double y, boolean flg)
{
if (flg) return y; // init mode
if (x == y) return x + 0.69314718055; // log(2)
double vmin = (x<=y) ? x : y;
double vmax = (x>=y) ? x : y;
if (vmax > vmin + 50) {
return vmax;
} else {
return vmax + Math.log (Math.exp (vmin - vmax) + 1.0);
}
}
}