BeginPackage["Reliability`"]; 

stepPlot::usage = "Plots a list of points as a step function (like \
Histogram[] \
but only showing the outline of the top of the bars). It is useful for \
plotting \
things like discrete probability distributions. It takes one option, \
PlotWith -> {f1,f2, . . .}, where the fi are other functions to be plotted. \
The \
other functions are plotted in a different color. This is useful in comparing \
\
the step plot to a fitted function."; 

PlotWith::usage = "PlotWith is an option for stepPlot, empiricalCDF, and \
empiricalHazard that specifies a list of additional functions to be plotted \
in \
the same range. These must be pure functions, e.g., Sin[#]&, not Sin[x]."; 

aggregate::usage = "Aggregates a list of data points for use by functions \
such \
as empiricalCDF. It takes one argument of the form {{1,t1},. . .,{i,ti},. . \
.,{n,tn}} \
representing n observations occurring at times t1, . . ., ti, . . ., tn. The \
output \
is a new list of the form {. . .,{ti,i},. . .} representing a discrete \
single-valued \
function of t. If there are multiple observations with the same value of t, \
the aggregation \
leaves only the observation with the largest i, for a given t.";

empiricalCDF::usage = "Plots an empirical CDF. It takes one argument of the \
form \
{{1,t1},. . .,{i,ti},. . .,{n,tn}} representing n observations occurring at \
times \
t1, . . ., ti, . . ., tn. Returns a list of points {t, p} defining the CDF. \
It takes \
two paramaters: Plot -> True causes the CDF to be plotted; PlotWith -> \
{f1,f2, . . .} \
gives a list of other functions to be plotted, as in stepPlot (q.v.).";

empiricalHazard::usage = "Plots an empirical hazard rate. It takes one \
argument of the form \
{{1,t1},. . .,{i,ti},. . .,{n,tn}} representing n observations occurring at \
times \
t1, . . ., ti, . . ., tn. Returns a list of points {t, p} defining the hazard \
function. It takes \
two paramaters: Plot -> True causes the hazard function to be plotted; \
PlotWith -> {f1,f2, . . .} \
gives a list of other functions to be plotted, as in stepPlot (q.v.).";

empiricalReliability::usage = "Returns a set of points {t, p} defining an \
empirical reliability \
or 'survivor' function R[t] = 1 - F[t] where F[t} is the CDF. It takes one \
argument of the form \
{{1,t1},. . .,{i,ti},. . .,{n,tn}} representing n observations occurring at \
times \
t1, . . ., ti, . . ., tn. There are no plotting options.";

ksStatistic::usage = "Returns the Kolmogorov-Smirnov statistic for goodness \
of fit between \
the two parameters supplied: a dataset representing an empirical CDF, of the \
form \
{ . . . {t, p}, . . .}, and a pure function representing a CDF against which \
the \
goodness of fit is to be tested."; 

weibullPDF::usage = "weibullPDF[a,b,t] returns a Weibull probability density \
function \
with scale parameter a, shape parameter b, and independent variable t.";

weibullCDF::usage = "weibullCDF[a,b,t] returns a Weibull cumulative \
distribution function \
with scale parameter a, shape parameter b, and independent variable t."; 

weibullHazard::usage = "weibullPDF[a,b,t] returns a Weibull hazard function \
\
with scale parameter a, shape parameter b, and independent variable t."; 

fitWeibull::usage = "fitWeibull[] fits (and returns) a Weibull CDF from its \
input, \
a set or points representing an empirical reliability function (e.g., the \
output from \
empiricalReliability[]. The returned value is a list {shapeParameter, \
scaleParameter, \
CDFfunction}. The function returned is a pure function, so, e.g., \
fitWeibull[data][[3]][x] \
will be a function of x.";

ksWeibullGof::usage = "Takes a list representing a failure dataset as input. \
Using a table \
of critical values for the K-S statistic for a fitted Weibull distribution, \
the goodness of \
fit is tested directly from the failure dataset. ksWeibullGof[] returns the \
raw data in a \
list: sample size, Weibull shape and scale parameters, a boolean reject \
indicating whether \
the null hypothesis is to be rejected, the significance level alpha, the \
value of the K-S \
statistic, and the fitted Weibull CDF as a pure function. Compare \
printKsWeibullGof[]."; 

