stat := module()

    export Anova1, Anova2m, Anova2s, AnovaBeta, BernoulliCDF, BernoulliPDF, 
        BernoulliS, BetaCDF, BetaP, BetaPDF, BetaS, BinomialCDF, BinomialPDF, 
        BinomialS, BivariateNormalPDF, BivariateNormalS, BoxWhisker, CauchyPDF, 
        CauchyS, ChisquareCDF, ChisquareFit, ChisquareP, ChisquarePDF, ChisquareS, 
        ClassFreq, ConfIntAvLen, ConfIntMean, ConfIntPlot, ConfIntProp, ConfIntSuc, 
        ConfIntVar, Contingency, ContinuousS, Convolution, Correlation, Craps, 
        Die, DiscUniformS, DiscreteS, DoubleExponentialPDF, ExponentialCDF, 
        ExponentialP, ExponentialPDF, ExponentialS, ExponentialSumS, FCDF, FP, 
        FPDF, FS, FastCraps, Freq, GBCDF, GBDPDF, GBP, GammaCDF, GammaP, GammaPDF, 
        GammaS, GeometricCDF, GeometricPDF, GeometricS, GraphRandWalk, GraphRandWalk2, 
        Histogram, HypergeometricPDF, HypergeometricS, KSFit, Kurtosis, LinReg, 
        Locate, LogisticPDF, LogisticS, MarginalRelFreq, Max, Mean, Median, MedianTest, 
        Min, N01S, NegBinomialCDF, NegBinomialPDF, NegBinomialS, NormTransVarS, 
        NormalCDF, NormalMeanS, NormalMedianS, NormalP, NormalPDF, NormalS, NormalSumS, 
        NormalVarianceS, Ogive, Percentile, PlotDiscCDF, PlotEmpCDF, PlotEmpPDF, 
        PlotPolyReg, PlotRunningAverage, PoissonCDF, PoissonPDF, PoissonS, PolyReg, 
        ProbHist, QQ, QQFit, RNG, RandWalk, RandWalk2, Range, RegAnal, RegBand, 
        Residuals, RunTest, RunningSum, ScatPlot, ScatPlotLine, SignTest, Skewness, 
        StDev, StemLeaf, TCDF, TP, TPDF, TS, TimePlot, UniformCDF, UniformMeanS, 
        UniformMedianS, UniformP, UniformPDF, UniformS, UniformSumS, Variance, 
        WeibullPDF, Wilcoxon, Wilcoxon2, XbarChart, MakeRandom, rng;

option package;


ChisquareFit := proc(L::list, expr::algebraic, Classes::list, Type::string)

local N, i, X, NC, NumC, var, XA, j, Q, Fr, p, pval, C_Dadj,
PlotSet, P1, P2, lastarg;
global color, blue, red, linestyle;

option `May 23, 1994 -- May 23, 1994`;
description "To perform a Chi square goodness of fi\
t test for the data in L, the distribution with CD\
F expr, and class boundaries given by Classes. The \
fourth argument, Type, must be the string D or C to\
indicate a discrete or a continuous type.";

    if type(L, listlist) then
        N := sum(L[i][2], i = 1 .. nops(L));
        X := [seq(L[i][1] $ L[i][2], i = 1 .. nops(L))]
    else N := nops(L); X := L
    fi;
    NC := nops(Classes);
    NumC := NC - 1;
    var := op(select(type, indets(expr), name));
    if nops([var]) = 0 or 1 < nops([var]) then 
        ERROR("The second argument must have exactly one varaible in it")
    fi;
    XA := convert(X, array);
    if nargs = 4 then C_Dadj := .5
    else C_Dadj := .1*10^(-7)
    fi;
    for i to N do
        if XA[i] < Classes[1] then XA[i] := Classes[1] fi;
        if Classes[NC] <= X[i] then
            XA[i] := Classes[NC] - C_Dadj
        fi
    od;
    X := sort(convert(XA, list));
    Fr := array(1 .. NumC);
    for i to NumC do Fr[i] := 0 od;
    j := 1;
    i := 1;
    while j <= N do
        if Classes[i] <= X[j] and X[j] < Classes[i + 1] then
            Fr[i] := Fr[i] + 1; j := j + 1
        else i := i + 1
        fi
    od;
    if nargs = 4 then C_Dadj := 1 else C_Dadj := 0 fi;
    p := array(1 .. NumC);
    p[1] := evalf(subs(var = Classes[2] - C_Dadj, expr));
    for i from 2 to NumC - 1 do p[i] := evalf(
        subs(var = Classes[i + 1] - C_Dadj, expr)
         - subs(var = Classes[i] - C_Dadj, expr))
    od;
    p[NumC] :=
        evalf(1 - subs(var = Classes[NumC] - C_Dadj, expr));
    Q := 0;
    for i to NumC do Q := Q + (Fr[i] - N*p[i])^2/(N*p[i]) od;
    for i to NumC do
        if N*p[i] < 4 then
            printf("Warning: Class %3.0f has an expected frequency of %f", i, N*p[i]);
            lprint()
        fi
    od;
    if nargs = 4 then C_Dadj := .5 else C_Dadj := 0 fi;
    if nargs = 4 then lastarg := NumC else lastarg := NC fi;
    PlotSet := {[[Classes[NumC] - C_Dadj, 0], [Classes[NumC] - C_Dadj,
        Fr[NumC]/(N*(Classes[NC] - Classes[NumC]))], [Classes[lastarg] + C_Dadj,
        Fr[NumC]/(N*(Classes[NC] - Classes[NumC]))],
        [Classes[lastarg] + C_Dadj, 0]]};
    for i to NumC - 1 do PlotSet := PlotSet union 
        {[[Classes[i] - C_Dadj, 0], [Classes[i] - C_Dadj,
        Fr[i]/(N*(Classes[i + 1] - Classes[i]))], [Classes[i + 1] - C_Dadj,
        Fr[i]/(N*(Classes[i + 1] - Classes[i]))],
        [Classes[i + 1] - C_Dadj, 0]]}
    od;
    P1 := plot(PlotSet, color = blue, linestyle = 3);
    if nargs = 4 then
        PlotSet := {[[Classes[NumC] - C_Dadj, 0], [Classes[NumC] - C_Dadj,
            p[NumC]/(Classes[NC] - Classes[NumC])], [Classes[NumC] + C_Dadj,
            p[NumC]/(Classes[NC] - Classes[NumC])],
            [Classes[NumC] + C_Dadj, 0]]};
        for i to NumC - 1 do PlotSet := PlotSet union 
            {[[Classes[i] - C_Dadj, 0], [Classes[i] - C_Dadj,
            p[i]/(Classes[i + 1] - Classes[i])], [Classes[i + 1] - C_Dadj,
            p[i]/(Classes[i + 1] - Classes[i])],
            [Classes[i + 1] - C_Dadj, 0]]}
        od;
        P2 := plot(PlotSet, color = red)
    else
        diff(expr, var);
        P2 := plot(%, var = Classes[1] .. Classes[NC])
    fi;
    with(plots, display);
    print(display({P1, P2}));
    pval := 1 - evalf(ChisquareCDF(NumC - 1, Q));
    printf(
        "The chi-square goodness of fit statistic is %7.4f",
        Q);
    lprint();
    printf("The number of classes is %3f", NumC);
    lprint();
    printf("The p-value for %3f", NumC - 1);
    printf(" degrees of freedom is %6.4f", pval);
    lprint();
    lprint()

end:

KSFit := proc(L::list, expr::algebraic, R::range)
local var, SL, N, KS, i, F, DOdd, DEven, Loc, P1, P2;
option `May 23, 1994 -- Zaven A. Karian`;
description "To perform a Kolmogorov-Smirnov goodness \
of fit test for the data in L, the distribution \
with CDF expr, on the range of values specified by\
R.";

    var := op(select(type, indets(expr), name));
    if nops([var]) = 0 or 1 < nops([var]) then 
        ERROR("The second argument must have exactly one varaible in \
        it")
    fi;
    SL := sort(L);
    N := nops(SL);
    KS := 0;
    for i to N do
        F := subs(var = SL[i], expr);
        DOdd := evalf(abs((i - 1)/N - F));
        DEven := evalf(abs(i/N - F));
        if KS < max(DOdd, DEven) then
            KS := max(DOdd, DEven); Loc := i
        fi
    od;
    P1 := PlotEmpCDF(SL, R);
    P2 := plot(expr, var = R, color = blue);
    with(plots, display);
    print(display({P1, P2}));
    printf("The K-S goodness of fit statistic of %7.4f", KS);
    printf(" occured at %6.4f", SL[Loc]);
    lprint();
    lprint()

end:


Anova1, proc(D::list)
local m, n, tot, SSE, SST, SSTO, N, PV, XBar, i, j, 
      GrandMean, MST, MSE;
option `July 7, 1993 -- Zaven A. Karian`;
description "One factor ANOVA";

    m := nops(D);
    n := 0;
    tot := 0;
    SSE := 0;
    SST := 0;
    SST := 0;
    N := array(1 .. m);
    XBar := array(1 .. m);
    for i to m do
        N[i] := nops(D[i]);
        XBar[i] := Mean(D[i]);
        n := n + N[i]
    od;
    for i to m do for j to N[i] do
            tot := tot + D[i][j];
            SSE := SSE + (D[i][j] - XBar[i])^2
        od
    od;
    GrandMean := tot/n;
    for i to m do SST := SST + N[i]*(XBar[i] - GrandMean)^2
    od;
    SSTO := SSE + SST;
    MST := SST/(m - 1);
    MSE := SSE/(n - m);
    PV := 1 - FCDF(m - 1, n - m, MST/MSE);
    lprint(`           Sum of`);
    lprint(`           Squares   Degrees of  Mean Sq`);
    lprint(`Source      (SS)      Freedom     (MS)      F-\
        Ratio   p-value`);
    lprint(`----------------------------------------------\
        ---------------`);
    printf("Treat.  %10.4f   %5d   %10.4f   %8.4f   %8.4f",
        SST, m - 1, MST, MST/MSE, PV);
    lprint();
    printf("Error   %10.4f   %5d   %10.4f", SSE, n - m, MSE);
    lprint();
    printf("Total   %10.4f   %5d", SSTO, n - 1);
    lprint();
    lprint(`----------------------------------------------\
        ---------------`)

end:

Anova2m := proc(D::listlist)
local a, b, c, RowMean, ColMean, CellMean, GrandMean, SSTO,
SSA, SSB, SSAB, SSE, i, j, k, MSA, MSB, MSAB, MSE, PVR, PVC,
PVInt;
option `July 5, 1993 -- Zaven A. Karian`;
description "Two factro ANOVA (multiple obs. per cell)";
    a := nops(D);
    b := nops(D[1]);
    c := nops(D[1][1]);
    RowMean := array(1 .. a);
    ColMean := array(1 .. b);
    CellMean := array(1 .. a, 1 .. b);
    GrandMean := 0;
    SSTO := 0;
    SSA := 0;
    SSB := 0;
    SSAB := 0;
    SSE := 0;
    for i to a do
        RowMean[i] := 0;
        for j to b do for k to c do
                RowMean[i] := RowMean[i] + D[i][j][k]
            od
        od;
        RowMean[i] := RowMean[i]/(b*c)
    od;
    for j to b do
        ColMean[j] := 0;
        for i to a do for k to c do
                ColMean[j] := ColMean[j] + D[i][j][k]
            od
        od;
        ColMean[j] := ColMean[j]/(a*c)
    od;
    for i to a do for j to b do
            CellMean[i, j] := 0;
            for k to c do
                CellMean[i, j] := CellMean[i, j] + D[i][j][k]
            od;
            CellMean[i, j] := CellMean[i, j]/c
        od
    od;
    for i to a do for j to b do for k to c do
                GrandMean := GrandMean + D[i][j][k]
            od
        od
    od;
    GrandMean := GrandMean/(a*b*c);
    for i to a do for j to b do for k to c do
                SSTO := SSTO + (D[i][j][k] - GrandMean)^2;
                SSE := SSE + (D[i][j][k] - CellMean[i, j])^2
            od
        od
    od;
    for i to a do SSA := SSA + (RowMean[i] - GrandMean)^2 od;
    for j to b do SSB := SSB + (ColMean[j] - GrandMean)^2 od;
    SSA := b*c*SSA;
    SSB := a*c*SSB;
    SSAB := SSTO - SSA - SSB - SSE;
    MSA := SSA/(a - 1);
    MSB := SSB/(b - 1);
    MSAB := SSAB/((a - 1)*(b - 1));
    MSE := SSE/(a*b*(c - 1));
    PVR := 1 - FCDF(a - 1, a*b*(c - 1), MSA/MSE);
    PVC := 1 - FCDF(b - 1, a*b*(c - 1), MSB/MSE);
    PVInt := 1 - FCDF((a - 1)*(b - 1), a*b*(c - 1), MSAB/MSE)
        ;
    lprint(`           Sum of`);
    lprint(`           Squares  Degrees of   Mean Sq`);
    lprint(`Source      (SS)     Freedom      (MS)     F-R\
        atio    p-value`);
    lprint(`----------------------------------------------\
        ---------------`);
    printf("Row(A)  %10.4f   %5d   %10.4f   %8.4f   %8.4f",
        SSA, a - 1, MSA, MSA/MSE, PVR);
    lprint();
    printf("Col(B)  %10.4f   %5d   %10.4f   %8.4f   %8.4f",
        SSB, b - 1, MSB, MSB/MSE, PVC);
    lprint();
    printf("Int(AB) %10.4f   %5d   %10.4f   %8.4f   %8.4f",
        SSAB, (a - 1)*(b - 1), MSAB, MSAB/MSE, PVInt);
    lprint();
    printf("Error   %10.4f   %5d   %10.4f", SSE, a*b*(c - 1),
        MSE);
    lprint();
    printf("Total   %10.4f   %5d", SSTO, a*b*c - 1);
    lprint();
    lprint(`----------------------------------------------\
        ---------------`)

end:


