ポアソン分布
事故の確率や、不良品の確率の推定に広く利用されている計算です。
応用例については、インターネットで "ポアソン分布" で検索して下さい。
Γ(k+1,λ)は、不完全ガンマ関数です。
λの値が大きくなると、正規分布に近づきます。
ｋの値は、１，２，３・・・ の整数です。
unit Main;
interface
uses
Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.StdCtrls, Vcl.ExtCtrls,
VclTee.TeeGDIPlus, Vcl.Buttons, VCLTee.TeEngine, VCLTee.Series,
VCLTee.TeeProcs, VCLTee.Chart, system.Math, system.UITypes;
type
TForm1 = class(TForm)
LabeledEdit1: TLabeledEdit;
LabeledEdit2: TLabeledEdit;
LabeledEdit3: TLabeledEdit;
RadioGroup1: TRadioGroup;
Chart1: TChart;
Series1: TBarSeries;
BitBtn1: TBitBtn;
Memo1: TMemo;
Series2: TPointSeries;
procedure FormCreate(Sender: TObject);
procedure BitBtn1Click(Sender: TObject);
private
{ Private 宣言 }
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
var
Epsilon: double;
const
// doubleで計算できるガンマ関数の最大値
// exp の計算で制限されます
max_double_arg: double = 169;
//ベルヌーイ数B(2)、B(4)、B(6)、...、B(16)
B : array[1..8] of double = ( 1.0 / 6, // B(2)
-1.0 / 30,
1.0 / 42,
-1.0 / 30,
5.0 / 66,
-691.0 / 2730,
7.0 / 6,
-3617.0 / 510); // B(16)
m: integer = sizeof(B) div sizeof(double);
// ガンマ関数
function Gamma(x: double): double;
var
i : integer;
sum : double;
xx, v : double;
xj : double;
lng : double;
ln_sqrt_2pi: double;
begin
ln_sqrt_2pi := ln(sqrt(2 * pi)); // = ln(2 * pi) / 2;
sum := 0;
v := 1;
while x < m do begin
v := v * x;
x := x + 1;
end;
xx := x * x;
xj := x;
lng := ln_sqrt_2pi - x + (x - 0.5) * ln(x) - ln(v);
for i := 1 to m do begin
sum := sum + B[i] / (i * 2 * (2 * i - 1)) / xj;
xj := xj * xx;
end;
result := exp(sum + lng);
end;
// γ(a, x) 第１種不完全ガンマ関数
function incomplete_gamma1(a, x: double): double;
var
k : integer;
term, previous : double;
begin
result := 0;
if x = 0 then exit;
term := exp(a * ln(x) - x) / a;
result := term;
for k := 1 to 1000 do begin
term := term * x / (a + k);
previous := result;
result := result + term;
if result = previous then begin
exit;
end;
end;
Form1.Memo1.Lines.add ('p_gamma(): 収束しません');
end;
// Γ(a, x) 第２種不完全ガンマ関数
function incomplete_gamma2(a, x: double): double;
var
k : integer;
w, temp, previous : double;
la, lb : double;
begin
if x < 1 + a then begin
result := Gamma(a) - incomplete_gamma1(a, x);
exit;
end;
la := 1; // Laguerreの多項式
lb := 1 + x - a;
w := exp(a * ln(x) - x);
result := w / lb;
for k := 2 to 1000 do begin
temp := ((k - 1 - a) * (lb - la) + (k + x) * lb) / k;
la := lb;
lb := temp;
w := w * (k - 1 - a) / k;
temp := w / (la * lb);
previous := result;
result := result + temp;
if result = previous then exit;
end;
Form1.Memo1.Lines.add ('q_gamma: 収束しません');
end;
// 確率質量の計算
function probability_mass(l: double; k:double): double;
begin
result := power(l, k) * exp(-l) / Gamma(k + 1);
end;
// 計算の実行
procedure TForm1.BitBtn1Click(Sender: TObject);
var
i, k: integer;
l, m, p, x: double;
min, max: integer;
pl, ppl : double;
Q, s : double;
dp : double;
begin
val(LabeledEdit1.Text, l, i);
if i <> 0 then begin
MessageDlg('λの値に間違いがあります。', mtInformation,
[mbOk], 0, mbOk);
exit;
end;
val(LabeledEdit2.Text, k, i);
if i <> 0 then begin
MessageDlg('kの値に間違いがあります。', mtInformation,
[mbOk], 0, mbOk);
exit;
end;
val(LabeledEdit3.Text, p, i);
if i <> 0 then begin
MessageDlg('Pの値に間違いがあります。', mtInformation,
[mbOk], 0, mbOk);
exit;
end;
if p >= 1 then begin
MessageDlg('Pの値が大きすぎます。', mtInformation,
[mbOk], 0, mbOk);
exit;
end;
min := 0;
max := round(l) * 2;
if max < 13 then max := 13;
if k > trunc(max_double_arg) then begin
MessageDlg('kの値が大きすぎて計算できません。', mtInformation,
[mbOk], 0, mbOk);
exit;
end;
if (k > max) then begin
if k <= max_double_arg then max := trunc(max_double_arg);
if k < max_double_arg - 1 then max := k + 1;
end;
Series1.Clear;
Series2.Clear;
memo1.Clear;
// 確率質量 累積確率 の計算
pl := 0;
x := 0;
s := 0;
for i := min to max do begin
m := probability_mass(l, i); // 確率質量
if RadioGroup1.ItemIndex = 0 then begin
if i < k then Series1.AddXY(i, m, '', clRed);
if i = k then Series1.AddXY(i, m, '', clFuchsia);
if i > k then Series1.AddXY(i, m, '', clBlue);
end;
if i <= k then pl := pl + m; // 下側累積確率
if i = k then x := m;
if RadioGroup1.ItemIndex = 1 then begin
s := s + m; // 累積確率
Series2.AddXY(i, s);
end;
end;
Q := 1 - pl + x; // 上側累積確率
memo1.Lines.Add('パーセント点(k=' + intTostr(k) + ')');
memo1.Lines.Add(' 確率質量 f ' + floatTostr(x));
memo1.Lines.Add('下側累積確率 P ' + floatTostr(pl));
memo1.Lines.Add('上側累積確率 Q ' + floatTostr(Q));
// 下側累積確率からパーセント点の計算
// 逆関数が難しいので、近似法により計算
x := trunc(max_double_arg);
for i := 0 to trunc(max_double_arg) do begin
ppl := incomplete_gamma2(i + 1, l) / Gamma(i + 1); // 下側累積確率
if ppl >= p then begin
x := i;
break;
end;
end;
dp := 1;
repeat
x := x + dp;
ppl := incomplete_gamma2(x, l) / Gamma(x); // 下側累積確率
s := ppl - p;
if s > 0 then begin
x := x - dp;
dp := dp / 2;
end;
until (abs(s) <= Epsilon) or (dp <= Epsilon);
x := x - 1;
memo1.Lines.Add('下側累積確率(P = ' + floatTostr(p) + ')');
// 改行なし出力
memo1.Lines.Text := memo1.Lines.Text
+ ' パーセント点 k ' + floatTostr(x);
// memo1.Lines.Add(' パーセント点 k ' + floatTostr(x));
k := trunc(x);
if RadioGroup1.ItemIndex = 2 then begin
for i := 0 to max do begin
m := probability_mass(l, i);
if i < k then Series1.AddXY(i, m, '', clRed);
if i = k then Series1.AddXY(i, m, '', clFuchsia);
if i > k then Series1.AddXY(i, m, '', clBlue);
end;
end;
end;
// 初期設定
// epsilonの設定
procedure TForm1.FormCreate(Sender: TObject);
var
m, d: double;
begin
LabeledEdit1.EditLabel.Caption := 'λ ≧0' + #13#10 + '母数(平均値)';
LabeledEdit2.EditLabel.Caption := 'k=1,2,3,…' + #13#10 + 'パーセント点';
LabeledEdit3.EditLabel.Caption := 'P ≧0' + #13#10 + '下側累積確率';
Epsilon := 1;
d := 1;
repeat
Epsilon := Epsilon / 2;
m := d + Epsilon;
until d = m;
Epsilon := Epsilon * 2;
end;
end.
Poisson_distribution.zip