printKsWeibullGof::usage = "Takes a list representing a failure dataset as \
input. Prints \
the goodness-of-fit information from ksWeibullGof[], and returns the fitted \
Weibull CDF \
as a pure function.";


Begin["`Private`"]; 

Needs["Statistics`DescriptiveStatistics`"]; 
Needs["Statistics`DataSmoothing`"]; 
Needs["Statistics`LinearRegression`"]; 
Needs["Statistics`HypothesisTests`"]; 
Needs["Statistics`ContinuousDistributions`"]; 
Off[General::spell1]; 
Off[DesignedRegress::tsos]; 
Off[DesignedRegress::rsqr]; 
Off[DesignedRegress::rank]; 
Off[Solve::ifun]; 

stepPlot[points_?ListQ, opts___?OptionQ] := 
   Module[{otherFuns := PlotWith /. Flatten[{{opts}, Options[stepPlot]}], 
     stepFun, vals, plotLim, p1, p2}, 
    stepFun[x_] := (
       vals = Select[Sort[points], (First[#1]<= x)&]; 
       If[vals === {}, 0, Last[Last[vals]]]
       ); 
    plotLim := Max[First /@ points] + 1; 
    p1 := Plot[stepFun[x], {x, 0, plotLim}, PlotRange->All, \
DisplayFunction->Identity]; 
    If[otherFuns === {}, 
      Show[p1, DisplayFunction -> $DisplayFunction], 
      (p2 := Plot[Evaluate[Map[Apply[#,{x}]&,otherFuns]],{x, 0, plotLim}, \
PlotRange->All, DisplayFunction->Identity]; 
       Show[p1, p2 /. Line[l___]->{Hue[1], Line[l]}, \
DisplayFunction->$DisplayFunction])]; 
 ]; 
Options[stepPlot] = {PlotWith -> {}}; 

aggregate[data_?ListQ]:=
  Module[{newData={},n},
    n=Length[data];
    For[i=1,i<n,i++,
      While[data[[i]][[2]]==data[[i+1]][[2]],i++];
      newData=Append[newData,data[[i]]];
      If[i==n-1,newData=Append[newData,data[[n]]]]];
    Map[{Last[#],First[#]}&,newData]];
    
empiricalCDF[data_?ListQ,opts___?OptionQ]:=
  Module[
    {plotQ:=Plot/.Flatten[{{opts},Options[empiricalCDF]}],
      otherFuns:=PlotWith/.Flatten[{{opts},Options[stepPlot]}],
      n, pFormula, distTable},
    n=Length[data];
    pFormula=(# - 0.3)/(n + 0.4)&; (* Benard formula *) 
    distTable=Map[{First[#],pFormula[Last[#]]}&,aggregate[data]];
    If[plotQ,stepPlot[distTable,PlotWith->otherFuns]; distTable,
        distTable]];
Options[empiricalCDF]={Plot->False,PlotWith->{}};

empiricalHazard[data_?ListQ,opts___?OptionQ]:=
    Module[
      {plotQ:=Plot/.Flatten[{{opts},Options[empiricalCDF]}],
        otherFuns:=PlotWith/.Flatten[{{opts},Options[stepPlot]}],
        n,aggData},
      n=Length[data]; 
      aggData=aggregate[data];  (* 
        aggData format is now {. . ., {t, i},. . .} *)
      hazardTable={};
      For[i=1,i<Length[aggData],i++,
        hazardTable=
          Append[hazardTable,{aggData[[i]][[2]],
              1/((n-aggData[[i]][[2]]+0.7)(aggData[[i+1]][[1]]-
                        aggData[[i]][[1]]))}]];
      If[plotQ,stepPlot[hazardTable,PlotWith->otherFuns];hazardTable,
        hazardTable]];
Options[empiricalHazard]={Plot->False,PlotWith->{}};

empiricalReliability[data_?ListQ]:= 
   Map[{#[[1]],1-#[[2]]}&,empiricalCDF[data,Plot->False]];

ksStatistic[data_,cdf_]:=
  Module[{dataProbs,cdfProbs, dPlus, dMinus},
    dataProbs:=Map[Last[#]&,data];
    cdfProbs:=Map[cdf[First[#]]&,data]; 
    dPlus:=Max[Abs[dataProbs-cdfProbs]];
    dMinus:=Max[Abs[Prepend[Drop[dataProbs,-1], 0]-cdfProbs]];
    Max[dPlus, dMinus]]
    
weibullPDF[a_, b_, t_] := (b/a)*(t/a)^(b - 1)*Exp[-(t/a)^b]; 
weibullCDF[a_, b_, t_] := 1 - Exp[-(t/a)^b]; 
weibullHazard[a_, b_, t_] := (b/a^b)*t^(b - 1); 

fitWeibull[data_?ListQ]:=
  Module[{logLogData,fit,intercept,shapeParm,scaleParm,a},
    logLogData:=Map[{Log[First[#]],Log[-Log[Last[#]]]}&,data];
         fit:=Fit[logLogData,{1,t},t]; 
         intercept:=fit[[1]]; 
         shapeParm :=fit[[2]][[1]];   (* Slope of linearized equation *) 
         scaleParm= a/.Solve[intercept\[Equal]-shapeParm* Log[a],a][[1]];
         {shapeParm,scaleParm,weibullCDF[scaleParm,shapeParm,#]&}];

ksWeibullGof[(data_)?ListQ] := 
Module[
{ksWeibullTable, sampleSize, sampleSizeList, fit, cdf, i, shapeParm, \
scaleParm, empCdf, ksStat, sigList,reject,alpha}, 
ksWeibullTable := {
  {10, {0.9, 0.76}, {0.95, 0.819}, {0.99, 0.944}}, 
  {20, {0.9, 0.779}, {0.95, 0.843}, {0.99, 0.973}}, 
  {50, {0.9, 0.79}, {0.95, 0.856}, {0.99, 0.988}}, 
  {100, {0.9, 0.803}, {0.95, 0.874}, {0.99, 1.007}}}; 
  sampleSize := Length[data];
  sampleSizeList := ksWeibullTable[[Length[ksWeibullTable]]]; 
  empCdf := empiricalCDF[data, Plot -> False]; 
  empCdfSize := Length[empCdf]; 
  fit := fitWeibull[empiricalReliability[data]]; 
  shapeParm := fit[[1]];
  scaleParm := fit[[2]]; 
  cdf := fit[[3]]; 
  (* Select appropriate row from K-S table for sample size *)
  For[
    i = Length[ksWeibullTable], i > 1, i--,
    If[
      sampleSize < ksWeibullTable[[i]][[1]], 
      sampleSizeList = ksWeibullTable[[i - 1]]
    ]
  ]; 
  ksStat := ksStatistic[empCdf, cdf]; 
  adjKsStat := Sqrt[sampleSize]*ksStat; 
  sigList = Null; 
  (* Select nearest higher K-S statistic value from table *)
  For[
    i = 2, i <= Length[sampleSizeList], i++, 
    If[
      adjKsStat < sampleSizeList[[i]][[2]], 
      sigList := sampleSizeList[[i]]; Break[]
    ]
  ]; 
    reject:=sigList===Null; 
    If[!reject, alpha:=1-sigList[[1]]];
    {sampleSize,shapeParm,scaleParm,reject,alpha,ksStat,cdf}
];

printKsWeibullGof[(data_)?ListQ]:=
   Module[
    {result,sampleSize,shapeParm,scaleParm,reject,ksStat,cdf},
    result:=ksWeibullGof[data];
    sampleSize:=result[[1]];
    shapeParm:=result[[2]];
    scaleParm:=result[[3]];
    reject:=result[[4]];
    alpha:=result[[5]];
    ksStat:=result[[6]];
    cdf:=result[[7]];
    If[sampleSize<10,
      Print[StringJoin["Sample size (",ToString[sampleSize],
          ") is less than 10, result may be inaccurate"]],
      Print[StringJoin["Sample size is ",ToString[sampleSize]]]
      ];
    Print[StringJoin["Weibull shape parameter = ",ToString[shapeParm]],
      StringJoin[", scale parameter = ",ToString[scaleParm]]];
    If[reject,
      Print[
        StringJoin["Kolmogorov-Smirnov statistic is ",ToString[ksStat],
          ". \nNull hypothesis can be rejected at significance level ",
          "\[Alpha]"," = 0.10"]],
      Print[StringJoin["Kolmogorov-Smirnov statistic is ",ToString[ksStat],
          ". \nNull hypothesis can be accepted at significance level ",
          "\[Alpha]"," = ",ToString[PaddedForm[alpha,{3,2}]]]]
      ];
    cdf
    ];

End[]; 
EndPackage[];