Anova2s := proc(D::listlist)
local a, b, tot, SSTO, SSA, SSB, SSE, RowMean, ColMean, i, j,
GrandMean, MSA, MSB, MSE, PVR, PVC;
option `July 5, 1993 -- Zaven A. Karian`;
description "Two factor ANOVA (one obs. per cell)";
    a := nops(D);
    b := nops(D[1]);
    tot := 0;
    SSTO := 0;
    SSA := 0;
    SSB := 0;
    SSE := 0;
    RowMean := array(1 .. a);
    ColMean := array(1 .. b);
    for j to b do ColMean[j] := 0 od;
    for i to a do for j to b do tot := tot + D[i][j] od od;
    GrandMean := tot/(a*b);
    for i to a do RowMean[i] := Mean(D[i]) od;
    for j to b do
        for i to a do ColMean[j] := ColMean[j] + D[i][j] od;
        ColMean[j] := ColMean[j]/a
    od;
    for i to a do for j to b do
            SSTO := SSTO + (D[i][j] - GrandMean)^2
        od
    od;
    for i to a do SSA := SSA + (RowMean[i] - GrandMean)^2 od;
    for j to b do SSB := SSB + (ColMean[j] - GrandMean)^2 od;
    SSA := b*SSA;
    SSB := a*SSB;
    SSE := SSTO - SSA - SSB;
    MSA := SSA/(a - 1);
    MSB := SSB/(b - 1);
    MSE := SSE/((a - 1)*(b - 1));
    PVR := 1 - FCDF(a - 1, (a - 1)*(b - 1), MSA/MSE);
    PVC := 1 - FCDF(b - 1, (a - 1)*(b - 1), MSB/MSE);
    lprint(`           Sum of`);
    lprint(`           Squares   Degrees of  Mean Sq`);
    lprint(`Source      (SS)      Freedom     (MS)     F-R\
        atio    p-value`);
    lprint(`----------------------------------------------\
        ---------------`);
    printf("Row(A)  %10.4f   %5d   %10.4f   %8.4f   %8.4f",
        SSA, a - 1, MSA, MSA/MSE, PVR);
    lprint();
    printf("Col(B)  %10.4f   %5d   %10.4f   %8.4f   %8.4f",
        SSB, b - 1, MSB, MSB/MSE, PVC);
    lprint();
    printf("Error   %10.4f   %5d   %10.4f", SSE,
        (a - 1)*(b - 1), MSE);
    lprint();
    printf("Total   %10.4f   %5d", SSTO, a*b - 1);
    lprint();
    lprint(`----------------------------------------------\
        ---------------`)

end:


AnovaBeta := proc()
local n, k, xbar, xlist, ylist, alphahat, betahat, sigma2hat,
betahat1, betahat2, ssto, ssr, sse, msr, mse, f, p;
option `Mike, Andrew, and Christine were here.`;
description "Test beta = 0.";
    if nargs < 1 or 2 < nargs then
        ERROR("Must have list of lists, or (list,list)")
    fi;
    n := nops(args[1]);
    if nargs = 1 then
        xlist := [seq(args[1][k][1], k = 1 .. n)];
        ylist := [seq(args[1][k][2], k = 1 .. n)]
    else xlist := args[1]; ylist := args[2]
    fi;
    xbar := Mean(xlist);
    alphahat := Mean(ylist);
    k := 'k';
    betahat1 := sum(ylist[k]*(xlist[k] - xbar), k = 1 .. n);
    betahat2 := sum((xlist[k] - xbar)^2, k = 1 .. n);
    betahat := betahat1/betahat2;
    ssto := sum((ylist[k] - alphahat)^2, k = 1 .. n);
    ssr := sum(betahat^2*(xlist[k] - xbar)^2, k = 1 .. n);
    sse := ssto - ssr;
    msr := ssr;
    mse := sse/(n - 2);
    f := msr/mse;
    p := 1 - FCDF(1, n - 2, f);
    lprint(`----------------------------------------------\
        -----------------`);
    lprint(`              Sum of`);
    lprint(`              Squares  Degr of    Mean Sq`);
    lprint(`  Source       (SS)    Freedom      (MS)     F\
        -Ratio   p-value`);
    lprint(`----------------------------------------------\
        -----------------`);
    printf("Regression %10.4f %5.0f %13.4f %10.4f  %8.5f",
        ssr, 1, msr, f, p);
    lprint();
    printf("  Error %13.4f %5.0f %13.4f", sse, n - 2, mse);
    lprint();
    printf("  Total %13.4f %5.0f", ssto, n - 1);
    lprint();
    lprint(`----------------------------------------------\
        -----------------`)

end:


MedianTest := proc(X::list(realcons), Y::list(realcons))
local S, x, y, i, k;
option
`Joshua D. Levy  ~  June 17, 1994, modified 6/1/96-zak`;
description "Median test";
    if nargs = 2 then
        if type(nops(X) + nops(Y), even) then
            k := 1/2*nops(X) + 1/2*nops(Y)
        else ERROR("sum of numbers of observations not even")
        fi
    elif type(args[3], posint) then k := args[3]
    else ERROR("bad argurment")
    fi;
    S := sort([op(map(x, evalf(X))), op(map(y, evalf(Y)))],
        (a, b) -> evalb(op(a) < op(b)));
    x := 1;
    y := 0;
    S := eval(S);
    sum(S[i], i = 1 .. k)

end:


RunTest := proc(X::list(realcons), Y::{realcons, list(realcons)})
local S, x, y, i, runs;
option
`Joshua D. Levy  ~  June 17, 1994, modified 6/1/96-zak`;
description "Run test. Number of runs in sorted, combined l\
ist of X and Y -or- number of runs above or below given nu\
mber";
    if type(Y, realcons) then
        S := map(unapply('sign'(x - evalf(Y)), x), X)
    else
        S := sort(
            [op(map(x, evalf(X))), op(map(y, evalf(Y)))],
            (a, b) -> evalb(op(a) < op(b)));
        x := 0;
        y := 1;
        S := eval(S)
    fi;
    runs := 1;
    for i to nops(S) - 1 do
        if S[i] <> S[i + 1] then runs := runs + 1 fi
    od;
    runs

end:


SignTest := proc(C::boolean)
local i, k, vars;
option
`Joshua D. Levy  ~  June 17, 1994, modified 6/1/96-zak`;
description
"Number of elements in list(s) for which C is true";
    if not type([args[2 .. nargs]], list(name = list)) then
        ERROR("all arguments beyond first must be of the f\
        orm name=list")
    fi;
    vars := select(type, indets(C), name);
    if vars = {} then
        ERROR("test expression contains no variable")
    fi;
    if vars <> indets([args[2 .. nargs]]) then ERROR("test\
         expression variable(s) must be defined in argumen\
        ts")
    fi;
    subs(true = 1, false = 0, convert([seq(evalb(evalf(subs(
        seq(lhs(args[i]) = rhs(args[i])[k], i = 2 .. nargs),
        C))), k = 1 .. nops(rhs(args[2])))], `+`))

end:


Wilcoxon := proc(L::list(realcons), m::realcons)
local S, n, w, x, i, j, k;
option `Joshua D. Levy  ~  June 17, 1994`;
description "Wilcoxon test for one median.";
    S := sort([seq(evalf(x - m), x = L)],
        (a, b) -> evalb(abs(a) < abs(b)));
    for i while S[i] = 0 do  od;
    S := S[i .. nops(S)];
    n := nops(S);
    w := 0;
    for i to n do
        for j from i to n - 1 while abs(S[j]) = abs(S[j + 1])
        do
        od;
        if j = i then w := w + sign(S[i])*i
        else
            w :=
                w + 1/2*(i + j)*sum('sign(S[k])', k = i .. j)
                ;
            i := j
        fi
    od;
    w

end:


Wilcoxon2 := proc(X::list(realcons), Y::list(realcons))
local S, n, x, y, w, i, j, k;
option
`Joshua D. Levy  ~  June 17, 1994, modified 6/1/96-zak`;
description "Wilcoxon test for equality of two medians. Ret\
urns sum of ranks of Y in the sorted combination of X and \
Y";
    S := sort([op(map(x, evalf(X))), op(map(y, evalf(Y)))],
        (a, b) -> evalb(op(a) < op(b)));
    n := nops(S);
    w := 0;
    for i to n do
        for j from i to n - 1 while op(S[j]) = op(S[j + 1])
        do
        od;
        if j = i then
            w := w + eval(subs(x = 0, y = 1, S[i]))*i
        else
            w := w + 1/2*(i + j)*eval(
                subs(x = 0, y = 1, sum(S[k], k = i .. j)));
            i := j
        fi
    od;
    w

end:


Contingency := proc(L::list(list(nonnegint)))
local h, k, n, npr, npc, q, c, i, j;
option `Joshua D. Levy  ~  June 21, 1994`;
description "Contingency table";
    k := nops(L);
    h := nops(L[1]);
    n := sum('convert(L[i], `+`)', i = 1 .. k);
    npr := [seq(sum(L[i][j], j = 1 .. h), i = 1 .. k)];
    npc := [seq(sum(L['i'][j], 'i' = 1 .. k), j = 1 .. h)];
    q := 0;
    lprint();
    for i to k do
        for j to h do printf("%7.0f ", L[i][j]) od;
        printf("  %7.0f\n ", npr[i]);
        for j to h do
            c := npr[i]*npc[j]/n;
            q := q + evalf((L[i][j] - c)^2/c);
            printf("(%6.2f)", c)
        od;
        lprint()
    od;
    lprint();
    for j to h do printf("%7.0f ", npc[j]) od;
    printf(" %7.0f", n);
    lprint();
    lprint();
    printf("Chi-square statsistic (%.0f degrees of freedom\
        ): %5g\n", (h - 1)*(k - 1), q);
    printf("p-value: %.6g\n",
        1 - ChisquareCDF((h - 1)*(k - 1), q));
    lprint()

end:


Convolution := proc(X1::list, X2::list)
local dim1, dim2, l1, l2, j1, j2, i, L1, L2, N1, N2, A, k,
SUM, j, f1, f2, t;
option `September 7, 1993 -- Zaven A. Karian`;
description
"Computation of probabilities via the convolution formula";
    j1 := sum(X1[2*i], i = 1 .. 1/2*nops(X1));
    j2 := sum(X2[2*i], i = 1 .. 1/2*nops(X2));
    if .1*10^(-6) < abs(j1 - 1) or .1*10^(-6) < abs(j2 - 1)
    then ERROR("Probabilities must add up to 1")
    fi;
    for i to 1/2*nops(X1) do
        if X1[2*i - 1] < 0 then ERROR(
            "Random variable values must be non-negative")
        fi
    od;
    for i to 1/2*nops(X2) do
        if X2[2*i - 1] < 0 then ERROR(
            "Random variable values must be non-negative")
        fi
    od;
    dim1 := nops(X1) + 2*X1[1];
    dim2 := nops(X2) + 2*X2[1];
    l1 := array(1 .. dim1);
    l2 := array(1 .. dim2);
    j1 := 1;
    j2 := 1;
    for i from 0 to X1[1] - 1 do
        l1[j1] := 0; l1[j1 + 1] := 0; j1 := j1 + 2
    od;
    for i to nops(X1) do l1[j1] := X1[i]; j1 := j1 + 1 od;
    for i from 0 to X2[1] - 1 do
        l2[j2] := 0; l2[j2 + 1] := 0; j2 := j2 + 2
    od;
    for i to nops(X2) do l2[j2] := X2[i]; j2 := j2 + 1 od;
    L1 := [seq(l1[i], i = 1 .. j1 - 1)];
    L2 := [seq(l2[i], i = 1 .. j2 - 1)];
    N1 := 1/2*nops(L1);
    N2 := 1/2*nops(L2);
    A := array(1 .. 2*N1 + 2*N2 - 2);
    for k from 0 to N1 + N2 - 2 do
        SUM := 0;
        for j from 0 to k do
            if j < 0 or N1 - 1 < j then f1 := 0
            else f1 := L1[2*j + 2]
            fi;
            if k - j < 0 or N2 - 1 < k - j then f2 := 0
            else f2 := L2[2*k - 2*j + 2]
            fi;
            SUM := SUM + f1*f2
        od;
        A[2*k + 1] := k;
        A[2*k + 2] := SUM
    od;
    j := 0;
    while A[j + 2] = 0 do j := j + 2 od;
    [seq(A[t], t = j + 1 .. 2*N1 + 2*N2 - 2)]

end:


Craps := proc()
local A, D, j, t, point, roll, game;
option `August 23, 1993 -- Zaven A. Karian`;
description "simulation of a single craps game";
    A := array(1 .. 100);
    j := 2;
    D := Die(6, 2);
    point := D[1] + D[2];
    A[j] := point;
    j := j + 1;
    if point = 7 or point = 11 then game := 1
    elif point = 2 or point = 3 or point = 12 then game := -1
    else game := 0
    fi;
    while game = 0 do
        D := Die(6, 2);
        roll := D[1] + D[2];
        A[j] := roll;
        j := j + 1;
        if roll = point then game := 1
        elif roll = 7 then game := -1
        fi
    od;
    if game = -1 then A[1] := 0 else A[1] := 1 fi;
    [seq(A[t], t = 1 .. j - 1)]

end:


FastCraps := proc(n::posint)
local A, D, i, point, roll, game;
option `August 23, 1993 -- Zaven A. Karian`;
description "simulation of n craps games";
    A := array(1 .. n);
    for i to n do
        D := Die(6, 2);
        point := D[1] + D[2];
        if point = 7 or point = 11 then game := 1
        elif point = 2 or point = 3 or point = 12 then
            game := -1
        else game := 0
        fi;
        while game = 0 do
            D := Die(6, 2);
            roll := D[1] + D[2];
            if roll = point then game := 1
            elif roll = 7 then game := -1
            fi
        od;
        if game = -1 then A[i] := 0 else A[i] := 1 fi
    od;
    [seq(A[i], i = 1 .. n)]

end:


GraphRandWalk := proc(pn::numeric, ps::numeric, pe::numeric, steps::posint)
local A, loc, r, j, L;
option `July 5, 1993 -- Zaven A. Karian`;
description "Graphic display of random walk";
    A := array(1 .. 2*steps);
    loc := array(1 .. 2);
    loc[1] := 0;
    loc[2] := 0;
    for j to steps do
        r := rng();
        if r <= pn then loc[2] := loc[2] + 1
        elif r <= pn + ps then loc[2] := loc[2] - 1
        elif r <= pn + ps + pe then loc[1] := loc[1] + 1
        else loc[1] := loc[1] - 1
        fi;
        A[2*j - 1] := loc[1];
        A[2*j] := loc[2]
    od;
    L := [seq([A[2*j - 1], A[2*j]], j = 1 .. steps)];
    plot(L)

end:


GraphRandWalk2 := proc(pn::numeric, pe::numeric, steps::posint)
local A, loc, r1, r2, j, L;
option `July 5, 1993 -- Zaven A. Karian`;
description "Graphic display of random walk";
    A := array(1 .. 2*steps);
    loc := array(1 .. 2);
    loc[1] := 0;
    loc[2] := 0;
    for j to steps do
        r1 := rng();
        r2 := rng();
        if r1 <= pn then loc[2] := loc[2] + 1
        else loc[2] := loc[2] - 1
        fi;
        if r2 <= pe then loc[1] := loc[1] + 1
        else loc[1] := loc[1] - 1
        fi;
        A[2*j - 1] := loc[1];
        A[2*j] := loc[2]
    od;
    L := [seq([A[2*j - 1], A[2*j]], j = 1 .. steps)];
    plot(L)

end:


