library Estimate;

{可靠性參數估計}

uses
  SysUtils,
  Classes,
  Math,
  u_EstimateInt in 'u_EstimateInt.pas';

{$R *.RES}


{指數分布密度函數}
function ExpDistf(lambda,t:Double):Double;StdCall;
begin
  result:=lambda*Exp(-lambda*t);
end;

{指數分布累計函數}
function ExpDistFT(lambda,t:Double):Double;StdCall;
begin
  result:=1-Exp(-lambda*t);
end;

{指數分布點估計}
function ExpDistPointEsti(WorkTime:TTimeList):Double;StdCall;
var
  i,r:integer;
  total:Double;
begin
  r:=0;
  total:=0;
  for i:=0 to length(WorkTime) do
  begin
    total:=total+WorkTime[i].Time;
    if WorkTime[i].Failed then Inc(r);
  end;
  if r>0 then
    result:=total/r
  else
    result:=-1;
end;

{威布爾分布密度函數}
function WeibullDistf(beta,eta,gamma,t:Double):Double;StdCall;
var
  temp1,temp2:Double;
begin
  temp1:=(t-gamma)/eta;
  if temp1=0 then temp1:=0.001;
  if temp1>0 then
  begin
    temp2:=Power(temp1,beta-1);
    result:=(beta/eta)*temp2*Exp(-temp2*temp1);
  end
  else
    result:=0;
end;

{威布爾分布累計函數}
function WeibullDistFT(beta,eta,gamma,t:Double):Double;StdCall;
var
  temp:Double;
begin
  temp:=(t-gamma)/eta;
  if temp>0 then
    result:=1-Exp(-Power((t-gamma)/eta,beta))
  else
    result:=0;
end;

{威布爾分布累計函數反函數}
function WeibullDistFA(beta,eta,gamma,F:Double):Double;StdCall;
begin
 result:=gamma+eta*Power(-ln(1-F),1/beta);
end;

{二參數威布爾分布點估計-最小二乘估計}
function Weibull2PointEstiLeast(WorkTime:TTimeList; var beta, eta:Double):Boolean;StdCall;
var
  TimeR:TSampleRList;
  PointList:TPointList;
  i,n:integer;
  a,b,r:Double;
begin
  n:=length(WorkTime);
  SetLength(TimeR,n);
  n:=TruncationR(WorkTime,TimeR);
  SetLength(TimeR,n);
  SetLength(PointList,n);
  for i:=0 to n-1 do
  begin
    PointList[i].X:=ln(TimeR[i].Time);
    PointList[i].Y:=ln(-ln(TimeR[i].R));
  end;
  r:=LeastSquaresFit(PointList,a,b);
  beta:=b;
  eta:=exp(-a/b);
  Result:=(r>0.999);
end;

{截尾數據可靠度計算}
function TruncationR(var WorkTime:TTimeList;CalResult:TSampleRList):integer;StdCall;
var
  i,j,n:integer;
  temp:TSample;
  flag:Boolean;
  r:Single;
begin
  //按時間排序
  flag:=True;
  j:=length(WorkTime)-2;
  while flag do
  begin
    flag:=False;
    for i:=0 to j do
    begin
      if WorkTime[i].Time>WorkTime[i+1].Time then
      begin
        temp:=WorkTime[i];
        WorkTime[i]:=WorkTime[i+1];
        WorkTime[i+1]:=temp;
        flag:=True;
      end;
    end;
    Dec(j);
  end;
  //標記記錄編號，計算R
  j:=0;
  r:=0;
  n:=length(WorkTime);
  for i:=0 to n-1 do
  begin
    if WorkTime[i].Failed then
    begin
      CalResult[j].Time:=WorkTime[i].Time;
      r:=r+(n+1-r)/(n+1-i);
      CalResult[j].R:=1-(r-0.3)/(n+0.4); //使用中位秩公式計算
      Inc(j);
    end;
  end;
  Result:=j;
end;

{最小二乘法直線擬合y=a+bx}
function LeastSquaresFit(PointList:TPointList;var a,b:Double):Double;StdCall;
var
  i,n:integer;
  sum_x,sum_y,sum_xy,sum_x2:Double;
  avg_x,avg_y:Double;
  delta_x,delta_y:Double;
  sum_delta_xy,sum_delta_x2,sum_delta_y2:Double;
begin
  n:=length(PointList);
  sum_x:=0;
  sum_y:=0;
  sum_xy:=0;
  sum_x2:=0;
  for i:=0 to n-1 do
  begin
    sum_x:=sum_x+PointList[i].X;
    sum_y:=sum_y+PointList[i].Y;
    sum_xy:=sum_xy+PointList[i].X*PointList[i].Y;
    sum_x2:=sum_x2+Power(PointList[i].X,2);
  end;
  a:=(sum_xy*sum_x-sum_y*sum_x2)/(Power(sum_x,2)-n*sum_x2);
  b:=(sum_x*sum_y-n*sum_xy)/(Power(sum_x,2)-n*sum_x2);
  avg_x:=sum_x/n;
  avg_y:=sum_y/n;
  sum_delta_xy:=0;
  sum_delta_x2:=0;
  sum_delta_y2:=0;
  for i:=0 to n-1 do
  begin
    delta_x:=PointList[i].X-avg_x;
    delta_y:=PointList[i].Y-avg_y;
    sum_delta_xy:=sum_delta_xy+delta_x*delta_y;
    sum_delta_x2:=sum_delta_x2+Power(delta_x,2);
    sum_delta_y2:=sum_delta_y2+Power(delta_y,2);
  end;
  result:=sum_delta_xy/(Sqrt(sum_delta_x2)*Sqrt(sum_delta_y2));
end;

{引出函數名}
exports
  ExpDistf, ExpDistFT, ExpDistPointEsti,
  WeibullDistf, WeibullDistFT, WeibullDistFA,
  Weibull2PointEstiLeast,
  TruncationR,
  LeastSquaresFit;
end.