---
layout: katex
title: ポアソン分布の正規分布による近似
date: 2020-01-04 15:55:45 +0000
tags:
- math
- javascript
---
ポアソン分布の平均mが大きくなると正規分布に近似される様子をグラフで確認するプログラムを作成した。

期待値 m 分散 m の[ポアソン分布](https://ja.wikipedia.org/wiki/%E3%83%9D%E3%82%A2%E3%82%BD%E3%83%B3%E5%88%86%E5%B8%83)の確率分布関数

[[ P(x)=\frac{m^x e^{-m}}{x!} ]]

は、mが大きくなると期待値 m 分散 m の[正規分布](https://ja.wikipedia.org/wiki/%E6%AD%A3%E8%A6%8F%E5%88%86%E5%B8%83)の確率密度関数

[[ f(x)=\frac{1}{\sqrt{2\pi m}} \exp \left( -\frac{(x-m)^2}{2 m} \right) ]]

に近似される。このページでは、その様子をグラフで確認する。

mの値をテキストボックスに直接入力（小数の入力可能）する。ボタンで1ずつ増減できる。

<hr>

m = <input name="m" id="m" type="text" value="3" size="4" onkeyup="update()">
<input type="button" value="-" onclick="decM();">
<input type="button" value="+" onclick="incM();">

<!-- -------------------------------------------------------------------------------------------- -->
<canvas id="canvas" width="600" height="400"
  style="max-width: 100%; height: auto; max-height: 100%">
このブラウザはHTML5のCanvas要素に対応していないためグラフを表示できません。
</canvas>
<script src="/js/graph.js"></script>
<script>
'use strict';
// Maximum of imput parameter
const maxM = 300;
// Figure size
const sizeY = 300;
// Location of the legend
const legendX = 360;
const legendY = 70;
update();

function decM() {
    var m = document.getElementById("m").value;
    m = parseInt(m, 10) - 1;
    if (m < 1 || isNaN(m)) {
        m = 1;
    }
    if (m > maxM) {
        m = maxM;
    }
    document.getElementById("m").value = m;
    update();
}

function incM() {
    var m = document.getElementById("m").value;
    m = parseInt(m, 10) + 1;
    if (isNaN(m)) {
        m = 1;
    }
    if (m > maxM) {
        m = maxM;
    }
    document.getElementById("m").value = m;
    update();
}

function update() {
    // Initialize canvas
    var c, ctx, textM, m, maxX, maxNorm, pZero;
    c = document.getElementById('canvas');
    ctx = c.getContext('2d');
    ctx.clearRect(0, 0, c.width, c.height);
    ctx.font = "20px serif"; // Font of the text
    ctx.lineWidth = 1; // Line width

    // Get parameter
    textM = document.getElementById("m").value;
    m = Number(textM);

    // Check range of m
    if (m > maxM) {
        m = NaN;
        ctx.fillStyle = "red";
        ctx.fillText("このプログラムでは m ≦ " + maxM + " としてください", 80, 300);
    }
    if (m <= 0) {
        m = NaN;
        ctx.fillStyle = "red";
        ctx.fillText("m > 0 でなければなりません", 150, 300);
    }

    // Set Cartesian coodinate system for the graph (GC)
    // Origin of GC with respect to canvas coordinate = (ctx.originX, ctx.originY)
    ctx.originX = 50;
    ctx.originY = sizeY + 70;
    // Unit vector of GC with respect to canvas coordinate = (ctx.unitX, ctx.unitY)
    maxX = m * 2.5;
    if (maxX < 20) {
        maxX = 20;
    }
    ctx.unitX = Math.floor(500 / maxX);
    if (ctx.unitX < 1) {
        ctx.unitX = 1;
    }
    maxNorm = 1 / Math.sqrt(2 * Math.PI * m);
    pZero = Math.pow(Math.E, -m);
    ctx.unitY = -Math.floor(sizeY / Math.max(maxNorm, pZero));

    // Draw graphs
    if (m > 0) {
        ctx.strokeStyle = "black";
        ctx.fillStyle = "black";
        drawAxis(ctx);
        ctx.scaleX = Math.pow(10, Math.floor(2.6 - Math.log10(ctx.unitX)));
        drawScaleX(ctx);
        ctx.scaleY = 1 / Math.pow(10, Math.floor(Math.log10(-ctx.unitY) - 1.5));
        ctx.offsetScaleY = 50;
        drawScaleY(ctx);
        ctx.strokeStyle = "red";
        ctx.fillStyle = "red";
        plotInt(ctx, poisson, m);
        ctx.strokeStyle = "blue";
        ctx.fillStyle = "blue";
        draw(ctx, normDist, m);
        // Legend
        ctx.beginPath();
        ctx.fillStyle = "red";
        ctx.arc(legendX + 15, legendY, 4, 0, Math.PI * 2);
        ctx.fill();
        ctx.fillStyle = "black";
        ctx.fillText("ポアソン分布", legendX + 40, legendY + 5);
        ctx.beginPath();
        ctx.strokeStyle = "blue";
        ctx.moveTo(legendX, legendY + 30);
        ctx.lineTo(legendX + 30, legendY + 30);
        ctx.stroke();
        ctx.fillText("正規分布", legendX + 40, legendY + 35);
    }
}

// Mathematical functions
function poisson(k, m) {
    if (k < 100) {
        return Math.pow(m, k) * Math.pow(Math.E, -m) / factorial(k);
    }
    var logP = k * Math.log(m) - m - logfact(k);
    return Math.pow(Math.E, logP);
}

function normDist(x, m) {
    return Math.pow(Math.E, -(x - m) * (x - m) / (2 * m)) / Math.sqrt(2 * Math.PI *
        m);
}

function factorial(n) {
    if (n < 2) {
        return 1;
    } else {
        return n * factorial(n - 1);
    }
}

function logfact(n) {
    if (n < 2) {
        return 0;
    } else {
        return Math.log(n) + logfact(n - 1);
    }
}
</script>
<!-- -------------------------------------------------------------------------------------------- -->