MarginalRelFreq := proc(A::listlist)
local NA, rmin, rmax, cmin, cmax, i, r, c, RCOUNT, CCOUNT,
rcount, ccount, RC, CC;
option `July 5, 1993 -- Zaven A. Karian`;
description
"Marginal relative frequencies of the list-of-lists, A";
    NA := nops(A);
    rmin := A[1][2];
    rmax := A[1][2];
    cmin := A[1][1];
    cmax := A[1][1];
    for i from 2 to NA do
        if A[i][2] < rmin then rmin := A[i][2] fi;
        if rmax < A[i][2] then rmax := A[i][2] fi;
        if A[i][1] < cmin then cmin := A[i][1] fi;
        if cmax < A[i][1] then cmax := A[i][1] fi
    od;
    if rmin < cmin then cmin := rmin else rmin := cmin fi;
    if cmax < rmax then cmax := rmax else rmax := cmax fi;
    RCOUNT := array(rmin .. rmax);
    CCOUNT := array(cmin .. cmax);
    for r from rmin to rmax do
        rcount := 0;
        for i to NA do
            if A[i][2] = r then rcount := rcount + 1 fi
        od;
        RCOUNT[r] := rcount
    od;
    for c from cmin to cmax do
        ccount := 0;
        for i to NA do
            if A[i][1] = c then ccount := ccount + 1 fi
        od;
        CCOUNT[c] := ccount
    od;
    RC := [seq(RCOUNT[r]/NA, r = rmin .. rmax)];
    CC := [seq(CCOUNT[c]/NA, c = cmin .. cmax)];
    [CC, RC]

end:


RandWalk := proc(pn::numeric, ps::numeric, pe::numeric,steps::posint, n::posint)
local A, loc, i, j, r;
option `July 5, 1993 -- Zaven A. Karian`;
description "Path of random walk";
    A := array(1 .. n);
    loc := array(1 .. 2);
    for i to n do
        loc[1] := 0;
        loc[2] := 0;
        for j to steps do
            r := rng();
            if r <= pn then loc[2] := loc[2] + 1
            elif r <= pn + ps then loc[2] := loc[2] - 1
            elif r <= pn + ps + pe then loc[1] := loc[1] + 1
            else loc[1] := loc[1] - 1
            fi
        od;
        A[i] := [loc[1], loc[2]]
    od;
    [seq(A[i], i = 1 .. n)]

end:


RandWalk2 := proc(pn::numeric, pe::numeric, steps::posint, n::posint)
local r1, r2, A, loc, j, i;
option `July 5, 1993 -- Zaven A. Karian`;
description "Path of random walk";
    A := array(1 .. n);
    loc := array(1 .. 2);
    for i to n do
        loc[1] := 0;
        loc[2] := 0;
        for j to steps do
            r1 := rng();
            r2 := rng();
            if r1 <= pn then loc[2] := loc[2] + 1
            else loc[2] := loc[2] - 1
            fi;
            if r2 <= pe then loc[1] := loc[1] + 1
            else loc[1] := loc[1] - 1
            fi
        od;
        A[i] := [loc[1], loc[2]]
    od;
    [seq(A[i], i = 1 .. n)]

end:


ConfIntAvLen := proc(LL::listlist)
local n, Len, i;
option `June 13, 1994 -- Zaven A. Karian`;
description "Average length of confidence intervals";
    n := nops(LL);
    Len := 0;
    for i to n do Len := Len + LL[i][2] - LL[i][1] od;
    Len/n

end:


ConfIntMean := proc()
local LL, z, L, A, i, n, m, a, s, S, alpha, ConfLev, ptr;
option `February 11, 1994 -- Zaven A. Karian`;
description
"Compute end-points of confidence intervals for mu";
    if not type(args[2], numeric) or args[2] < 0 or
    100 < args[2] then ERROR("Second argument must be a pe\
        rcentage; i.e. between 0 and 100")
    fi;
    if nargs = 3 and args[3] <= 0 then ERROR(
        "Third argument, when present, must be positive")
    fi;
    LL := args[1];
    ConfLev := args[2];
    if nargs = 3 then S := sqrt(args[3]) fi;
    n := nops(LL);
    m := nops(LL[1]);
    alpha := 1/2 + 1/200*ConfLev;
    if nargs = 3 or 24 < m then z := NormalP(0, 1, alpha)
    else z := TP(m - 1, alpha)
    fi;
    A := array(1 .. 2*n);
    ptr := 1;
    for i to n do
        L := LL[i];
        a := Mean(L);
        if nargs = 2 then s := StDev(L) else s := S fi;
        A[ptr] := evalf(a - z*s/sqrt(m));
        A[ptr + 1] := evalf(a + z*s/sqrt(m));
        ptr := ptr + 2
    od;
    [seq([A[2*i - 1], A[2*i]], i = 1 .. n)]

end:


ConfIntPlot := proc(LL::listlist, Value::realcons)
local n, A, i;
option `February 11, 1994 -- Zaven A. Karian`;
description "Graphic display of confidence intervals";
    n := nops(LL);
    A := {};
    for i to n do
        A := A union {[[LL[i][1], i], [LL[i][2], i]]}
    od;
    if nargs = 2 then A := A union {[[Value, 1], [Value, n]]}
    fi;
    plot(A, color = blue)

end:


ConfIntProp := proc(LL::listlist, ConfLev::realcons)
local L, A, i, j, n, m, alpha, z, p, f, ptr;
option `February 11, 1994 -- Zaven A. Karian`;
description "Compute confidence intervals for p";
    if ConfLev < 0 or 100 < ConfLev then ERROR("Second arg\
        ument must be a percentage; i.e. between 0 and 100\
        ")
    fi;
    n := nops(LL);
    m := nops(LL[1]);
    alpha := 1/2 + 1/200*ConfLev;
    z := NormalP(0, 1, alpha);
    A := array(1 .. 2*n);
    ptr := 1;
    for i to n do
        L := LL[i];
        p := sum(L[j], j = 1 .. m)/m;
        f := z*sqrt(p*(1 - p)/m);
        A[ptr] := evalf(p - f);
        A[ptr + 1] := evalf(p + f);
        ptr := ptr + 2
    od;
    [seq([A[2*i - 1], A[2*i]], i = 1 .. n)]

end:


ConfIntSuc := proc(L::listlist, Value::realcons)
local n, count, i, val, LL;
option `June 13, 1994 -- Zaven A. Karian`;
description
"number of confidence intervals that contain value";
    val := evalf(Value);
    LL := evalf(L);
    n := nops(LL);
    count := 0;
    for i to n do
        if val <= LL[i][2] and LL[i][1] <= val then
            count := count + 1
        fi
    od;
    count

end:


ConfIntVar := proc(LL::listlist, ConfLev::realcons)
local L, A, i, n, m, ss, ptr, chi1, chi2, CHI1, CHI2;
option `February 11, 1994 -- Zaven A. Karian`;
description "Compute confidence intervals for sigma**2";
    if ConfLev < 0 or 100 < ConfLev then ERROR("Second arg\
        ument must be a percentage; i.e. between 0 and 100\
        ")
    fi;
    n := nops(LL);
    m := nops(LL[1]);
    chi1 := 1/2 - 1/200*ConfLev;
    chi2 := 1 - chi1;
    CHI1 := (m - 1)/ChisquareP(m - 1, chi1);
    CHI2 := (m - 1)/ChisquareP(m - 1, chi2);
    A := array(1 .. 2*n);
    ptr := 1;
    for i to n do
        L := LL[i];
        ss := Variance(L);
        A[ptr] := ss*CHI2;
        A[ptr + 1] := ss*CHI1;
        ptr := ptr + 2
    od;
    [seq([A[2*i - 1], A[2*i]], i = 1 .. n)]

end:


Correlation := proc()
local X, Y, N, SX, SY, XY, SumX, SumY, i, CXY;
option `August 12, 1993 -- Zaven A. Karian`;
description "Mean correlation of the lists X and Y";
    if nargs = 1 and not type(args[1], listlist) then ERROR(
        "When one argument is used, it must be a list of l\
        ists")
    fi;
    if nargs = 2 and
    not (type(args[1], list) and type(args[2], list)) then
        ERROR(
        "When two arguments are used, both must be lists")
    fi;
    if nargs = 2 and nops(args[1]) <> nops(args[2]) then
        ERROR("The two lists must be of the same length")
    fi;
    if nargs = 1 then
        X := [seq(args[1][i][1], i = 1 .. nops(args[1]))];
        Y := [seq(args[1][i][2], i = 1 .. nops(args[1]))]
    else X := args[1]; Y := args[2]
    fi;
    N := nops(X);
    SX := StDev(X);
    SY := StDev(Y);
    XY := 0;
    SumX := 0;
    SumY := 0;
    for i to N do
        XY := XY + X[i]*Y[i];
        SumX := SumX + X[i];
        SumY := SumY + Y[i]
    od;
    CXY := (N*XY - SumX*SumY)/(N*(N - 1));
    evalf(CXY/(SY*SX))

end:


LinReg := proc()
local X, Y, N, XY, SumX, SumY, i, CXY, slope, const, x;
option `August 12, 1993 -- Zaven A. Karian`;
description "Slope and intercept of regression line";
    if nargs = 2 and not type(args[1], listlist) then ERROR(
        "When tow arguments are used the first must be a l\
        ist of lists")
    fi;
    if nargs = 3 and
    not (type(args[1], list) and type(args[2], list)) then
        ERROR("When three arguments are used, the first tw\
        o must be lists")
    fi;
    if nargs = 3 and nops(args[1]) <> nops(args[2]) then
        ERROR("The two lists must be of the same length")
    fi;
    if nargs = 2 then
        X := [seq(args[1][i][1], i = 1 .. nops(args[1]))];
        Y := [seq(args[1][i][2], i = 1 .. nops(args[1]))];
        x := args[2]
    else X := args[1]; Y := args[2]; x := args[3]
    fi;
    N := nops(X);
    XY := 0;
    SumX := 0;
    SumY := 0;
    for i to N do
        XY := XY + X[i]*Y[i];
        SumX := SumX + X[i];
        SumY := SumY + Y[i]
    od;
    CXY := evalf((N*XY - SumX*SumY)/(N*(N - 1)));
    slope := CXY/Variance(X);
    const := Mean(Y) - Mean(X)*slope;
    const + slope*x

end:


PlotPolyReg := proc()
local X, Y, A, i, k, n, B, C, j, s, BI, eqn, x, pts, c, m;
option `January 27, 1994 -- Zaven A. Karian`;
description
"Produce scatter plot of X-Y pairs with poly. reg. line";
    if nargs = 2 and not type(args[1], listlist) then ERROR(
        "When two arguments are used, the first must be a \
        list of lists")
    fi;
    if nargs = 3 and
    not (type(args[1], list) and type(args[2], list)) then
        ERROR("When three arguments are used, the first tw\
        o must be lists")
    fi;
    if nargs = 3 and nops(args[1]) <> nops(args[2]) then
        ERROR("The two lists must have the same length")
    fi;
    Digits := 50;
    if nargs = 2 then
        X := [seq(args[1][i][1], i = 1 .. nops(args[1]))];
        Y := [seq(args[1][i][2], i = 1 .. nops(args[1]))];
        m := args[2]
    else X := args[1]; Y := args[2]; m := args[3]
    fi;
    n := nops(X);
    B := array(1 .. m + 1, 1 .. m + 1);
    C := array(1 .. m + 1, 1 .. 1);
    B[1, 1] := n;
    for i to m + 1 do for j from i to m + 1 do
            if j = 1 then j := j + 1 fi;
            s := 0;
            s := convert([seq(X[k]^(i + j - 2), k = 1 .. n)],
                `+`);
            B[i, j] := s;
            B[j, i] := s
        od
    od;
    C[1, 1] := convert(Y, `+`);
    for i from 2 to m + 1 do
        s := 0;
        s := convert([seq(Y[k]*X[k]^(i - 1), k = 1 .. n)],
            `+`);
        C[i, 1] := s
    od;
    BI := linalg[inverse](B);
    A := linalg[multiply](BI, C);
    eqn := 0;
    for i to m + 1 do eqn := eqn + A[i, 1]*x^(i - 1) od;
    pts := ScatPlot(X, Y);
    c := plot(eqn, x = Min(X) .. Max(X));
    with(plots, display);
    display({pts, c})

end:


PolyReg := proc()
local X, Y, A, i, n, B, C, j, s, BI, eqn, x, k, m;
option `January 27, 1994 -- Zaven A. Karian`;
description
"Produce scatter plot of X-Y pairs with poly. reg. line";
    if nargs = 3 and not type(args[1], listlist) then ERROR(
        "When three arguments are used, the first must be \
        a list of lists")
    fi;
    if nargs = 4 and
    not (type(args[1], list) and type(args[2], list)) then
        ERROR("When four arguments are used, the first two\
         must be lists")
    fi;
    if nargs = 4 and nops(args[1]) <> nops(args[2]) then
        ERROR("The two lists must have the same length")
    fi;
    Digits := 50;
    if nargs = 3 then
        X := [seq(args[1][i][1], i = 1 .. nops(args[1]))];
        Y := [seq(args[1][i][2], i = 1 .. nops(args[1]))];
        m := args[2];
        x := args[3]
    else
        X := args[1];
        Y := args[2];
        m := args[3];
        x := args[4]
    fi;
    n := nops(X);
    B := array(1 .. m + 1, 1 .. m + 1);
    C := array(1 .. m + 1, 1 .. 1);
    B[1, 1] := n;
    for i to m + 1 do for j from i to m + 1 do
            if j = 1 then j := j + 1 fi;
            s := 0;
            s := convert([seq(X[k]^(i + j - 2), k = 1 .. n)],
                `+`);
            B[i, j] := s;
            B[j, i] := s
        od
    od;
    C[1, 1] := convert(Y, `+`);
    for i from 2 to m + 1 do
        s := 0;
        s := convert([seq(Y[k]*X[k]^(i - 1), k = 1 .. n)],
            `+`);
        C[i, 1] := s
    od;
    BI := linalg[inverse](B);
    A := linalg[multiply](BI, C);
    eqn := 0;
    for i to m + 1 do eqn := eqn + A[i, 1]*x^(i - 1) od;
    evalf(eqn, 10)

end:


RegAnal := proc()
local n, k, xbar, xlist, ylist, CL, c1, c2, c3, alphahat,
betahat, sigma2hat, betahat1, betahat2, alphaconf, betaconf,
sigma2conf;
option `Mike, Andrew, and Christine were here.`;
description "Regression Analysis";
    if nargs <= 1 or 3 < nargs then ERROR("Must have list \
        of lists and confidence level list, or (list,list,\
        confidence list)")
    fi;
    n := nops(args[1]);
    if nargs = 2 then
        xlist := [seq(args[1][k][1], k = 1 .. n)];
        ylist := [seq(args[1][k][2], k = 1 .. n)];
        CL := args[2]
    else xlist := args[1]; ylist := args[2]; CL := args[3]
    fi;
    xbar := Mean(xlist);
    alphahat := Mean(ylist);
    k := 'k';
    betahat1 := sum(ylist[k]*(xlist[k] - xbar), k = 1 .. n);
    betahat2 := sum((xlist[k] - xbar)^2, k = 1 .. n);
    betahat := betahat1/betahat2;
    sigma2hat := sum(
        (ylist[k] - alphahat - betahat*(xlist[k] - xbar))^2,
        k = 1 .. n)/n;
    alphaconf[1] := alphahat -
        TP(n - 2, 1 - 1/2*CL[1])*sqrt(1.0*sigma2hat/(n - 2));
    alphaconf[2] := alphahat +
        TP(n - 2, 1 - 1/2*CL[1])*sqrt(1.0*sigma2hat/(n - 2));
    betaconf[1] := betahat - TP(n - 2, 1 - 1/2*CL[2])*
        sqrt(1.0*n*sigma2hat/((n - 2)*betahat2));
    betaconf[2] := betahat + TP(n - 2, 1 - 1/2*CL[2])*
        sqrt(1.0*n*sigma2hat/((n - 2)*betahat2));
    sigma2conf[1] :=
        n*sigma2hat/ChisquareP(n - 2, 1 - 1/2*CL[3]);
    sigma2conf[2] := n*sigma2hat/ChisquareP(n - 2, 1/2*CL[3])
        ;
    c1 := 1 - CL[1];
    c2 := 1 - CL[2];
    c3 := 1 - CL[3];
    lprint(`----------------------------------------------\
        ---------------`);
    lprint(`                 Point      Confidence       C\
        onfidence`);
    lprint(`Parameter      Estimates       Level          \
        Interval`);
    lprint(`----------------------------------------------\
        ---------------`);
    printf("alpha     %13.4f   %10.2f       [%1.4f,%1.4f]\n",
        alphahat, c1, alphaconf[1], alphaconf[2]);
    lprint();
    printf("beta      %13.4f   %10.2f       [%1.4f,%1.4f]\n",
        betahat, c2, betaconf[1], betaconf[2]);
    lprint();
    printf("sigma^2   %13.4f   %10.2f       [%1.4f,%1.4f]\n",
        sigma2hat, c3, sigma2conf[1], sigma2conf[2]);
    lprint(`----------------------------------------------\
        ---------------`)

end:


RegBand :=proc()
local X, Y, N, A, Minx, Maxx, Miny, Maxy, XR, YR, EX, EY, i,
x, MX, alpha, beta, Sig2, SSQ, T, Option, ConfLev, Term1,
Term2, upper, lower, p1, p2, p3, p4, L, Ymin, Ymax;
option `August 21, 1993 -- Zaven A. Karian`;
description
"Produce scatter plot of X-Y pairs with the reg. line";
    if nargs = 3 and not type(args[1], listlist) then ERROR(
        "When one argument is used, it must be a list of l\
        ists")
    fi;
    if nargs = 4 and
    not (type(args[1], list) and type(args[2], list)) then
        ERROR(
        "When two arguments are used, both must be lists")
    fi;
    if nargs = 4 and nops(args[1]) <> nops(args[2]) then
        ERROR("The two lists must have the same length")
    fi;
    if nargs = 3 then
        X := [seq(args[1][i][1], i = 1 .. nops(args[1]))];
        Y := [seq(args[1][i][2], i = 1 .. nops(args[1]))];
        ConfLev := args[2];
        Option := args[3]
    else
        X := args[1];
        Y := args[2];
        ConfLev := args[3];
        Option := args[4]
    fi;
    N := nops(X);
    A := array(1 .. 2*N);
    MX := Mean(X);
    for i to N do A[2*i - 1] := X[i]; A[2*i] := Y[i] od;
    Minx := Min(X);
    Maxx := Max(X);
    Miny := Min(Y);
    Maxy := Max(Y);
    EX := 1/20*Maxx - 1/20*Minx;
    EY := 1/20*Maxy - 1/20*Miny;
    XR := Minx - EX .. Maxx + EX;
    L := LinReg(X, Y, x);
    beta := op(1, op(2, L));
    alpha := op(1, L) + beta*MX;
    Sig2 := 0;
    for i to N do
        Sig2 := Sig2 + (Y[i] - alpha - beta*(X[i] - MX))^2
    od;
    Sig2 := Sig2/N;
    SSQ := (N - 1)*Variance(X);
    if N < 33 then T := TP(N - 2, .01*ConfLev)
    else T := NormalP(0, 1, .01*ConfLev)
    fi;
    Term1 := alpha + beta*(x - MX);
    if Option = "C" then Term2 :=
        T*sqrt(N*Sig2*(1/N + (x - MX)^2/SSQ)/(N - 2))
    else Term2 :=
        T*sqrt(N*Sig2*(1 + 1/N + (x - MX)^2/SSQ)/(N - 2))
    fi;
    upper := Term1 + Term2;
    lower := Term1 - Term2;
    Ymin := min(Miny - EY, subs(x = Minx, lower),
        subs(x = Maxx, lower));
    Ymax := max(Maxy + EY, subs(x = Minx, upper),
        subs(x = Maxx, upper));
    YR := Ymin .. Ymax;
    p1 := plot([seq([A[2*i - 1], A[2*i]], i = 1 .. N)], XR,
        YR, style = POINT, color = black);
    p2 := plot(L, x = XR, YR, color = red);
    p3 := plot(upper, x = XR, YR, color = blue);
    p4 := plot(lower, x = XR, YR, color = blue);
    with(plots, display);
    display({p1, p2, p3, p4})

end:


Residuals := proc()
local X, Y, N, Minx, Maxx, XR, YR, p2, P, EX, EY, i, x, L,
alpha, beta, Res, p, f, m;
option `August 21, 1993 -- Zaven A. Karian`;
description
"Produce scatter plot of X-Y pairs with the reg. line";
    if nargs = 2 and not type(args[1], listlist) then ERROR(
        "When two arguments are used, the first two must b\
        e a list of lists")
    fi;
    if nargs = 3 and
    not (type(args[1], list) and type(args[2], list)) then
        ERROR("When three arguments are used, the first tw\
        o must be lists")
    fi;
    if nargs = 3 and nops(args[1]) <> nops(args[2]) then
        ERROR("The two lists must have the same length")
    fi;
    if nargs = 2 then
        X := [seq(args[1][i][1], i = 1 .. nops(args[1]))];
        Y := [seq(args[1][i][2], i = 1 .. nops(args[1]))];
        m := args[2]
    else X := args[1]; Y := args[2]; m := args[3]
    fi;
    N := nops(X);
    if m = 1 then L := LinReg(X, Y, x)
    else L := PolyReg(X, Y, m, x)
    fi;
    f := student[makeproc](L, x);
    Res := [seq(Y[i] - f(X[i]), i = 1 .. N)];
    Minx := Min(X);
    Maxx := Max(X);
    EX := 1/20*Maxx - 1/20*Minx;
    XR := Minx - EX .. Maxx + EX;
    EY := .05000000000*Max(Res) - .05000000000*Min(Res);
    YR := Min(Res) - EY .. Max(Res) + EY;
    P := {};
    for i to N do P := P union {[X[i], Res[i]]} od;
    p2 := plot(P, XR, YR, style = POINT, color = red);
    print(p2);
    lprint(`             Independent`);
    lprint(`              Variable      Residual`);
    lprint(`             -----------------------`);
    for i to N do
        printf(" %20.4f %13.4f", X[i], Res[i]); lprint()
    od

end:


ScatPlot := proc()
local X, Y, N, A, Minx, Maxx, Miny, Maxy, XR, YR, EX, EY, i;
option `January 27, 1994 -- Zaven A. Karian`;
description "Produce a scatter plot of X-Y pairs";
    if nargs = 1 and not type(args[1], listlist) then ERROR(
        "When one argument is used it must be a list of li\
        sts")
    fi;
    if nargs = 2 and
    not (type(args[1], list) and type(args[2], list)) then
        ERROR(
        "When two arguments are used, both must be lists")
    fi;
    if nargs = 2 and nops(args[1]) <> nops(args[2]) then
        ERROR("The two lists must have the same length")
    fi;
    if nargs = 1 then
        X := [seq(args[1][i][1], i = 1 .. nops(args[1]))];
        Y := [seq(args[1][i][2], i = 1 .. nops(args[1]))]
    else X := args[1]; Y := args[2]
    fi;
    N := nops(X);
    A := array(1 .. 2*N);
    for i to N do A[2*i - 1] := X[i]; A[2*i] := Y[i] od;
    Minx := Min(X);
    Maxx := Max(X);
    Miny := Min(Y);
    Maxy := Max(Y);
    EX := 1/20*Maxx - 1/20*Minx;
    EY := 1/20*Maxy - 1/20*Miny;
    XR := Minx - EX .. Maxx + EX;
    YR := Miny - EY .. Maxy + EY;
    plot([seq([A[2*i - 1], A[2*i]], i = 1 .. N)], XR, YR,
        style = POINT)

end:


ScatPlotLine, proc()
local X, Y, N, A, Minx, Maxx, Miny, Maxy, XR, YR, EX, EY, i,
x, p1, p2, L;
option `August 21, 1993 -- Zaven A. Karian`;
description
"Produce scatter plot of X-Y pairs with the reg. line";
    if nargs = 1 and not type(args[1], listlist) then ERROR(
        "When one argument is used, it must be a list of l\
        ists")
    fi;
    if nargs = 2 and
    not (type(args[1], list) and type(args[2], list)) then
        ERROR(
        "When two arguments are used, both must be lists")
    fi;
    if nargs = 2 and nops(args[1]) <> nops(args[2]) then
        ERROR("The two lists must have the same length")
    fi;
    if nargs = 1 then
        X := [seq(args[1][i][1], i = 1 .. nops(args[1]))];
        Y := [seq(args[1][i][2], i = 1 .. nops(args[1]))]
    else X := args[1]; Y := args[2]
    fi;
    N := nops(X);
    A := array(1 .. 2*N);
    for i to N do A[2*i - 1] := X[i]; A[2*i] := Y[i] od;
    Minx := Min(X);
    Maxx := Max(X);
    Miny := Min(Y);
    Maxy := Max(Y);
    EX := 1/20*Maxx - 1/20*Minx;
    EY := 1/20*Maxy - 1/20*Miny;
    XR := Minx - EX .. Maxx + EX;
    YR := Miny - EY .. Maxy + EY;
    L := LinReg(X, Y, x);
    p1 := plot([seq([A[2*i - 1], A[2*i]], i = 1 .. N)], XR,
        YR, style = POINT);
    p2 := plot(L, x = XR, YR);
    with(plots, display);
    display({p1, p2})

end:


BoxWhisker := proc()
local N, i, A, D, m, M, Minn, Maxx, P25, P50, P75, eh;
option `August 19, 1993 -- Zaven A. Karian`;
description "Box-and whisker display";
    N := nargs;
    Minn := Min(args[1]);
    Maxx := Max(args[1]);
    for i to N do
        if not type(args[i], list) then
            ERROR("All arguments must be lists")
        fi
    od;
    A := {};
    for i to N do
        D := args[i];
        m := Min(D);
        M := Max(D);
        if m < Minn then Minn := m fi;
        if Maxx < M then Maxx := M fi;
        P25 := Percentile(D, .25);
        P50 := Percentile(D, .50);
        P75 := Percentile(D, .75);
        A := A union {[[P75, i - .4], [P75, i + .4]],
            [[P25, i + .4], [P75, i + .4]],
            [[P25, i - .4], [P25, i + .4]],
            [[P50, i - .4], [P50, i + .4]],
            [[m, i], [P25, i]], [[P75, i], [M, i]],
            [[P25, i - .4], [P75, i - .4]]}
    od;
    eh := 1/20*Maxx - 1/20*Minn;
    plot(A, Minn - eh .. Maxx + eh, 0 .. N + 1, color = blue)

end:


Histogram := proc()
local YYY, YY, Y, n, l, u, v, i, j, k, p, t, U, V, count,
truemin, truemax, nint;
option `July 7, 1993 -- Zaven A. Karian`;
description "Graphic display of histogram";
    YY := sort(evalf(args[1]));
    n := nops(YY);
    if nargs = 1 then
        truemin := YY[1]; truemax := YY[n]; nint := 8
    else
        truemin := lhs(args[2]);
        truemax := rhs(args[2]);
        nint := args[3]
    fi;
    count := 0;
    for i to n do
        if YY[i] < truemin or truemax < YY[i] then
            count := count + 1
        fi
    od;
    if 0 < count then print("WARNING: There are", count,
        " data points out of plot range.")
    fi;
    j := 1;
    YYY := array(1 .. n);
    u := array(1 .. nint + 1);
    v := array(1 .. nint + 1);
    for i to n do
        if truemin <= YY[i] and YY[i] <= truemax then
            YYY[j] := YY[i]; j := j + 1
        fi
    od;
    n := j - 1;
    Y := [seq(YYY[k], k = 1 .. n)];
    p := array(1 .. 8*nint);
    l := (truemax - truemin)/nint;
    u[1] := truemin;
    v[1] := 0;
    for i from 2 to nint + 1 do
        u[i] := u[i - 1] + l; v[i] := 0
    od;
    if u[nint + 1] < Y[n] then u[nint + 1] := Y[n] fi;
    i := 2;
    for j to n do
        if Y[j] <= u[i] then v[i - 1] := v[i - 1] + 1
        else i := i + 1; j := j - 1
        fi
    od;
    t := 1;
    for i to nint do
        U[t] := u[i];
        V[t] := 0;
        t := t + 1;
        U[t] := U[t - 1];
        V[t] := v[i];
        t := t + 1;
        U[t] := U[t - 1] + l;
        V[t] := v[i];
        t := t + 1;
        U[t] := U[t - 1];
        V[t] := 0;
        t := t + 1
    od;
    for i to t - 1 do
        p[2*i - 1] := U[i]; p[2*i] := V[i]/(l*n)
    od;
    plot([seq([p[2*i - 1], p[2*i]], i = 1 .. t - 1)],
        color = red, linestyle = 3)

end:


Ogive, proc(L::list, R::range, n::posint)
local Min, Max, freq, A, X, LEP, REP, j, inc, i, N;
option `July 5, 1993 -- Zaven A. Karian`;
description "Graphic display of ogive for L";
    Min := lhs(R);
    Max := rhs(R);
    freq := array(1 .. n);
    A := array(1 .. 2*n + 2);
    N := nops(L);
    X := sort(L);
    LEP := Min;
    inc := (Max - Min)/n;
    REP := LEP + inc;
    for j to n do freq[j] := 0 od;
    for j to n do
        for i to N do
            if X[i] < REP then freq[j] := freq[j] + 1 fi
        od;
        LEP := REP;
        REP := LEP + inc
    od;
    A[1] := Min;
    A[2] := 0;
    for j to n do
        A[2*j + 1] := A[2*j - 1] + inc;
        A[2*j + 2] := freq[j]/N
    od;
    [seq([A[2*i - 1], A[2*i]], i = 1 .. n + 1)];
    plot(%, color = blue, linestyle = 3)

end:


PlotDiscCDF := proc(expr::algebraic, R::range)
local A, MIN, p, i, var, Min, Max, newp;
option `August 8, 1993 -- Zaven A. Karian`;
description "Graphic display of distribution expr";
    Min := lhs(R);
    Max := rhs(R);
    A := {};
    MIN := Min;
    p := 0;
    var := op(select(type, indets(expr), name));
    for i to Max - Min + 1 do
        newp := eval(subs(var = MIN, expr));
        if p < 1 and 0 <= newp then
            p := p + eval(subs(var = MIN, expr))
        fi;
        A := A union {[[MIN, p], [MIN + 1, p]]};
        MIN := MIN + 1
    od;
    plot(A, color = blue)

end:


PlotEmpCDF := proc(L::list, R::range)
local i, A, N, X;
option `August 8, 1993 -- Zaven A. Karian`;
description "Graphic display of empirical distribution L";
    N := nops(L);
    X := sort(L);
    A := {[[min(lhs(R), X[1]), 0], [X[1], 0]]};
    for i to N - 1 do
        if X[i] < X[i + 1] then
            A := A union {[[X[i], i/N], [X[i + 1], i/N]]}
        fi
    od;
    A := A union {[[X[N], 1], [max(rhs(R), X[N]), 1]]};
    plot(A, color = red, linestyle = 3)

end:


PlotEmpPDF := proc()
local YYY, YY, Y, n, l, u, v, i, j, k, p, t, U, V, count,
truemin, truemax, nint;
option `June 1, 1996 -- Zaven A. Karian`;
description "Graphic display of histogram";
    YY := sort(evalf(args[1]));
    n := nops(YY);
    if nargs = 1 then
        truemin := YY[1]; truemax := YY[n]; nint := 8
    else
        truemin := lhs(args[2]);
        truemax := rhs(args[2]);
        nint := args[3]
    fi;
    count := 0;
    for i to n do
        if YY[i] < truemin or truemax < YY[i] then
            count := count + 1
        fi
    od;
    if 0 < count then print("WARNING: There are", count,
        " data points out of plot range.")
    fi;
    j := 1;
    YYY := array(1 .. n);
    u := array(1 .. nint + 1);
    v := array(1 .. nint + 1);
    for i to n do
        if truemin <= YY[i] and YY[i] <= truemax then
            YYY[j] := YY[i]; j := j + 1
        fi
    od;
    n := j - 1;
    Y := [seq(YYY[k], k = 1 .. n)];
    p := array(1 .. 2*nint + 4);
    l := (truemax - truemin)/nint;
    u[1] := truemin;
    v[1] := 0;
    for i from 2 to nint + 1 do
        u[i] := u[i - 1] + l; v[i] := 0
    od;
    if u[nint + 1] < Y[n] then u[nint + 1] := Y[n] fi;
    i := 2;
    for j to n do
        if Y[j] <= u[i] then v[i - 1] := v[i - 1] + 1
        else i := i + 1; j := j - 1
        fi
    od;
    U[1] := truemin - 1/2*l;
    V[1] := 0;
    t := 2;
    for i to nint do
        U[t] := u[i] + 1/2*l; V[t] := v[i]; t := t + 1
    od;
    U[nint + 2] := truemax + 1/2*l;
    V[nint + 2] := 0;
    for i to nint + 2 do
        p[2*i - 1] := U[i]; p[2*i] := V[i]/(l*n)
    od;
    [seq([p[2*i - 1], p[2*i]], i = 1 .. nint + 2)];
    plot(%, color = red, linestyle = 3)

end:


PlotRunningAverage := proc(L::list)
local N, SUM, A, i;
option `August 23, 1993 -- Zaven A. Karian`;
description "Plotting the running average of a list";
    N := nops(L);
    SUM := L[1];
    A := array(1 .. 2*N);
    A[1] := 1;
    A[2] := L[1];
    for i from 2 to N do
        SUM := SUM + L[i]; A[2*i - 1] := i; A[2*i] := SUM/i
    od;
    plot([seq([A[2*i - 1], A[2*i]], i = 1 .. N)], color = red)

end:


ProbHist := proc()
local n, L, lp, ap, A, MIN, dim, i, p, var, Min, Max;
option `January 26, 1994 -- Zaven A. Karian`;
description "Graphic display of probability histogram";
    if nargs = 1 and not type(args[1], list) then ERROR("I\
        f one argument is used, it must be a probability l\
        ist")
    fi;
    if nargs = 2 and not type(args[1], algebraic) then ERROR(
        "If two arguments are used, the first must be an e\
        xpression")
    fi;
    if nargs = 2 and not type(args[2], range) then ERROR("\
        If two arguments are used, the second must be a ra\
        nge")
    fi;
    if nargs = 1 then
        n := nops(args[1]);
        L := copy(args[1]);
        A := array(1 .. 4*n);
        lp := 1;
        ap := 1;
        for i to 1/2*n do
            A[ap] := L[lp] - .5;
            A[ap + 1] := 0;
            A[ap + 2] := A[ap];
            A[ap + 3] := L[lp + 1];
            A[ap + 4] := L[lp] + .5;
            A[ap + 5] := L[lp + 1];
            A[ap + 6] := A[ap + 4];
            A[ap + 7] := 0;
            lp := lp + 2;
            ap := ap + 8
        od;
        plot([seq([A[2*i - 1], A[2*i]], i = 1 .. 2*n)],
            color = blue)
    else
        Min := lhs(args[2]);
        Max := rhs(args[2]);
        dim := 8*trunc(Max - Min + 1/1000000000) + 8;
        A := array(1 .. dim);
        MIN := Min;
        var := op(select(type, indets(args[1]), name));
        for i to Max - Min + 1 do
            p := evalf(subs(var = MIN, args[1]));
            A[8*i - 7] := MIN - .5;
            A[8*i - 6] := 0;
            A[8*i - 5] := MIN - .5;
            A[8*i - 4] := p;
            A[8*i - 3] := MIN + .5;
            A[8*i - 2] := p;
            A[8*i - 1] := MIN + .5;
            A[8*i] := 0;
            MIN := MIN + 1
        od;
        [seq([A[2*i - 1], A[2*i]], i = 1 .. 4*Max - 4*Min + 4)
            ];
        plot(%, color = blue)
    fi

end:


QQ := proc()
local X, Y, N, M, XX, YY, MIN, P, NY, NEWY, NEWYL, i, r, f;
option `August 12, 1993 -- Zaven A. Karian`;
description "Produce the Q-Q plot of the lists X and Y";
    if nargs = 1 and not type(args[1], listlist) then ERROR(
        "When one argument is used, it must be a list of l\
        ists")
    fi;
    if nargs = 2 and
    not (type(args[1], list) and type(args[2], list)) then
        ERROR(
        "When two arguments are used, both must be lists")
    fi;
    if nargs = 1 then
        X := [seq(args[1][i][1], i = 1 .. nops(args[1]))];
        Y := [seq(args[1][i][2], i = 1 .. nops(args[1]))]
    else X := args[1]; Y := args[2]
    fi;
    N := nops(X);
    M := nops(Y);
    if N <= M then XX := sort(X); YY := sort(Y)
    else XX := sort(Y); YY := sort(X)
    fi;
    MIN := min(N, M);
    NEWY := array(1 .. MIN);
    P := [seq(i/(MIN + 1), i = 1 .. MIN)];
    NY := nops(YY);
    for i to MIN do
        r := trunc(P[i]*(NY + 1));
        f := P[i]*(NY + 1) - r;
        if f = 0 then NEWY[i] := YY[r]
        else NEWY[i] := YY[r] + f*(YY[r + 1] - YY[r])
        fi
    od;
    NEWYL := [seq(NEWY[i], i = 1 .. MIN)];
    ScatPlot(XX, NEWYL)

end:


QQFit := proc(L::list, expr::algebraic, R::range)
local
var, DataList, DistList, NQuantiles, Range, PlotRange, i;
    var := op(select(type, indets(expr), name));
    if nops([var]) = 0 or 1 < nops([var]) then ERROR("The \
        second argument must have exactly one varaible in \
        it")
    fi;
    if nargs = 2 then
        Range := -infinity .. infinity;
        PlotRange := Min(L) .. Max(L)
    elif nargs = 3 then Range := R; PlotRange := R
    else ERROR(
        "This procedure requires either 2 or 3 arguments.")
    fi;
    NQuantiles := nops(L);
    DataList := sort(L);
    DistList := [seq(
        fsolve(i/(NQuantiles + 1) = expr, var, Range),
        i = 1 .. NQuantiles)];
    plot({
        plot(x, x = PlotRange), ScatPlot(DistList, DataList)})

end:


TimePlot, proc()
local N, B, n, i, j, t, S, A;
option `August 19, 1993 -- Zaven A. Karian`;
description "Time series plot of a list";
    N := nargs;
    n := nops(args[1]);
    S := {};
    for i to N do
        if not type(args[i], list) then
            ERROR("All arguments must be lists")
        fi;
        if n < nops(args[i]) then n := nops(args[i]) fi
    od;
    B := array(1 .. 2*n);
    for i to N do
        A := args[i];
        for j to nops(A) do B[2*j - 1] := j; B[2*j] := A[j]
        od;
        S := S union
            {[seq([B[2*t - 1], B[2*t]], t = 1 .. nops(A))]}
    od;
    plot(S)

end:


XbarChart := proc(X::list, SampleSize::posint, sbar::numeric)
local N, XX, S, Sbar, n, Xbarbar, A3, i, offset, UCL, LCL,
pts, P1, P2;
option `January 24, 1994 -- Zaven A. Karian`;
description "Produce the control chart for the mean ";
    N := nops(X);
    if type(X, listlist) then
        XX := [seq(Mean(X[i]), i = 1 .. N)];
        S := [seq(StDev(X[i]), i = 1 .. N)];
        Sbar := Mean(S);
        n := nops(X[1])
    else XX := copy(X); Sbar := sbar; n := SampleSize
    fi;
    Xbarbar := Mean(XX);
    A3 := 3*sqrt(n - 1)*GAMMA(1/2*n - 1/2)/(
        sqrt(n)*sqrt(2)*GAMMA(1/2*n));
    offset := A3*Sbar;
    UCL := [[0, Xbarbar + offset], [N + 1, Xbarbar + offset]]
        ;
    LCL := [[0, Xbarbar - offset], [N + 1, Xbarbar - offset]]
        ;
    pts := [seq([i, XX[i]], i = 1 .. N)];
    P1 := plot({UCL, LCL}, color = red);
    P2 := plot(pts, color = blue);
    with(plots, display);
    display({P1, P2})

end:


BetaP := proc(b3::algebraic, b4::algebraic, p::algebraic)
local eqn, C, x, y;
option `August 12, 1993 -- Zaven A. Karian`;
description "Beta percentile";
    if not (type(b3, numeric) and type(b4, numeric) and
    type(p, numeric)) then RETURN('procname(args)')
    fi;
    if b3 <= -1 then ERROR("First argument must be >= 1") fi;
    if b4 <= -1 then ERROR("Second argument must be >= 1") fi
    ;
    if p < 0 or 1 < p then
        ERROR("Third argument must be between 0 and 1")
    fi;
    C := GAMMA(b3 + b4 + 2)/(GAMMA(b3 + 1)*GAMMA(b4 + 1));
    eqn := p = int(C*x^b3*(1 - x)^b4, x = 0 .. y);
    fsolve(eqn, y, y = 0 .. 1)

end:


ChisquareP := proc(r::algebraic, p::algebraic)
local eqn, a, x, y;
option `August 12, 1993 -- Zaven A. Karian`;
description "Chi square percentile";
    if not (type(r, numeric) and type(p, numeric)) then
        RETURN('procname(args)')
    fi;
    if p < 0 or 1 < p then
        ERROR("Second argument must be between 0 and 1")
    fi;
    a := 1/(GAMMA(1/2*r)*2^(1/2*r));
    eqn := p = int(a*x^(1/2*r - 1)*exp(- 1/2*x), x = 0 .. y);
    fsolve(eqn, y, y = 0 .. infinity)

end:


ExponentialP := proc(t::algebraic, p::algebraic)
local eqn, x, y;
option `August 12, 1993 -- Zaven A. Karian`;
description "Exponential percentile";
    if not (type(t, numeric) and type(p, numeric)) then
        RETURN('procname(args)')
    fi;
    if p < 0 or 1 < p then
        ERROR("Second argument must be between 0 and 1")
    fi;
    if t <= 0 then ERROR("First argument must be positive")
    fi;
    eqn := p = int(exp(- x/t)/t, x = 0 .. y);
    fsolve(eqn, y, y = 0 .. infinity)

end:


FP := proc(df1::algebraic, df2::algebraic, p::algebraic)
local eqn, num, den, x, y;
option `August 12, 1993 -- Zaven A. Karian`;
description "F percentile";
    if not (type(df1, numeric) and type(df2, numeric) and
    type(p, numeric)) then RETURN('procname(args)')
    fi;
    if p < 0 or 1 < p then
        ERROR("Third argument must be between 0 and 1")
    fi;
    if df1 < 0 or df2 < 0 then
        ERROR("First two arguments must be positive")
    fi;
    num := x^(1/2*df1 - 1)*GAMMA(1/2*df1 + 1/2*df2)*
        (df1/df2)^(1/2*df1);
    den := GAMMA(1/2*df1)*GAMMA(1/2*df2)*
        (1 + df1*x/df2)^(1/2*df1 + 1/2*df2);
    eqn := p = int(num/den, x = 0 .. y);
    fsolve(eqn, y, y = 0 .. infinity)

end:


GBP := proc(b1::algebraic, b2::algebraic, b3::algebraic,b4::algebraic, p::algebraic)
local eqn, C, CC, x, y;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generalized beta percentile";
    if not (type(b1, numeric) and type(b2, numeric) and
    type(b3, numeric) and type(b4, numeric) and
    type(p, numeric)) then RETURN('procname(args)')
    fi;
    if b3 <= -1 then ERROR("Third argument must be >= 1") fi;
    if b4 <= -1 then ERROR("Fourth argument must be >= 1") fi
    ;
    if p < 0 or 1 < p then
        ERROR("Fifth argument must be between 0 and 1")
    fi;
    C := GAMMA(b3 + b4 + 2)/(GAMMA(b3 + 1)*GAMMA(b4 + 1));
    CC := (1/b2)^(b3 + b4 + 1);
    eqn := p =
        int(C*CC*(x - b1)^b3*(b1 + b2 - x)^b4, x = b1 .. y);
    fsolve(eqn, y)

end:


GammaP := proc(a::algebraic, t::algebraic, p::algebraic)
local eqn, b, x, y;
option `August 12, 1993 -- Zaven A. Karian`;
description "Gamma percentile";
    if not (type(a, numeric) and type(t, numeric) and
    type(p, numeric)) then RETURN('procname(args)')
    fi;
    if p < 0 or 1 < p then
        ERROR("Third argument must be between 0 and 1")
    fi;
    if a <= 0 then ERROR("First argument must be positive")
    fi;
    if t <= 0 then ERROR("Second argument must be positive")
    fi;
    b := 1/(GAMMA(a)*t^a);
    eqn := p = int(b*x^(a - 1)*exp(- x/t), x = 0 .. y);
    fsolve(eqn, y, y = 0 .. infinity)

end:


NormalP := proc(mu::algebraic, var::algebraic, p::algebraic)
local eqn, t, X, x, guess;
option `August 12, 1993 -- Zaven A. Karian`;
description "Normal percentile";
    if not (type(mu, numeric) and type(var, numeric) and
    type(p, numeric)) then RETURN('procname(args)')
    fi;
    if p < 0 or 1 < p then
        ERROR("Third argument must be between 0 and 1")
    fi;
    if var <= 0 then
        ERROR("Second argument must be positive")
    fi;
    t := x/sqrt(2);
    guess := 5.063291139*p^.135 - 5.063291139*(1 - p)^.135;
    eqn := p = evalf(.5 + .5*erf(t));
    X := fsolve(eqn, x, guess - .5 .. guess + .5);
    sqrt(var)*X + mu

