二项分布和超几何分布
///
/// 经典算法单元
////
/// 主要包括:二项分布和超几何分布
/// 2005-04-18
//////经典算法,计算可靠度置信下限
/////function Classical_R(gama:double;n,num,f,kind:integer):double;
////kind=0 二项分布 kind=1 超几何分布
//gama-置信水平 num-批量 n-样本量 f-失效数 f只能为0和1
unit ClassicalU;
interface
uses
Windows, Messages, SysUtils, Classes,Forms,Dialogs,math;
/////计算一次循环的值
function Factorial(f,n:integer;x:extended):Extended;
//二项分布可靠度置信下限计算
function Binomial(n,f:integer;r:double):double;
///二项分布的样本量
function Binomial_N(f:integer;r,gama:double):integer;
////超几何分布
function Hypergeometri(gama:double;n,num,f:integer): double;
///超几何分布子过程
Function cfv(R,c: double;N,m,NF: integer): extended;
//////经典算法
function Classical_R(gama:double;n,num,f,kind:integer):double;
implementation
/////求样本量
function Binomial_N(f:integer;r,gama:double):integer;
////n 样本量 f-失效数 r-置信度
var
OldR,NewR:Extended;//1-r 的值
status0,status1:double;
Step:integer;//步长
Eps:double;//精度要求
dn:double;
n:integer;
i,j,Direction,OldD:integer;//1－正向 -1 反向 方向改变时 步长/10
begin
OldR:=1.00-gama;
NewR:=0.0;
result:=f;
Step:=1;
Eps:=1.0e-6;
Direction:=1;
OldD:=1;
if (f<0) then
begin
Application.MessageBox('输入数据不合理!', PChar(Application.Title), MB_OK +
MB_ICONINFORMATION);
Exit;
end;
/////f=0 单独的计算公式
if f=0 then
begin
dn:=ln(1-gama)/ln(r);
if Frac(dn)>0 then
result:=trunc(dn)+1
else
result:=trunc(dn);
exit;
end;
/////f>0 的情况
j:=0;
status0:=0;
status1:=0;
result:=f;
while abs(OldR-NewR)>eps do
begin
NewR:=0.0;
for i:=0 to f do
begin
NewR:=NewR+Factorial(i,result,r)
end;
//////如果方向改变 步长/10
if OldR-NewR>0 then
begin
Direction:=-1;
if Direction<>OldD then
begin
status0:=abs(oldr-newr);
end;
OldD:=-1;
end
else
begin
Direction:=1;
if Direction<>OldD then
begin
status1:=abs(oldr-newr);
end;
OldD:=1;
end;
result:=Result+step*Direction;
if (status1<status0)and(status1>0) then
begin
if result>200 then
result:=result-1;
exit;
end;
if (status0<status1)and(status0>0) then
begin
if result>200 then
result:=result+1;
exit;
end;
/////防止死循环//////
j:=j+1;
if j>100000000 then
begin
Exit;
end;
end;
end;
function Classical_R(gama:double;n,num,f,kind:integer):double;
////经典算法,计算可靠度置信下限
////kind=0 二项分布 kind=1 超几何分布
begin
result:=0.0;
if kind=0 then////二项分布
begin
// if num<10*n then exit;
result:=Binomial(n,f,gama);
end;
if kind=1 then
begin
// if num>=10*n then exit;
result:=Hypergeometri(gama,n,num,f);
end;
end;
Function cfv(R,c: double;N,m,NF: integer): extended;
/////超几何分布的子过程 c-置信度
////当 f=0,f=1 or f=0+1 时 Nf=1/2/3
///r-可靠度置信下限 N-样本量 m-批量
var
s,a,b,xm,y,d,rm,g:double;
i,j,k,mr,l:integer;
begin
xm:=m*1.0;
y:=xm*r;
g:=y;
mr:=trunc(g);
rm:=mr*1.0;
s:=1.0;
if n<=3 then
begin
if n=1 then ///样本量为1的情况
begin
s:=c-r-(1.0-r)*(nf-1);
result:=s;
exit;
end;
if n=2 then ///样本量为2的情况
begin
if xm<>1.0 then
s:=c-(R*(y-1.0)+2.0*y*(1.0-r)*(nf-1))/(xm-1.0)
else
s:=0.0;
result:=s;
Exit;
end;
if n=3 then ///样本量为3的情况
begin
if (xm<>1.0)and(xm<>2.0) then
s:=c-(y-1.0)*R*(y-2.0+3.0*xm*(1.0-r)*(nf-1))/(xm-1.0)/(xm-2.0)
else
s:=0;
result:=s;
exit;
end;
end;
////////n>3 的情况////////
l:=trunc(n-c)-2;
for k:=1 to l do
begin
d:=k*1.0;
if xm<>d then
s:=s*(rm-d)/(xm-d)*1.0
else
s:=0;
end;
if (m-n+2<>0)and (m-n+1<>0) then
s:=s/((m-n+2)*(m-n+1))*1.0
else
s:=0.0;
j:=mr-n+1;
if nf=2 then
j:=j+n*(m-mr);
i:=j+1;
if nf=2 then
i:=i-n;
b:=(mr-n+2)*j*1.0;
a:=s*b;
b:=s*rm*i*1.0;
if mr+1=m then
b:=1.0;
s:=c-((b-a)*(y-rm)+a)*r;
result:=s;
end;
function Hypergeometri(gama:double;n,num,f:integer): double;
//可靠度单侧置信下限(超几何分布)
//num-批量 n-样本量 f-失效数 gama-置信水平
var
a,b:double;
eps:double;
p,q:Double;
dr,s:Double;
i:integer;
begin
eps:=1.0e-7;
gama:=1-gama;
f:=IfThen(f=0,1,2);
a:=n/num+1.0e-6;
p:=cfv(a,gama,n,num,f);
b:=0.9999999;
q:=cfv(b,gama,n,num,f);
if p*q>0 then
begin
dr:=0.0;
Exit;
end;
s:=q;
i:=0;
while (abs(s)>eps) do
begin
dr:=0.618*(b-a)+a;
s:=cfv(dr,gama,n,num,f);
if s*p>0.0 then
a:=dr;
if s*p<=0.0 then
b:=dr;
if I>10000 then
begin
dr:=0.0;
break;;
end;
inc(i);
end;
result:=dr;
end;
/////计算一次循环的值
function Factorial(f,n:integer;x:extended):Extended;
/////f 当前失效数，n 样本量 x 当前递增（减)R可靠度值
/////计算一次循环的值
var
i:integer;
begin
if (n<=0) or (f<0) or (x<=0.0) then
begin
result:=0.0;
exit;
end;
if f=0 then
begin
result:=power(x,n);
exit;
end;
result:=1.0;
//n 和 n-x 相约以后
for i:=n-f+1 to n do
begin
result:=Result*i;
end;
result:=result*power(x,n-f);
for i:=1 to f do
begin
result:=Result/i;
result:=Result*(1-x);
end;
end;
//二项分布可靠度置信下限计算
function Binomial(n,f:integer;r:double):double;
////n 样本量 f-失效数 r-置信度
////二项分布可靠度置信下限的计算
var
OldR,NewR:Extended;//1-r 的值
Step:double;//步长
Eps:double;//精度要求
i,j,Direction,OldD:integer;//1－正向 -1 反向 方向改变时 步长/10
begin
OldR:=1.0-r;
NewR:=0.0;
Result:=0.1;
Step:=0.1;
Eps:=1.0e-6;
Direction:=1;
OldD:=1;
if (N<=0) or (n<f) or(n>99999) then
begin
Application.MessageBox('输入数据不合理!', PChar(Application.Title), MB_OK +
MB_ICONINFORMATION);
Exit;
end;
/////f=0 单独的计算公式
if f=0 then
begin
result:=power(OldR,1.0/n);
Exit;
end;
/////f>0 的情况
j:=0;
while abs(OldR-NewR)>eps do
begin
NewR:=0.0;
for i:=0 to f do
begin
NewR:=NewR+Factorial(i,n,Result)
end;
//////如果方向改变 步长/10
if OldR-NewR>0 then
begin
Direction:=1;
if Direction<>OldD then
begin
Step:=Step/10;
end;
OldD:=1;
end
else
begin
Direction:=-1;
if Direction<>OldD then
begin
Step:=Step/10;
end;
OldD:=-1;
end;
result:=Result+step*Direction;
/////防止死循环//////
j:=j+1;
if j>1000000 then
Exit;
end;
end;
end.
