ポアソン尤度最大化法と最小二乗法による直線の傾き推定

11524 ワード
統計学 統計学 テキストリンク

背景

データがポアソン分布に従っている場合、最小二乗法よりもポアソン尤度最大化法を採用する方がよさそうですが、あまり比較というものをしたことがなかったので、簡単な例で比較してみることにしました。

問題設定

$x$における真の$y$は$y=ax$で与えられるものとします。$y$はポアソン分布に従うものとし、$x$の誤差は考えません。

推定法

まず、最小二乗法を用いてみます。結果だけ示すと、

a = \frac{\sum_ix_iy_i}{\sum_ix_i^2}

となります。
ポアソン尤度最大化法の場合は

a = \frac{\sum_iy_i}{\sum_ix_i}

となります。
素朴に、各点毎に傾きを算出して平均する方法とも比較してみましょう。

a = \sum_i\frac{y_i}{x_i}

コード

#include <random>
#include <iostream>
#include <vector>

double var(const std::vector<double>& a)
{
    int Size = a.size();
    double sum_a = 0;
    double sum_aa = 0;
    for (int i = 0; i < Size; i++)
    {
        sum_a += a[i];
        sum_aa += a[i]*a[i];
    }
    double ret = sum_aa / (double)Size - sum_a * sum_a / (double)Size / (double)Size;
    return ret;
}

int main()
{
    std::mt19937_64 mt(0);
    std::uniform_real_distribution<> u0(1, 5);
    std::uniform_real_distribution<> u1(10, 15);

    int N = 50;//num of data
    double a = 1;

    int Iter = 10000;//試行回数
    std::vector<double> A_LS;
    std::vector<double> A_ML;
    std::vector<double> A_AV;
    for (int iter = 0; iter < Iter; iter++)
    {
        double sum_xy = 0;
        double sum_xx = 0;
        double sum_y = 0;
        double sum_x = 0;
        double av_a = 0;
        for (int n = 0; n < N; n++)
        {
            double x = (n < N / 2 ? u0(mt) : u1(mt));
            std::poisson_distribution<> poi(a * x);
            double y = poi(mt);
            sum_x += x;
            sum_y += y;
            sum_xx += x * x;
            sum_xy += x * y;
            av_a += y / x;
        }
        double a_ls = sum_xy / sum_xx;
        double a_ml = sum_y / sum_x;
        av_a /= (double)N;
        A_LS.push_back(a_ls);
        A_ML.push_back(a_ml);
        A_AV.push_back(av_a);
    //  std::cout << av_a << " " << a_ls << " " << a_ml << std::endl;
    }
    std::cout << sqrt(var(A_AV)) << " " << sqrt(var(A_LS)) << " " << sqrt(var(A_ML)) << std::endl;

    return 0;
}

結果

真値は$a=1$です。データ数は$N=50$としました。
標準偏差は、最も素朴な方法で$0.069$、最小二乗法で$0.054$、ポアソン尤度最大化法で$0.051$となりました。
ちなみに、データが標準偏差$5$の正規分布に従うとしてみると、それぞれ$0.226$、$0.0768$、$0.0911$となりました。

まとめ

傾き推定のような単純な問題でも多少の差は出るようです。より複雑な問題ではもっと違いが出る場合もあるでしょう。
ただし、違いが出ないケースというのもありうるので、その場合はわざわざ面倒なポアソン尤度最大化法を採用する必要はありません。
実際には勾配法等を用いて尤度方程式を解くケースが多く、より単純な最小二乗法の方が安定に解が求まるというケースも考えられます。

C#の静的メソッドと非静的メソッドの比較
C#特性と反射