end:


TP := proc(df::algebraic, p::algebraic)
local eqn, C1, C2, C3, x, y;
option `August 12, 1993 -- Zaven A. Karian`;
description "Students t percentile";
    if not (type(df, numeric) and type(p, numeric)) then
        RETURN('procname(args)')
    fi;
    if p < 0 or 1 < p then
        ERROR("Second argument must be between 0 and 1")
    fi;
    C1 := GAMMA(1/2*df + 1/2);
    C2 := GAMMA(1/2*df)*sqrt(df*Pi);
    C3 := (1 + x^2/df)^(1/2*df + 1/2);
    eqn := p = int(C1/(C2*C3), x = -infinity .. y);
    fsolve(eqn, y)

end:


UniformP := proc(R::range, p::algebraic)
local a, b, x;
option `August 12, 1993 -- Zaven A. Karian`;
description "uniform percentile";
    a := lhs(R);
    b := rhs(R);
    if not (type(a, numeric) and type(b, numeric) and
    type(p, numeric)) then RETURN('procname(args)')
    fi;
    if p < 0 or 1 < p then
        ERROR("Second argument must be between 0 and 1")
    fi;
    if b < a then ERROR("First argument must be a range") fi;
    fsolve((x - a)/(b - a) = p, x)

end:


BernoulliCDF := proc(p::algebraic, y::algebraic)
local i;
option `August 8, 1993 -- Zaven A. Karian`;
description "Bernoulli CDF";
    if type(y, numeric) and y < 0 then 0
    elif type(y, numeric) and 1 < y then 1
    else evalf(sum(p^i*(1 - p)^(1 - i), i = 0 .. y))
    fi

end:


BetaCDF := proc(b3::algebraic, b4::algebraic, y::algebraic)
local C, x, A;
option `August 8, 1993 -- Zaven A. Karian`;
description "Beta CDF";
    if type(y, numeric) and y < 0 then A := 0
    elif type(y, numeric) and 1 < y then A := 1
    else
        C := GAMMA(b3 + b4 + 2)/(GAMMA(b3 + 1)*GAMMA(b4 + 1))
            ;
        A := evalf(int(C*x^b3*(1 - x)^b4, x = 0 .. y))
    fi;
    A

end:


BinomialCDF := proc(N::algebraic, p::algebraic, y::algebraic)
local i;
option `August 8, 1993 -- Zaven A. Karian`;
description "Binomial CDF";
    if type(y, numeric) and y < 0 then 0
    elif type(y, numeric) and type(N, numeric) and N < y then
        1
    else evalf(
        sum(binomial(N, i)*p^i*(1 - p)^(N - i), i = 0 .. y))
    fi

end:


ChisquareCDF := proc(r::algebraic, y::algebraic)
local a, x;
option `August 8, 1993 -- Zaven A. Karian`;
description "Chi square CDF";
    if type(y, numeric) and y < 0 then 0
    else
        a := 1/(GAMMA(1/2*r)*2^(1/2*r));
        evalf(int(a*x^(1/2*r - 1)*exp(- 1/2*x), x = 0 .. y))
    fi

end:


ExponentialCDF := proc(t::algebraic, y::algebraic)
local x;
option `August 8, 1993 -- Zaven A. Karian`;
description "Exponential CDF";
    if type(y, numeric) and y < 0 then 0
    else evalf(int(exp(- x/t)/t, x = 0 .. y))
    fi

end:


FCDF := proc(df1::algebraic, df2::algebraic, y::algebraic)
local num, den, x, A;
option `August 8, 1993 -- Zaven A. Karian`;
description "F distribution CDF";
    if type(y, numeric) and y < 0 then A := 0
    else
        num := x^(1/2*df1 - 1)*GAMMA(1/2*df1 + 1/2*df2)*
            (df1/df2)^(1/2*df1);
        den := GAMMA(1/2*df1)*GAMMA(1/2*df2)*
            (1 + df1*x/df2)^(1/2*df1 + 1/2*df2);
        A := evalf(int(num/den, x = 0 .. y))
    fi;
    A

end:


GBCDF := proc(b1::algebraic, b2::algebraic,b3::algebraic, b4::algebraic, y::algebraic)
local C, CC, x, A;
option `August 8, 1993 -- Zaven A. Karian`;
description "Generalized beta CDF";
    if type(b1, numeric) and type(y, numeric) and y < b1 then
        A := 0
    elif type(b2, numeric) and type(y, numeric) and b2 < y
    then A := 1
    else
        C := GAMMA(b3 + b4 + 2)/(GAMMA(b3 + 1)*GAMMA(b4 + 1))
            ;
        CC := (1/b2)^(b3 + b4 + 1);
        A := evalf(int(C*CC*(x - b1)^b3*(b1 + b2 - x)^b4,
            x = b1 .. y))
    fi;
    A

end:


GammaCDF := proc(a::algebraic, t::algebraic, y::algebraic)
local b, x, A;
option `August 8, 1993 -- Zaven A. Karian`;
description "Gamma CDF";
    if type(y, numeric) and y < 0 then A := 0
    else
        b := 1/(GAMMA(a)*t^a);
        A := evalf(int(b*x^(a - 1)*exp(- x/t), x = 0 .. y))
    fi;
    A

end:


GeometricCDF := proc(p::algebraic, y::algebraic)
local i;
option `August 8, 1993 -- Zaven A. Karian`;
description "Geometric CDF";
    if type(y, numeric) and y < 1 then 0
    else evalf(sum(p*(1 - p)^(i - 1), i = 1 .. y))
    fi

end:


NegBinomialCDF := proc(r::algebraic, p::algebraic, y::algebraic)
local i;
option `August 8, 1993 -- Zaven A. Karian`;
description "Negative binomial CDF";
    if type(y, numeric) and type(r, numeric) and y < r then 0
    else evalf(sum(
        binomial(i - 1, r - 1)*p^r*(1 - p)^(i - r),
        i = r .. y))
    fi

end:


NormalCDF := proc(mu::algebraic, var::algebraic, y::algebraic)
local t;
option `July 7, 1993 -- Zaven A. Karian`;
description "Normal CDF";
    t := (y - mu)/sqrt(2*var); evalf(.5 + .5*erf(t))

end:


PoissonCDF := proc(l::algebraic, y::algebraic)
local i;
option `Augusr 8, 1993 -- Zaven A. Karian`;
description "Poisson CDF";
    if type(y, numeric) and y < 0 then 0
    else evalf(sum(l^i*exp(-l)/i!, i = 0 .. y))
    fi

end:


TCDF := proc(df::algebraic, y::algebraic)
local C1, C2, C3, x;
option `July 7, 1993 -- Zaven A. Karian`;
description "Students t CDF";
    C1 := GAMMA(1/2*df + 1/2);
    C2 := GAMMA(1/2*df)*sqrt(df*Pi);
    C3 := (1 + x^2/df)^(1/2*df + 1/2);
    evalf(int(C1/(C2*C3), x = -infinity .. y))

end:


UniformCDF := proc(R::range, y::algebraic)
local l, r;
option `August 8, 1993 -- Zaven A. Karian`;
description "Uniform CDF";
    l := lhs(R);
    r := rhs(R);
    if type(l, numeric) and type(y, numeric) and y < l then 0
    elif type(r, numeric) and type(y, numeric) and r < y then
        1
    else evalf((y - l)/(r - l))
    fi

end:


BernoulliPDF := proc(p::algebraic, x::algebraic)
option `August 9, 1993 -- Zaven A. Karian`;
description "Bernoulli density";
    if type(p, numeric) and (p <= 0 or 1 <= p) then ERROR(
        "First argument must be a variable or between 0 an\
        d 1")
    fi;
    if type(x, numeric) and
    (x < 0 or 1 < x or not type(x, integer)) then RETURN(0)
    fi;
    p^x*(1 - p)^(1 - x)

end:


BetaPDF := proc(b3::algebraic, b4::algebraic, x::algebraic)
local C;
option `August 11, 1993 -- Zaven A. Karian`;
description "Beta density";
    if type(b3, numeric) and b3 <= 0 then ERROR(
        "First argument must be a variable or positive")
    fi;
    if type(b4, numeric) and b4 <= 0 then ERROR(
        "Second argument must be a variable or positive")
    fi;
    if type(x, numeric) and (x <= 0 or 1 <= x) then RETURN(0)
    fi;
    C := GAMMA(b3 + b4 + 2)/(GAMMA(b3 + 1)*GAMMA(b4 + 1));
    C*x^b3*(1 - x)^b4

end:


BinomialPDF := proc(N::algebraic, p::algebraic, x::algebraic)
option `July 7, 1993 -- Zaven A. Karian`;
description "Binomial density";
     if type(N, numeric) and not type(N, integer) then ERROR(
        "First argument must be a variable or a positive i\
         nteger")
     fi;
     if type(N, numeric) and N <= 0 then ERROR("First argum\
         ent must be a variable or a  positive integer")
     fi;
     if type(p, numeric) and (p < 0 or 1 < p) then ERROR("S\
         econd argument must be an expression or between 0 \
         and 1")
     fi;
     if type(x, numeric) and type(N, numeric) and
     (x < 0 or N < x) then RETURN(0)
     fi;
     if type(x, numeric) and (not type(x, integer) or x < 0)
     then RETURN(0)
     fi;
    
      binomial(N, x)*p^x*(1 - p)^(N - x)
end:


BivariateNormalPDF:= proc(mux::algebraic,varx::algebraic, x::algebraic, 
    muy::algebraic,vary::algebraic, y::algebraic, rho::algebraic)
local A, B, C, Q;
option `August 11, 1993 -- Zaven A. Karian`;
description "Bivariate normal density";
    if type(varx, numeric) and varx <= 0 then ERROR(
        "Second argument must be a variable or positive")
    fi;
    if type(vary, numeric) and vary <= 0 then ERROR(
        "Fifth argument must be a variable or positive")
    fi;
    if type(rho, numeric) and (rho < -1 or 1 < rho) then
        ERROR("Seventh argument must be a variable or betw\
        een -1 and +1")
    fi;
    A := 2*Pi*sqrt(varx*vary*(1 - rho^2));
    B := (x - mux)^2/varx + (y - muy)^2/vary;
    C := -2*rho*(x - mux)*(y - muy)/sqrt(varx*vary);
    Q := (B + C)/(1 - rho^2);
    exp(-Q)/A

end:


CauchyPDF := proc(x::algebraic)
option `July 16, 1993 -- Zaven A. Karian`;
description "Cauchy density";
    1/(Pi*(1 + x^2))

end:


ChisquarePDF := proc(r::algebraic, x::algebraic)
local a;
option `August 11, 1993 -- Zaven A. Karian`;
description "Chi square density";
    if type(r, numeric) and r <= 0 then ERROR(
        "First argument must be a variable or positive")
    fi;
    if type(x, numeric) and x <= 0 then RETURN(0) fi;
    a := 1/(GAMMA(1/2*r)*2^(1/2*r));
    a*x^(1/2*r - 1)*exp(- 1/2*x)

end:


DoubleExponentialPDF := proc(t::algebraic, sig::algebraic, x::algebraic)
local A;
option `August 11, 1993 -- Zaven A. Karian`;
description "Double exponential density";
    if type(sig, numeric) and sig <= 0 then ERROR(
        "Second argument must be a variable or positive")
    fi;
    A := - abs(x - t)/sig;
    1/2*exp(A)/sig

end:


ExponentialPDF := proc(t::algebraic, x::algebraic)
option `August 11, 1993 -- Zaven A. Karian`;
description "Exponentiall density";
    if type(t, numeric) and t <= 0 then ERROR(
        "First argument must be a variable or positive")
    fi;
    if type(x, numeric) and x < 0 then RETURN(0) fi;
    exp(- x/t)/t

end:


FPDF := proc(df1::algebraic, df2::algebraic, x::algebraic)
local num, den;
option `August 11, 1993 -- Zaven A. Karian`;
description "F density";
    if type(df1, numeric) and not type(df1, posint) then
        ERROR("First argument must be a variable or a posi\
        tive integer")
    fi;
    if type(df2, numeric) and not type(df2, posint) then
        ERROR("Second argument must be a variable or a pos\
        itive integer")
    fi;
    if type(x, numeric) and x <= 0 then RETURN(0) fi;
    num := x^(1/2*df1 - 1)*GAMMA(1/2*df1 + 1/2*df2)*
        (df1/df2)^(1/2*df1);
    den := GAMMA(1/2*df1)*GAMMA(1/2*df2)*
        (1 + df1*x/df2)^(1/2*df1 + 1/2*df2);
    num/den

end:


GBDPDF := proc(b1::algebraic, b2::algebraic,b3::algebraic, b4::algebraic, x::algebraic)
local C, CC;
option `August 11, 1993 -- Zaven A. Karian`;
description "Generalized beta density";
    if type(b3, numeric) and b3 <= -1 then
        ERROR("Third argument must be a variable or > -1")
    fi;
    if type(b4, numeric) and b4 <= -1 then
        ERROR("Fourth argument must be a variable or > -1")
    fi;
    C := GAMMA(b3 + b4 + 2)/(GAMMA(b3 + 1)*GAMMA(b4 + 1));
    CC := (1/b2)^(b3 + b4 + 1);
    C*CC*(x - b1)^b3*(b1 + b2 - x)^b4

end:


GammaPDF := proc(a::algebraic, t::algebraic, x::algebraic)
local b;
option `August 11, 1993 -- Zaven A. Karian`;
description "Gamma density";
    if type(a, numeric) and a <= 0 then ERROR(
        "First argument must be a variable or positive")
    fi;
    if type(t, numeric) and t <= 0 then ERROR(
        "Second argument must be a variable or positive")
    fi;
    if type(x, numeric) and x < 0 then RETURN(0) fi;
    b := 1/(GAMMA(a)*t^a);
    b*x^(a - 1)*exp(- x/t)

end:


GeometricPDF := proc(p::algebraic, x::algebraic)
option `August 11, 1993 -- Zaven A. Karian`;
description "Geometric density";
    if type(p, numeric) and (p < 0 or 1 < p) then ERROR("F\
        irst argument must be a variable or between 0 and \
        1")
    fi;
    if type(x, numeric) and (not type(x, integer) or x < 1)
    then RETURN(0)
    fi;
    p*(1 - p)^(x - 1)

end:


HypergeometricPDF := proc(n1::algebraic, n2::algebraic, r::algebraic, x::algebraic)
option `January 27, 1994 -- Zaven A. Karian`;
description "Hypergeometric density";
    if type(x, numeric) and type(r, numeric) and r < x then
        ERROR("Fourth argument cannot be larger than the t\
        hird argument")
    fi;
    if type(x, numeric) and type(n1, numeric) and n1 < x then
        ERROR("Fourth argument cannot be larger than the f\
        irst argument")
    fi;
    if type(x, numeric) and type(r, numeric) and
    type(n2, numeric) and n2 < r - x then ERROR("2nd argum\
        ent cannot be smaller than the 3rd-4th arguments")
    fi;
    binomial(n1, x)*binomial(n2, r - x)/binomial(n1 + n2, r)

