You've already forked lazarus-ccr
21484fc4bd
git-svn-id: https://svn.code.sf.net/p/lazarus-ccr/svn@7925 8e941d3f-bd1b-0410-a28a-d453659cc2b4
108 lines
2.4 KiB
ObjectPascal
108 lines
2.4 KiB
ObjectPascal
program Project1;
|
|
|
|
uses
|
|
SysUtils,
|
|
spe,
|
|
MathUnit;
|
|
|
|
// Numerical recipes
|
|
|
|
function BetaCf(a,b,x: Double) : Double;
|
|
{ Kettenbruch-Entwicklung der unvollständigen Beta-Funktion }
|
|
const
|
|
ItMax = 100;
|
|
EPS = 1E-9;
|
|
var
|
|
tem, em, d : Double;
|
|
qap, qam, qab : Double;
|
|
bz, bpp, BP, bm : Double;
|
|
az, app, am, aold, ap : Double;
|
|
m : integer;
|
|
begin
|
|
am := 1.0;
|
|
bm := 1.0;
|
|
az := 1.0;
|
|
qab := a+b;
|
|
qap := a+1.0;
|
|
qam := a-1.0;
|
|
bz := 1.0 - qab*x/qap;
|
|
for m:=1 to ItMax do begin
|
|
em := m;
|
|
tem := em+em;
|
|
d := em*(b-m)*x/((qam+tem)*(a+tem));
|
|
ap := az + d*am;
|
|
BP := bz + d*bm;
|
|
d := -(a+em)*(qab+em)*x/((a+tem)*(qap+tem));
|
|
app := ap + d*az;
|
|
bpp := BP + d*bz;
|
|
aOld:= az;
|
|
am := ap/bpp;
|
|
bm := BP/bpp;
|
|
az := app/bpp;
|
|
bz := 1.0;
|
|
if abs(az-aOld) < EPS*abs(az) then begin
|
|
result := az;
|
|
Exit;
|
|
end;
|
|
end;
|
|
raise Exception.Create('Kettenbruchentwicklung der unvollst.Beta-Funktion divergiert.');
|
|
// a oder b zu groß, oder ItMax zu klein
|
|
end;
|
|
|
|
function BetaI(a,b,x: Double): Double;
|
|
{ berechnet die Unvollständige Beta-Funktion Ix(a,b) = Bx(a,b)/B(a,b),
|
|
/x a-1 b-1
|
|
wobei Bx(a,b) = | t * (1-t) dt
|
|
/0
|
|
}
|
|
var
|
|
bt : Double;
|
|
begin
|
|
if (x<0.0) then
|
|
raise EMathError.Create('Argument der unvollst.Beta-Funktion ist <0.');
|
|
if (x>1.0) then
|
|
raise EMathError.Create('Argument der unvollst.Beta-Funktion ist >1.');
|
|
if (x=0.0) or (x=1.0) then
|
|
bt := 0.0
|
|
else
|
|
bt := exp(GammaLn(a+b) - GammaLn(a) - GammaLn(b) + a*ln(x) + b*ln(1.0-x));
|
|
if x < (a+1.0)/(a+b+2.0)
|
|
then result := bt * BetaCf(a,b,x)/a
|
|
else result := 1.0 - bt*BetaCf(b,a,1.0-x)/b;
|
|
end;
|
|
|
|
|
|
const
|
|
xmin = 0.0;
|
|
xmax = 1.0;
|
|
dx = 0.2;
|
|
|
|
procedure Test(a, b: Double);
|
|
var
|
|
i: Integer;
|
|
x, y_lazStats, y_numlib, y_numrecip: Double;
|
|
begin
|
|
WriteLn('a = ', a:0:3, ' b = ', b:0:3);
|
|
WriteLn;
|
|
WriteLn('x':20, 'y(lazstats)':20, 'y(numlib)':20, 'y(Num.Recip)':20);
|
|
x := xmin;
|
|
while (x <= xmax) do begin
|
|
y_numlib := spe.betai(a, b, x);
|
|
y_lazstats := mathunit.betai(a, b, x);
|
|
y_numrecip := betai(a, b, x);
|
|
|
|
WriteLn(x:20:5, y_lazstats:20:5, y_numlib:20:5, y_numrecip:20:5);
|
|
x := x + dx;
|
|
end;
|
|
end;
|
|
|
|
begin
|
|
WriteLn('incomplete beta function');
|
|
WriteLn;
|
|
Test(0.5, 0.5);
|
|
Test(1.0, 2.0);
|
|
Test(2.0, 1.0);
|
|
ReadLn;
|
|
end.
|
|
|