Ｆ分布
Ｆ分布は統計学、確率論に用いられる連続確率分布で詳細は、ウィキペディア Ｆ分布 を参照してください。
又、Ｆ分布は重要な分布なので、インターネットで検索すれば、詳細な説明や計算式が見つかります。
プログラムは、ｃ言語辞典 アルゴリズムにあるものをDelph用に変換すると同時に、グラフの追加をしています。
同じ、値を求めるのに二種類の計算を組み込んであります。
全体の累積確率は、１なので、上側か下側を計算し、１から差し引けば、反対側の累積確率が計算出来ますが、此処ではそれぞれ独立に計算しています。
プログラム
unit Main;
interface
uses
Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes, Vcl.Graphics,
Vcl.Controls, Vcl.Forms, Vcl.Dialogs, VclTee.TeeGDIPlus, Vcl.StdCtrls, Vcl.ExtCtrls,
VCLTee.Series, VCLTee.TeEngine, VCLTee.TeeProcs, VCLTee.Chart, system.UITypes, system.Math;
type
TForm1 = class(TForm)
Chart1: TChart;
Memo1: TMemo;
LabeledEdit2: TLabeledEdit;
LabeledEdit1: TLabeledEdit;
LabeledEdit3: TLabeledEdit;
Button1: TButton;
Series1: TLineSeries;
Series2: TPointSeries;
Series3: TLineSeries;
CheckBox1: TCheckBox;
Series4: TLineSeries;
procedure Button1Click(Sender: TObject);
private
{ Private 宣言 }
public
{ Public 宣言 }
end;
var
Form1: TForm1;
implementation
{$R *.dfm}
// ************ loggamma(x) ************
const
N = 8;
B0 = 1; // 以下はBernoulli数
B1 = (-1.0 / 2.0);
B2 = ( 1.0 / 6.0);
B4 = (-1.0 / 30.0);
B6 = ( 1.0 / 42.0);
B8 = (-1.0 / 30.0);
B10 = ( 5.0 / 66.0);
B12 = (-691.0 / 2730.0);
B14 = ( 7.0 / 6.0);
B16 = (-3617.0 / 510.0);
function loggamma(x : double): double; // ガンマ関数の対数
var
LOG_2PI : Double;
v, w : Double;
begin
v := 1;
while x < N do begin
v := v * x;
x := x + 1;
end;
LOG_2PI := ln(2 * pi);
w := 1 / (x * x);
result := ((((((((B16 / (16 * 15)) * w + (B14 / (14 * 13))) * w
+ (B12 / (12 * 11))) * w + (B10 / (10 * 9))) * w
+ (B8 / ( 8 * 7))) * w + (B6 / ( 6 * 5))) * w
+ (B4 / ( 4 * 3))) * w + (B2 / ( 2 * 1))) / x
+ 0.5 * LOG_2PI - ln(v) - x + (x - 0.5) * ln(x);
end;
//**********************************************************
function beta(a, b: double): double;
begin
result := exp(loggamma(a) + loggamma(b) - loggamma(a + b));
end;
//**********************************************************
function p_beta(x, a, b : double): double;
var
k : integer;
p1, q1, p2, q2, d, previous : double;
begin
if a <= 0 then begin
result := MaxDouble;
exit;
end;
if b <= 0 then begin
if x < 1 then begin
result := 0;
exit;
end;
if x = 1 then begin
result := 1;
exit;
end;
result := MaxDouble;
exit;
end;
if x > (a + 1) / (a + b + 2) then begin
result := 1 - p_beta(1 - x, b, a);
exit;
end;
if x <= 0 then begin
result := 0;
exit;
end;
p1 := 0;
q1 := 1;
p2 := exp(a * ln(x) + b * ln(1 - x)
+ loggamma(a + b) - loggamma(a) - loggamma(b)) / a;
q2 := 1;
k := 0;
while k < 200 do begin
previous := p2;
d := - (a + k) * (a + b + k) * x
/ ((a + 2 * k) * (a + 2 * k + 1));
p1 := p1 * d + p2;
q1 := q1 * d + q2;
inc(k);
d := k * (b - k) * x / ((a + 2 * k - 1) * (a + 2 * k));
p2 := p2 * d + p1;
q2 := q2 * d + q1;
if q2 = 0 then begin
p2 := Maxdouble;
continue;
end;
p1 := p1 / q2;
q1 := q1 / q2;
p2 := p2 / q2;
q2 := 1;
if p2 = previous then begin
result := p2;
exit;
end;
end;
Form1.Memo1.Lines.Add('p_beta: 収束しません');
result := p2;
end;
function q_beta(x, a, b: double): double;
begin
result := 1 - p_beta(x, a, b);
end;
function p_f(f: double; df1, df2: integer): double; // F 分布の下側確率
begin
if f <= 0 then begin
result := 0;
exit;
end;
result := p_beta(df1 / (df1 + df2 / f), 0.5 * df1, 0.5 * df2);
// result := p_beta(df1 * f / (df1 * f + df2), 0.5 * df1, 0.5 * df2); // G(x)
end;
function q_f(f: double; df1, df2: integer): double; // F 分布の上側確率
begin
if f <= 0 then begin
result := 1;
exit;
end;
result := p_beta(df2 / (df2 + df1 * f), 0.5 * df2, 0.5 * df1);
end;
//**********************************************************
function q_x(f: double; df1, df2: integer): double; // 確率密度関数 g(x)
var
p1, p2, p3, p4 : double;
begin
if f = 0 then begin
if df1 = 1 then begin
result := infinity;
exit;
end;
if df1 div 3 = 0 then result := 1
else result := 0;
exit;
end;
p1 := 1 / beta(df1 / 2, df2 / 2);
p2 := df1 * f / (df1 * f + df2);
p3 := power(p2, df1 / 2);
p4 := power(1 - p2, df2 / 2);
result := p1 * p3 * p4 / f;
end;
//**********************************************************
function q_f2( f: double; df1, df2: integer): double; // 上側累積確率
var
i: integer;
cos2, sin2, prob, temp: double;
begin
if f <= 0 then begin
result := 1;
exit;
end;
if (df1 mod 2 <> 0) and (df2 mod 2 = 0) then begin
result := 1 - q_f(1 / f, df2, df1);
exit;
end;
cos2 := 1 / (1 + df1 * f / df2);
sin2 := 1 - cos2;
if df1 mod 2 = 0 then begin
prob := power(cos2, df2 / 2.0);
temp := prob;
i := 2;
while i < df1 do begin
temp := temp * (df2 + i - 2) * sin2 / i;
prob := prob + temp;
inc(i, 2);
end;
result := prob;
exit;
end;
prob := arctan(sqrt(df2 / (df1 * f)));
temp := sqrt(sin2 * cos2);
i := 3;
while i <= df1 do begin
prob := prob + temp;
temp := temp * (i - 1) * sin2 / i;
inc(i, 2);
end;
temp := temp * df1;
i := 3;
while i <= df2 do begin
prob := prob - temp;
temp := temp * (df1 + i - 2) * cos2 / i;
inc(i, 2);
end;
result := prob * 2 / PI;
end;
function p_f2( f: double; df1, df2: integer): double; // 下側累積確率
begin
if f <= 0 then result := 0
else result := q_f(1 / f, df2, df1);
end;
//**********************************************************
procedure TForm1.Button1Click(Sender: TObject);
var
ch : integer;
df1, df2 : integer;
x, dx, px: double;
Po, Qo : double;
fo: double;
I: Integer;
min, max: double;
begin
val(LabeledEdit1.Text, df1, Ch);
if ch <> 0 then begin
MessageDlg('V1値に間違いがあります。', mtInformation, [mbOk], 0, mbOk);
exit;
end;
if df1 <= 0 then begin
MessageDlg('V1の値をゼロより大きくしてください。', mtInformation, [mbOk], 0, mbOk);
exit;
end;
val(LabeledEdit2.Text, df2, Ch);
if ch <> 0 then begin
MessageDlg('V2値に間違いがあります。', mtInformation, [mbOk], 0, mbOk);
exit;
end;
if df2 <= 0 then begin
MessageDlg('V2の値をゼロより大きくしてください。', mtInformation, [mbOk], 0, mbOk);
exit;
end;
val(LabeledEdit3.Text, x, Ch);
if ch <> 0 then begin
MessageDlg('x値に間違いがあります。', mtInformation, [mbOk], 0, mbOk);
exit;
end;
if x < 0 then begin
MessageDlg('xの値がゼロ以下です。', mtInformation, [mbOk], 0, mbOk);
exit;
end;
Memo1.Clear;
fo := q_x(x, df1, df2);
Memo1.Lines.Add(' 確率密度 f ' + floattostr(fo));
if checkbox1.Checked = true then begin
Po := p_f2(x, df1, df2);
Qo := q_f2(x, df1, df2);
end
else begin
Po := p_f(x, df1, df2);
Qo := q_f(x, df1, df2);
end;
Memo1.Lines.Add('下側累積確率 P ' + floattostr(Po));
Memo1.Lines.Add('上側累積確率 Q ' + floattostr(Qo));
Series1.Clear;
Series2.Clear;
Series3.Clear;
Series4.Clear;
min := 0;
max := 5;
if x > max - 0.5 then max := x + 1;
dx := (max - min) / 200;
if fo < MAXSINGLE then begin
Series2.AddXY(x, fo);
Series3.AddXY(x, 0, '', clRed);
Series3.AddXY(x, fo, '', clRed);
end
else begin
fo := q_x(dx, df1, df2);
Series3.AddXY(0, 0, '', clRed);
Series3.AddXY(0, fo, '', clRed);
end;
for i := 0 to 200 do begin
px := i * dx + min;
fo := q_x(px, df1, df2);
if fo < MAXSINGLE then
if px <= x then begin
Series1.AddXY(px, fo,'', clRed);
Series4.AddXY(px, 0,'', clRed);
Series4.AddXY(px, fo,'', clRed);
Series4.AddXY(px, 0,'', clRed);
end
else begin
Series1.AddXY(px, fo,'', clBlue);
Series4.AddXY(px, 0,'', clBlue);
Series4.AddXY(px, fo,'', clBlue);
Series4.AddXY(px, 0,'', clBlue);
end;
end;
end;
end.
F-distribution.zip