end:


LogisticPDF := proc(x::algebraic)
option `July 16, 1993 -- Zaven A. Karian`;
description "Logistic density";
    exp(-x)/(1 + exp(-x))^2

end:


NegBinomialPDF := proc(r::algebraic, p::algebraic, x::algebraic)
option `July 7, 1993 -- Zaven A. Karian`;
description "Negative binomial density";
    if type(p, numeric) and (p < 0 or 1 < p) then ERROR("S\
        econd argument must be a variable or between 0 and\
         1")
    fi;
    if type(r, numeric) and not type(r, integer) then ERROR(
        "First argument must be a variable or an integer")
    fi;
    if type(x, numeric) and not type(x, integer) then
        RETURN(0)
    fi;
    if type(x, numeric) and type(r, numeric) and x < r then
        RETURN(0)
    fi;
    binomial(x - 1, r - 1)*p^r*(1 - p)^(x - r)

end:


NormalPDF := proc(mu::algebraic, var::algebraic, x::algebraic)
local a;
option `August 11, 1993 -- Zaven A. Karian`;
description "Normal density";
    if type(var, numeric) and var <= 0 then ERROR(
        "Second argument must be a variable or positive")
    fi;
    a := 1/sqrt(2*Pi*var);
    a*exp(- 1/2*(x - mu)^2/var)

end:


PoissonPDF := proc(l::algebraic, x::algebraic)
option `August 11, 1993 -- Zaven A. Karian`;
description "Poisson density";
    if type(l, numeric) and l <= 0 then ERROR(
        "First argument must be a variable or positive")
    fi;
    if type(x, numeric) and (not type(x, integer) or x < 0)
    then RETURN(0)
    fi;
    l^x*exp(-l)/x!

end:


TPDF := proc(df::algebraic, x::algebraic)
local C1, C2, C3;
option `August 11, 1993 -- Zaven A. Karian`;
description "Students t density";
    if type(df, numeric) and not type(df, posint) then ERROR(
        "First argument must be a variable or a positive i\
        nteger")
    fi;
    C1 := GAMMA(1/2*df + 1/2);
    C2 := GAMMA(1/2*df)*sqrt(df*Pi);
    C3 := (1 + x^2/df)^(1/2*df + 1/2);
    C1/(C2*C3)

end:


UniformPDF := proc(R::range, x::algebraic)
local a, b;
option `August 12, 1993 -- Zaven A. Karian`;
description "Uniform density";
    a := lhs(R);
    b := rhs(R);
    if type(x, numeric) and type(a, numeric) and x < a then
        RETURN(0)
    fi;
    if type(x, numeric) and type(b, numeric) and b < x then
        RETURN(0)
    fi;
    1/(b - a)

end:


WeibullPDF := proc(a::algebraic, b::algebraic, x::algebraic)
local A;
option `August 11, 1993 -- Zaven A. Karian`;
description "Weibull density";
    if type(a, numeric) and a <= 0 then ERROR(
        "First argument must be a variable or positive")
    fi;
    if type(b, numeric) and b <= 0 then ERROR(
        "Second argument must be a variable or positive")
    fi;
    if type(x, numeric) and x < 0 then RETORN(0) fi;
    A := a*(x/b)^(a - 1)/b;
    A*exp(-(x/b)^a)

end:


ExponentialSumS := proc(t::numeric, n::posint, m::posint)
local X, Y, i, j;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates sums from exponential";
    if t <= 0 then ERROR("First argument must be positive")
    fi;
    Y := array(1 .. m);
    for i to m do
        X := ExponentialS(t, n);
        Y[i] := sum(X[j], j = 1 .. n)
    od;
    [seq(Y[i], i = 1 .. m)]

end:


NormTransVarS := proc(a::numeric, b::numeric, n::posint, m::posint)
local X, Y, i;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates transformed variances from normal";
    if b <= 0 then ERROR("Second argument must be positive")
    fi;
    Y := array(1 .. m);
    for i to m do
        X := NormalS(a, b, n);
        Y[i] := ((n - 1)*Variance(X))/b
    od;
    [seq(Y[i], i = 1 .. m)]

end:


NormalMeanS := proc(a::numeric, b::numeric, n::posint, m::posint)
option `May 21, 1994 -- Zaven A. Karian`;
description "Generates means from normal";
    if b <= 0 then ERROR("Second argument must be positive")
    fi;
    NormalS(a, b/n, m)

end:


NormalMedianS := proc(a::numeric, b::numeric, n::posint, m::posint)
local X, Y, i;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates medians from normal";
    if b <= 0 then ERROR("Second argument must be positive")
    fi;
    Y := array(1 .. m);
    for i to m do
        X := NormalS(a, b, n); Y[i] := Percentile(X, .5)
    od;
    [seq(Y[i], i = 1 .. m)]

end:


NormalSumS := proc(a::numeric, b::numeric, n::posint, m::posint)
option `May 21, 1994 -- Zaven A. Karian`;
description "Generates sums from normal";
    if b <= 0 then ERROR("Second argument must be positive")
    fi;
    NormalS(a, n*b, m)

end:


NormalVarianceS := proc(a::numeric, b::numeric, n::posint, m::posint)
local X, Y, i;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates variances from normal";
    if b <= 0 then ERROR("Second argument must be positive")
    fi;
    Y := array(1 .. m);
    for i to m do X := NormalS(a, b, n); Y[i] := Variance(X)
    od;
    [seq(Y[i], i = 1 .. m)]

end:


UniformMeanS := proc(R::range, n::posint, m::posint)
local X, Y, i;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates means from continuous uniform";
    if not (type(lhs(R), numeric) and type(rhs(R), numeric))
     or rhs(R) < lhs(R) then
        ERROR("First argument must be a numeric range")
    fi;
    Y := array(1 .. m);
    for i to m do X := UniformS(R, n); Y[i] := Mean(X) od;
    [seq(Y[i], i = 1 .. m)]

end:


UniformMedianS := proc(R::range, n::posint, m::posint)
local X, Y, i;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates medians from continuous uniform";
    if not (type(lhs(R), numeric) and type(rhs(R), numeric))
     or rhs(R) < lhs(R) then
        ERROR("First argument must be a numeric range")
    fi;
    Y := array(1 .. m);
    for i to m do
        X := UniformS(R, n); Y[i] := Percentile(X, .5)
    od;
    [seq(Y[i], i = 1 .. m)]

end:


UniformSumS := proc(R::range, n::posint, m::posint)
local X, Y, i;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates sums from continuous uniform";
    if not (type(lhs(R), numeric) and type(rhs(R), numeric))
     or rhs(R) < lhs(R) then
        ERROR("First argument must be a numeric range")
    fi;
    Y := array(1 .. m);
    for i to m do
        X := UniformS(R, n); Y[i] := convert(X, `+`)
    od;
    [seq(Y[i], i = 1 .. m)]

end:

 
BernoulliS := proc(p::numeric, n::posint)
local P;
option `August 12, 1993 -- Zaven A. Karian`;
description "n random nos. from Bernoulli dist.";
    if p < 0 or 1 < p then
        ERROR("First argument must be between 0 and 1")
    fi;
    P := [0, 1 - p, 1, p];
    DiscreteS(P, n)

end:


BetaS := proc(n1::numeric, n2::numeric, n::posint)
local i, y, yy, A, X;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates sample from beta";
    if n1 <= 0 then ERROR("First argument must be positive")
    fi;
    if n2 <= 0 then ERROR("Second argument must be positive")
    fi;
    X := array(1 .. n);
    A := FS(2*n2 + 2, 2*n1 + 2, n);
    for i to n do
        yy := ((n2 + 1)*A[i])/(n1 + 1);
        y := 1/(yy + 1);
        X[i] := y
    od;
    [seq(X[i], i = 1 .. n)]

end:


BinomialS := proc(N::posint, p::numeric, n::posint)
local A, i, L;
option `August 12, 1993 -- Zaven A. Karian`;
description "n random nos. from B(n,p)";
    if p < 0 or 1 < p then
        ERROR("Second argument must be between 0 and 1")
    fi;
    A := array(1 .. 2*N + 2);
    for i from 0 to N do
        A[2*i + 1] := i;
        A[2*i + 2] := binomial(N, i)*p^i*(1 - p)^(N - i)
    od;
    L := [seq(A[i], i = 1 .. 2*N + 2)];
    DiscreteS(L, n)

end:


BivariateNormalS := proc(mux::numeric, varx::numeric,muy::numeric, 
    vary::numeric, rho::numeric, n::posint)
local i, u, v, a, b, cmua, cmuy, cvar, A;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates sample from bivariate normal";
    if varx <= 0 then
        ERROR("Second argument must be positive")
    fi;
    if vary <= 0 then
        ERROR("Fourth argument must be positive")
    fi;
    if rho < -1 or 1 < rho then
        ERROR("Fifth argument must be between -1 and +1")
    fi;
    cvar := vary*(1 - rho^2);
    cmua := rho*sqrt(vary/varx);
    A := array(1 .. n, 1 .. 2);
    for i to n do
        u := rng();
        v := rng();
        a := sqrt(-2*log(u))*sin(2*Pi*v);
        a := evalf(sqrt(varx)*a + mux);
        A[i, 1] := a;
        cmuy := muy + cmua*(a - mux);
        b := sqrt(-2*log(u))*cos(2*Pi*v);
        b := evalf(sqrt(cvar)*b + cmuy);
        A[i, 2] := b
    od;
    [seq([A[i, 1], A[i, 2]], i = 1 .. n)]

end:


CauchyS := proc(n::posint)
local i;
option `July 6, 1993 -- Zaven A. Karian`;
description "Generates sample from Cauchy";
    [seq(evalf(tan(Pi*rng() - 1/2*Pi)), i = 1 .. n)]

end:


ChisquareS := proc(r::posint, n::posint)
local i, j, X, A, t;
option `May 27, 1994 -- Zaven A. Karian`;
description "Generates sample from standard normal";
    A := array(1 .. n);
    X := N01S(r*n);
    j := 1;
    for i to n do
        A[i] :=
            convert([seq(X[t]^2, t = j .. j + r - 1)], `+`);
        j := j + r
    od;
    [seq(A[i], i = 1 .. n)]

end:


ContinuousS := proc(expr::algebraic, R::range, n::posint)
local i, var, x, F, CDF, t, ll, ul, pts, Y, j, Rand, a, b, c, d;
option `January 24, 1994 -- Zaven A. Karian`;
description "Gen sample (size n) with pdf given by expr on \
the interval R";
    var := op(select(type, indets(expr), name));
    ll := lhs(args[2]);
    ul := rhs(args[2]);
    F := int(expr, var = ll .. t);
    CDF := student[makeproc](F, t);
    pts := array(1 .. 72);
    Y := array(1 .. n);
    j := 3;
    for i from 2 to 35 do
        x := ll + 1/35*(i - 1)*(ul - ll);
        pts[j] := CDF(x);
        pts[j + 1] := x;
        j := j + 2
    od;
    pts[1] := 0;
    pts[2] := ll;
    pts[71] := 1;
    pts[72] := ul;
    Rand := RNG(n);
    for i to n do
        j := 1;
        while pts[2*j - 1] < Rand[i] do j := j + 1 od;
        a := pts[2*j - 3];
        b := pts[2*j - 2];
        c := pts[2*j - 1];
        d := pts[2*j];
        Y[i] := b + ((Rand[i] - a)*(d - b))/(c - a)
    od;
    [seq(Y[i], i = 1 .. n)]

end:


Die := proc(m::posint, n::posint)
local die, i;
option `July 5, 1993 -- Zaven A. Karian`;
description "n rolls of an m-sided die";
    die := rand(1 .. m); [seq(die(), i = 1 .. n)]

end:


DiscUniformS := proc(R::range, n::posint)
local i, a, b;
option `July 6, 1993 -- Zaven A. Karian`;
description "Generates sample from discrete uniform";
    a := lhs(R);
    b := rhs(R);
    [seq(a + trunc((b - a + 1)*rng()), i = 1 .. n)]

end:


DiscreteS := proc()
local N, A, B, ll, ul, ptr, var, F, i, r, j, L, expr, n;
option `July 5, 1993 -- Zaven A. Karian`;
description "n rand. obs. from dist. described by L";
    if nargs = 2 and not type(args[1], list) then ERROR("W\
        hen two arguments are used, the first must be a li\
        st")
    fi;
    if nargs = 2 and not type(args[2], posint) then ERROR(
        "When two arguments are used, the second must be a\
         positive integer")
    fi;
    if nargs = 3 and not type(args[1], algebraic) then ERROR(
        "When three arguments are used, the first must be \
        an expression")
    fi;
    if nargs = 3 and not type(args[2], range) then ERROR("\
        When three arguments are used, the second must be \
        a range")
    fi;
    if nargs = 3 and not type(args[3], posint) then ERROR(
        "When three arguments are used, the third must be \
        a positive integer")
    fi;
    if nargs = 3 then
        expr := args[1];
        n := args[3];
        var := op(select(type, indets(expr), name));
        ll := lhs(args[2]);
        ul := rhs(args[2]);
        B := array(1 .. 2*ul - 2*ll + 2);
        ptr := 1;
        for i from ll to ul do
            B[ptr] := i;
            B[ptr + 1] := subs(var = i, expr);
            ptr := ptr + 2
        od;
        L := convert(B, list)
    else L := args[1]; n := args[2]
    fi;
    N := 1/2*nops(L);
    A := array(1 .. n);
    F := array(1 .. N);
    F[1] := L[2];
    for i from 2 to N do F[i] := F[i - 1] + L[2*i] od;
    for i to n do
        r := rng();
        if r < F[1] then A[i] := L[1] fi;
        if F[N] < r then A[i] := L[2*N - 1] fi;
        for j from 2 to N do
            if r < F[j] and F[j - 1] <= r then
                A[i] := L[2*j - 1]
            fi
        od
    od;
    [seq(A[i], i = 1 .. n)]

end:


ExponentialS := proc(t::numeric, n::posint)
local i;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates sample from exponential";
    if t <= 0 then ERROR("First argument must be positive")
    fi;
    [seq(evalf(- t*log(1 - rng())), i = 1 .. n)]

end:


FS := proc(n1::posint, n2::posint, n::posint)
local i, A, B;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates sample from F";
    A := ChisquareS(n1, n);
    B := ChisquareS(n2, n);
    [seq(A[i]*n2/(n1*B[i]), i = 1 .. n)]

end:


GammaS := proc(a::numeric, t::numeric, n::posint)
local WW, V, X, ay, Zee, d, q, theta, k, j, Y, U1, U2, b, P;
    if a <= 0 then ERROR("First argument must be positive")
    fi;
    if t <= 0 then ERROR("Second argument must be positive")
    fi;
    X := array(1 .. n);
    k := 1;
    while k <= n do
        if a < 1 then
            b := evalf((exp(1) + a)/exp(1));
            U1 := RNG(1)[1];
            P := b*U1;
            if 1 < P then
                Y := evalf(-ln((b - P)/a));
                U2 := RNG(1)[1];
                if U2 <= evalf(Y^(a - 1)) then
                    X[k] := t*Y; k := k + 1
                fi
            else
                Y := P^(1/a);
                U2 := RNG(1)[1];
                if U2 <= evalf(exp(-Y)) then
                    X[k] := t*Y; k := k + 1
                fi
            fi
        else
            ay := evalf(1/sqrt(2*a - 1));
            b := evalf(a - ln(4));
            q := evalf(a + 1/ay);
            theta := 4.5;
            d := 1 + ln(theta);
            U1 := RNG(1)[1];
            U2 := RNG(1)[1];
            V := evalf(ay*ln(U1/(1 - U1)));
            Y := evalf(a*exp(V));
            Zee := U1^2*U2;
            WW := b + q*V - Y;
            if 0 <= evalf(WW + d - theta*Zee) then
                X[k] := t*Y; k := k + 1
            else
                if evalf(ln(Zee)) <= WW then
                    X[k] := t*Y; k := k + 1
                fi
            fi
        fi
    od;
    [seq(X[j], j = 1 .. n)]

end:


GeometricS := proc(p::numeric, n::posint)
local A, SumProb, i, Prob, UpperLimit, L;
option `August 12, 1993 -- Zaven A. Karian`;
description "n random nos. from geometric dist.";
    if p < 0 or 1 < p then
        ERROR("First argument must be between 0 and 1")
    fi;
    A := array(1 .. 500);
    SumProb := 0;
    i := 1;
    if p = 1 then A[1] := 1; A[2] := 1; UpperLimit := 2
    else
        while SumProb < .999 do
            Prob := p*(1 - p)^(i - 1);
            SumProb := SumProb + Prob;
            A[2*i - 1] := i;
            A[2*i] := Prob;
            i := i + 1
        od;
        A[2*i - 1] := i;
        A[2*i] := 1 - SumProb;
        UpperLimit := 2*i
    fi;
    L := [seq(A[i], i = 1 .. UpperLimit)];
    DiscreteS(L, n)

end:


HypergeometricS :=proc(n1::posint, n2::posint, r::posint, n::posint)
local lb, ub, A, ptr, x, L, i;
option `January 27, 1994 -- Zaven A. Karian`;
description "Hypergeometric density";
    lb := max(0, r - n2);
    ub := min(r, n1);
    A := array(1 .. 2*ub - 2*lb + 2);
    ptr := 1;
    for x from lb to ub do
        A[ptr] := x;
        A[ptr + 1] := binomial(n1, x)*binomial(n2, r - x)/
            binomial(n1 + n2, r);
        ptr := ptr + 2
    od;
    L := [seq(A[i], i = 1 .. 2*ub - 2*lb + 2)];
    DiscreteS(L, n)

end:


LogisticS := proc(n::posint)
local a, r, i, X;
option `July 6, 1993 -- Zaven A. Karian`;
description "Generates sample from logistic";
    X := array(1 .. n);
    for i to n do
        r := rng(); a := evalf(log(r/(1 - r))); X[i] := a
    od;
    [seq(X[i], i = 1 .. n)]

end:


N01S := proc(n::posint)
local j, i, u, v, X, P2;
option `July 6, 1993 -- Zaven A. Karian`;
description "Generates sample from standard normal";
    X := array(1 .. n + 1);
    j := 1;
    P2 := evalf(2*Pi);
    for i to 1/2*n + 1/2 do
        u := rng();
        v := rng();
        X[j] := evalf(sqrt(-2*log(u))*sin(P2*v));
        j := j + 1;
        X[j] := evalf(sqrt(-2*log(u))*cos(P2*v));
        j := j + 1
    od;
    [seq(X[i], i = 1 .. n)]

end:


NegBinomialS := proc(r::posint, p::numeric, n::posint)
local A, SumProb, i, Prob, UpperLimit, L;
option `August 12, 1993 -- Zaven A. Karian`;
description "n random nos. from negative binomial dist.";
    if p < 0 or 1 < p then
        ERROR("Second argument must be between 0 and 1")
    fi;
    A := array(1 .. 100);
    SumProb := 0;
    i := r;
    while SumProb < .999 do
        Prob :=
            evalf(binomial(i - 1, r - 1)*p^r*(1 - p)^(i - r))
            ;
        SumProb := SumProb + Prob;
        A[2*i - 2*r + 1] := i;
        A[2*i - 2*r + 2] := Prob;
        i := i + 1
    od;
    A[2*i - 2*r + 1] := i;
    A[2*i - 2*r + 2] := 1 - SumProb;
    UpperLimit := 2*i - 2*r + 2;
    L := [seq(A[i], i = 1 .. UpperLimit)];
    DiscreteS(L, n)

end:


NormalS := proc(mu::numeric, var::numeric, n::posint)
local j, i, P2, u, v, a, b, X, sqrt_var;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates sample from normal";
    if var <= 0 then
        ERROR("Second argument must be positive")
    fi;
    X := array(1 .. n + 1);
    j := 1;
    P2 := evalf(2*Pi);
    sqrt_var := sqrt(var);
    for i to 1/2*n + 1/2 do
        u := rng();
        v := rng();
        a := sqrt(-2*log(u))*sin(P2*v);
        X[j] := evalf(sqrt_var*a + mu);
        j := j + 1;
        b := sqrt(-2*log(u))*cos(P2*v);
        X[j] := evalf(sqrt_var*b + mu);
        j := j + 1
    od;
    [seq(X[i], i = 1 .. n)]

end:


PoissonS := proc(l::numeric, n::posint)
local A, SumProb, i, Prob, UpperLimit, L;
option `July 5, 1993 -- Zaven A. Karian`;
description "n random nos. from Poisson dist.";
    if l <= 0 then ERROR("First argument must be positive")
    fi;
    A := array(1 .. 100);
    SumProb := 0;
    i := 0;
    while SumProb < .999 do
        Prob := evalf(l^i*exp(-l)/i!);
        SumProb := SumProb + Prob;
        A[2*i + 1] := i;
        A[2*i + 2] := Prob;
        i := i + 1
    od;
    A[2*i + 1] := i;
    A[2*i + 2] := 1 - SumProb;
    UpperLimit := 2*i + 2;
    L := [seq(A[i], i = 1 .. UpperLimit)];
    DiscreteS(L, n)

end:


RNG := proc(n::posint)
local i;
option `July 5, 1993 -- Zaven A. Karian`;
description "n random numbers from [0,1]";
    [seq(Float(rand(), -12), i = 1 .. n)]

end:


TS := proc(nu::posint, n::posint)
local i, A, B;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates sample from Student's t";
    A := N01S(n);
    B := ChisquareS(nu, n);
    [seq(A[i]/sqrt(B[i]/nu), i = 1 .. n)]

end:


UniformS := proc(R::range, n::posint)
local a, b, i;
option `August 12, 1993 -- Zaven A. Karian`;
description "Generates sample from continuous uniform";
    a := lhs(R);
    b := rhs(R);
    if not (type(a, numeric) and type(b, numeric)) or b < a
    then ERROR(
        "First argument must be a proper numeric range")
    fi;
    [seq(evalf((b - a)*rng() + a), i = 1 .. n)]

end:


MakeRandom := proc()
global _seed;
option `July 5, 1993 -- Zaven A. Karian`;
description "n random numbers from [0,1)";
    _seed := round(123456789*time())

end:


rng := () -> Float(rand(), -12):


ClassFreq := proc(L::list, R::range, n::posint)
local N, i, A, SL, bdry, delta, ptr;
option `January 26, 1994 -- Zaven A. Karian`;
description
"Class Frequencies of L on the range R in n intervals";
    if nargs = 2 then N := 8 else N := n fi;
    A := array(1 .. N);
    SL := sort(L);
    delta := (rhs(R) - lhs(R))/N;
    ptr := 1;
    for i to N do
        A[i] := 0;
        bdry := i*delta + lhs(R);
        while ptr <= nops(SL) and SL[ptr] <= bdry do
            A[i] := A[i] + 1; ptr := ptr + 1
        od
    od;
    [seq(A[i], i = 1 .. N)]

end:


Freq := proc(L::list, R::range)
local count, N, i, Min, Max;
option `August 12, 1993 -- Zaven A. Karian`;
description "Frequencies of L on the range R";
    Min := lhs(R);
    Max := rhs(R);
    if not (type(Min, numeric) and type(Max, numeric)) or
    Max < Min then ERROR(
        "Second argument must be a propoer numerical range")
    fi;
    count := array(Min .. Max);
    N := nops(L);
    for i from Min to Max do count[i] := 0 od;
    for i to N do
        if Min <= L[i] and L[i] <= Max then
            count[L[i]] := count[L[i]] + 1
        fi
    od;
    [seq(count[i], i = Min .. Max)]

end:


Kurtosis := proc(L::list)
local sf, M, N, V, i, X;
option `August 19, 1993 -- Zaven A. Karian`;
description "Kurtosis or 3rd standardized moment";
    if type(args[1], listlist) then
        X := [seq(L[i][1] $ L[i][2], i = 1 .. nops(L))]
    else X := L
    fi;
    sf := 0;
    M := Mean(X);
    N := nops(X);
    V := ((N - 1)*Variance(X))/N;
    for i to N do sf := sf + (X[i] - M)^4 od;
    sf/(N*V^2)

end:


Locate := proc(X::list, item::algebraic)
local n, A, ptr, i, j;
option `January 24, 1994 -- Zaven A. Karian`;
description "Find the indeces in the list X of the item";
    n := nops(X);
    A := array(1 .. n);
    ptr := 1;
    for i to n do
        if X[i] = item then A[ptr] := i; ptr := ptr + 1 fi
    od;
    [seq(A[j], j = 1 .. ptr - 1)]

end:


Max := proc(L::list)
option `July 5, 1993 -- Zaven A. Karian`;
description "Maximum value of the list L";
    max(op(L))

end:


Mean := proc(L::list)
local S, N, i;
option `July 5, 1993 -- Zaven A. Karian`;
description "Arithmetic mean of the list L";
    if type(args[1], listlist) then
        S := 0;
        N := 0;
        for i to nops(L) do
            S := S + L[i][1]*L[i][2]; N := N + L[i][2]
        od
    else S := convert(L, `+`); N := nops(L)
    fi;
    S/N

end:

Median := proc(L::list)
option `May 26, 1994 -- Zaven A. Karian`;
description "The 1n of the list L";
    Percentile(L, .5)

end:


Min := proc(L::list)
option `July 5, 1993 -- Zaven A. Karian`;
description "Minimum value of the list L";
    min(op(L))

end:


Percentile := proc(L::list(numeric), p::numeric)
local n, LL, F, r, ab;
option `August 12, 1993 -- Zaven A. Karian`;
description "The 100*p-th percentile of the list L";
    if p < 0 or 1 < p then
        ERROR("Second argument must be between 0 and 1")
    fi;
    n := nops(L);
    if p < 1/(n + 1) or n/(n + 1) < p then
        ERROR("Percentile does not exist")
    fi;
    LL := sort(L);
    F := convert((n + 1)*p, fraction);
    r := trunc(F);
    ab := F - r;
    LL[r] + ab*(LL[r + 1] - LL[r])

end:


Range := proc(L::list)
option `July 5, 1993 -- Zaven A. Karian`;
description "Range of the list L";
    Max(L) - Min(L)

end:


RunningSum := proc(L::list)
local N, A, i;
option `August 23, 1993 -- Zaven A. Karian`;
description "Running sum of a list";
    N := nops(L);
    A := array(1 .. N);
    A[1] := L[1];
    for i from 2 to N do A[i] := A[i - 1] + L[i] od;
    [seq(A[i], i = 1 .. N)]

end:


Skewness := proc(L::list)
local sc, M, N, V, i, X;
option `August 19, 1993 -- Zaven A. Karian`;
description "Skewness or 3rd standardized moment";
    if type(args[1], listlist) then
        X := [seq(L[i][1] $ L[i][2], i = 1 .. nops(L))]
    else X := L
    fi;
    sc := 0;
    M := Mean(X);
    N := nops(X);
    V := ((N - 1)*Variance(X))/N;
    for i to N do sc := sc + (X[i] - M)^3 od;
    sc/(N*V^(3/2))

end:


StDev := proc(L::list)
option `July 5, 1993 -- Zaven A. Karian`;
description "Standard deviation of the list L";
    sqrt(Variance(L))

end:


StemLeaf := proc(L::list(numeric))
local n, stems, digits, Ls, range, logscale, x, fracwidth,
minwidth, maxwidth, M, Label, labels, D, i, width, stemnum,
start, format, leaf, stem, dispstem, tenstem;
option `Joshua D. Levy  ~  June 10, 1994`;
description "Ordered stem-and-leaf display";
    M := [1, 2, 5, 10];
    Label :=
        '[[" "], ["*", t, f, s, "+"], ["*", "+"], [" "]]';
    n := nops(L);
    Ls := sort(evalf(L));
    range := Ls[n] - Ls[1];
    if 1 < nargs and type(args[2], posint) then
        digits := args[2]
    else digits := min(length(op(1, range)), 2)
    fi;
    if 2 < nargs and type(args[3], posint) then
        stems := args[3]
    else stems := round(evalf(sqrt(n)))
    fi;
    logscale := digits - ilog10(range/stems) - 1;
    Ls := [seq(round(10^logscale*x), x = Ls)];
    fracwidth := (Ls[n] - Ls[1])/stems;
    minwidth := 10^ilog10(fracwidth);
    maxwidth := 10*minwidth;
    D := seq(evalf(abs(fracwidth - M[i]*minwidth)),
        i = 1 .. 4);
    member(min(D), [D], 'i');
    width := M[i]*minwidth;
    labels := nops(Label[i]);
    stemnum := iquo(Ls[1], width);
    start := width*stemnum;
    leaf := 1;
    format := cat(" %0", length(minwidth), ".0f");
    printf("\n   Stem   Leaf");
    for stem from start by width while stem <= Ls[n] do
        dispstem := iquo(stem, maxwidth, 'r');
        tenstem := maxwidth*dispstem;
        printf("\n%6.0f%c |", dispstem,
            Label[i][1 + irem(stemnum, labels)]);
        for leaf from leaf while
        leaf <= n and Ls[leaf] < stem + width do
            printf(format, Ls[leaf] - tenstem)
        od;
        stemnum := stemnum + 1
    od;
    if logscale <> 0 then printf("\n                      \
         %.0f\n(Multiply numbers by 10 .)", -logscale)
    fi;
    lprint()

end:


Variance := proc(L::list)
local ssq, s, N, i, X;
option `July 5, 1993 -- Zaven A. Karian`;
description "Variance of the list L";
    if type(args[1], listlist) then
        X := [seq(L[i][1] $ L[i][2], i = 1 .. nops(L))]
    else X := L
    fi;
    N := nops(X);
    s := convert(X, `+`);
    ssq := convert([seq(X[i]^2, i = 1 .. N)], `+`);
    (ssq - s^2/N)/(N - 1)

end:

end module: