{VERSION 4 0 "IBM INTEL NT" "4.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "" -1 256 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 257 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "Courier" 1 10 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE " Heading 1" -1 3 1 {CSTYLE "" -1 -1 "Times" 1 18 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 4 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 2" -1 4 1 {CSTYLE "" -1 -1 "Times" 1 14 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 2 1 0 1 0 2 2 0 1 }{PSTYLE "Bullet Item" -1 15 1 {CSTYLE "" -1 -1 "Tim es" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 3 3 1 0 1 0 2 2 15 2 } {PSTYLE "Title" -1 18 1 {CSTYLE "" -1 -1 "Times" 1 18 0 0 0 1 2 1 1 2 2 2 1 1 1 1 }3 1 0 0 12 12 1 0 1 0 2 2 19 1 }{PSTYLE "Author" -1 19 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 8 8 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 18 "" 0 "" {TEXT -1 34 "The Statistics Supplement to Maple" }}{PARA 19 "" 0 "" {TEXT -1 78 "by Zavan Karian, karian@den ison.edu, http://www.denison.edu/~karian/index.html" }}}{EXCHG {PARA 3 "" 0 "" {TEXT -1 39 "About the Statistics Supplement package" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 110 "This pac kage provides over 130 Maple functions written for exploring concepts \+ in probability and statistics. " }}{PARA 0 "" 0 "" {TEXT -1 82 "It is intended for instructors and students in Probability and Statistics c ourses." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 72 "This worksheet is for Maple 6. A version for Maple 5 is also avai lable." }}}{EXCHG {PARA 3 "" 0 "" {TEXT -1 35 "How to install and use \+ this package" }}{PARA 15 "" 0 "" {TEXT -1 86 "Scroll to the bottom of \+ this worksheet and follow the installation instructions there." }} {PARA 15 "" 0 "" {TEXT -1 108 "Once the package is installed, you can \+ use the package in any Maple worksheet by entering the Maple command \+ " }{TEXT 256 11 "with(stat) " }{TEXT -1 14 "after setting " }{TEXT 265 8 "libname " }{TEXT -1 57 "to the directory where the package was \+ saved (see below)." }}{PARA 15 "" 0 "" {TEXT -1 98 "Documentation on e ach of the functions in this package is provided in the Microsoft Word document " }{TEXT 257 23 "stats_documentation.doc" }{TEXT -1 1 "." }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 99253 "stat := module()\n\n \+ export Anova1, Anova2m, Anova2s, AnovaBeta, BernoulliCDF, BernoulliPDF , \n BernoulliS, BetaCDF, BetaP, BetaPDF, BetaS, BinomialCDF, B inomialPDF, \n BinomialS, BivariateNormalPDF, BivariateNormalS, BoxWhisker, CauchyPDF, \n CauchyS, ChisquareCDF, ChisquareFit, ChisquareP, ChisquarePDF, ChisquareS, \n ClassFreq, ConfIntAvL en, ConfIntMean, ConfIntPlot, ConfIntProp, ConfIntSuc, \n ConfI ntVar, Contingency, ContinuousS, Convolution, Correlation, Craps, \n \+ Die, DiscUniformS, DiscreteS, DoubleExponentialPDF, ExponentialC DF, \n ExponentialP, ExponentialPDF, ExponentialS, ExponentialS umS, FCDF, FP, \n FPDF, FS, FastCraps, Freq, GBCDF, GBDPDF, GBP , GammaCDF, GammaP, GammaPDF, \n GammaS, GeometricCDF, Geometri cPDF, GeometricS, GraphRandWalk, GraphRandWalk2, \n Histogram, \+ HypergeometricPDF, HypergeometricS, KSFit, Kurtosis, LinReg, \n \+ Locate, LogisticPDF, LogisticS, MarginalRelFreq, Max, Mean, Median, M edianTest, \n Min, N01S, NegBinomialCDF, NegBinomialPDF, NegBin omialS, NormTransVarS, \n NormalCDF, NormalMeanS, NormalMedianS , NormalP, NormalPDF, NormalS, NormalSumS, \n NormalVarianceS, \+ Ogive, Percentile, PlotDiscCDF, PlotEmpCDF, PlotEmpPDF, \n Plot PolyReg, PlotRunningAverage, PoissonCDF, PoissonPDF, PoissonS, PolyReg , \n ProbHist, QQ, QQFit, RNG, RandWalk, RandWalk2, Range, RegA nal, RegBand, \n Residuals, RunTest, RunningSum, ScatPlot, Scat PlotLine, SignTest, Skewness, \n StDev, StemLeaf, TCDF, TP, TPD F, TS, TimePlot, UniformCDF, UniformMeanS, \n UniformMedianS, U niformP, UniformPDF, UniformS, UniformSumS, Variance, \n Weibul lPDF, Wilcoxon, Wilcoxon2, XbarChart, MakeRandom, rng;\n\noption packa ge;\n\n\nChisquareFit := proc(L::list, expr::algebraic, Classes::list, Type::string)\n\nlocal N, i, X, NC, NumC, var, XA, j, Q, Fr, p, pval, C_Dadj,\nPlotSet, P1, P2, lastarg;\nglobal color, blue, red, linestyl e;\n\noption `May 23, 1994 -- May 23, 1994`;\ndescription \"To perform a Chi square goodness of fi\\\nt test for the data in L, the distribu tion with CD\\\nF expr, and class boundaries given by Classes. The \\ \nfourth argument, Type, must be the string D or C to\\\nindicate a di screte or a continuous type.\";\n\n if type(L, listlist) then\n \+ N := sum(L[i][2], i = 1 .. nops(L));\n X := [seq(L[i][1] $ \+ L[i][2], i = 1 .. nops(L))]\n else N := nops(L); X := L\n fi;\n \+ NC := nops(Classes);\n NumC := NC - 1;\n var := op(select(typ e, indets(expr), name));\n if nops([var]) = 0 or 1 < nops([var]) th en \n ERROR(\"The second argument must have exactly one varaibl e in it\")\n fi;\n XA := convert(X, array);\n if nargs = 4 th en C_Dadj := .5\n else C_Dadj := .1*10^(-7)\n fi;\n for i to \+ N do\n if XA[i] < Classes[1] then XA[i] := Classes[1] fi;\n \+ if Classes[NC] <= X[i] then\n XA[i] := Classes[NC] - C_ Dadj\n fi\n od;\n X := sort(convert(XA, list));\n Fr : = array(1 .. NumC);\n for i to NumC do Fr[i] := 0 od;\n j := 1; \n i := 1;\n while j <= N do\n if Classes[i] <= X[j] and \+ X[j] < Classes[i + 1] then\n Fr[i] := Fr[i] + 1; j := j + 1 \n else i := i + 1\n fi\n od;\n if nargs = 4 then \+ C_Dadj := 1 else C_Dadj := 0 fi;\n p := array(1 .. NumC);\n p[1] := evalf(subs(var = Classes[2] - C_Dadj, expr));\n for i from 2 to NumC - 1 do p[i] := evalf(\n subs(var = Classes[i + 1] - C_Dad j, expr)\n - subs(var = Classes[i] - C_Dadj, expr))\n od;\n p[NumC] :=\n evalf(1 - subs(var = Classes[NumC] - C_Dadj, e xpr));\n Q := 0;\n for i to NumC do Q := Q + (Fr[i] - N*p[i])^2/ (N*p[i]) od;\n for i to NumC do\n if N*p[i] < 4 then\n \+ printf(\"Warning: Class %3.0f has an expected frequency of %f\", i, N*p[i]);\n lprint()\n fi\n od;\n if nargs \+ = 4 then C_Dadj := .5 else C_Dadj := 0 fi;\n if nargs = 4 then last arg := NumC else lastarg := NC fi;\n PlotSet := \{[[Classes[NumC] - C_Dadj, 0], [Classes[NumC] - C_Dadj,\n Fr[NumC]/(N*(Classes[NC ] - Classes[NumC]))], [Classes[lastarg] + C_Dadj,\n Fr[NumC]/(N *(Classes[NC] - Classes[NumC]))],\n [Classes[lastarg] + C_Dadj, 0]]\};\n for i to NumC - 1 do PlotSet := PlotSet union \n \+ \{[[Classes[i] - C_Dadj, 0], [Classes[i] - C_Dadj,\n Fr[i]/(N*( Classes[i + 1] - Classes[i]))], [Classes[i + 1] - C_Dadj,\n Fr[ i]/(N*(Classes[i + 1] - Classes[i]))],\n [Classes[i + 1] - C_Da dj, 0]]\}\n od;\n P1 := plot(PlotSet, color = blue, linestyle = \+ 3);\n if nargs = 4 then\n PlotSet := \{[[Classes[NumC] - C_D adj, 0], [Classes[NumC] - C_Dadj,\n p[NumC]/(Classes[NC] - \+ Classes[NumC])], [Classes[NumC] + C_Dadj,\n p[NumC]/(Classe s[NC] - Classes[NumC])],\n [Classes[NumC] + C_Dadj, 0]]\}; \n for i to NumC - 1 do PlotSet := PlotSet union \n \+ \{[[Classes[i] - C_Dadj, 0], [Classes[i] - C_Dadj,\n p[i]/( Classes[i + 1] - Classes[i])], [Classes[i + 1] - C_Dadj,\n \+ p[i]/(Classes[i + 1] - Classes[i])],\n [Classes[i + 1] - C_ Dadj, 0]]\}\n od;\n P2 := plot(PlotSet, color = red)\n \+ else\n diff(expr, var);\n P2 := plot(%, var = Classes[ 1] .. Classes[NC])\n fi;\n with(plots, display);\n print(disp lay(\{P1, P2\}));\n pval := 1 - evalf(ChisquareCDF(NumC - 1, Q));\n printf(\n \"The chi-square goodness of fit statistic is %7. 4f\",\n Q);\n lprint();\n printf(\"The number of classes \+ is %3f\", NumC);\n lprint();\n printf(\"The p-value for %3f\", N umC - 1);\n printf(\" degrees of freedom is %6.4f\", pval);\n lp rint();\n lprint()\n\nend proc;\nKSFit := proc(L::list, expr::algeb raic, R::range)\nlocal var, SL, N, KS, i, F, DOdd, DEven, Loc, P1, P2; \noption `May 23, 1994 -- Zaven A. Karian`;\ndescription \"To perform \+ a Kolmogorov-Smirnov goodness \\\nof fit test for the data in L, the d istribution \\\nwith CDF expr, on the range of values specified by\\\n R.\";\n\n var := op(select(type, indets(expr), name));\n if nops ([var]) = 0 or 1 < nops([var]) then \n ERROR(\"The second argum ent must have exactly one varaible in \\\n it\")\n fi;\n \+ SL := sort(L);\n N := nops(SL);\n KS := 0;\n for i to N do\n \+ F := subs(var = SL[i], expr);\n DOdd := evalf(abs((i - 1 )/N - F));\n DEven := evalf(abs(i/N - F));\n if KS < max (DOdd, DEven) then\n KS := max(DOdd, DEven); Loc := i\n \+ fi\n od;\n P1 := PlotEmpCDF(SL, R);\n P2 := plot(expr, va r = R, color = blue);\n with(plots, display);\n print(display(\{ P1, P2\}));\n printf(\"The K-S goodness of fit statistic of %7.4f\" , KS);\n printf(\" occured at %6.4f\", SL[Loc]);\n lprint();\n \+ lprint()\n\nend proc;\n\n\nAnova1, proc(D::list)\nlocal m, n, tot, S SE, SST, SSTO, N, PV, XBar, i, j, \n GrandMean, MST, MSE;\noption `July 7, 1993 -- Zaven A. Karian`;\ndescription \"One factor ANOVA\"; \n\n m := nops(D);\n n := 0;\n tot := 0;\n SSE := 0;\n \+ SST := 0;\n SST := 0;\n N := array(1 .. m);\n XBar := array(1 .. m);\n for i to m do\n N[i] := nops(D[i]);\n XBar[ i] := Mean(D[i]);\n n := n + N[i]\n od;\n for i to m do f or j to N[i] do\n tot := tot + D[i][j];\n SSE := SSE + (D[i][j] - XBar[i])^2\n od\n od;\n GrandMean := to t/n;\n for i to m do SST := SST + N[i]*(XBar[i] - GrandMean)^2\n \+ od;\n SSTO := SSE + SST;\n MST := SST/(m - 1);\n MSE := SSE/ (n - m);\n PV := 1 - FCDF(m - 1, n - m, MST/MSE);\n lprint(` \+ Sum of`);\n lprint(` Squares Degrees of Mean Sq `);\n lprint(`Source (SS) Freedom (MS) F-\\\n \+ Ratio p-value`);\n lprint(`--------------------------------- -------------\\\n ---------------`);\n printf(\"Treat. %10. 4f %5d %10.4f %8.4f %8.4f\",\n SST, m - 1, MST, MST/MSE , PV);\n lprint();\n printf(\"Error %10.4f %5d %10.4f\", S SE, n - m, MSE);\n lprint();\n printf(\"Total %10.4f %5d\", \+ SSTO, n - 1);\n lprint();\n lprint(`---------------------------- ------------------\\\n ---------------`)\n\nend proc;\n\nAnova2 m := proc(D::listlist)\nlocal a, b, c, RowMean, ColMean, CellMean, Gra ndMean, SSTO,\nSSA, SSB, SSAB, SSE, i, j, k, MSA, MSB, MSAB, MSE, PVR, PVC,\nPVInt;\noption `July 5, 1993 -- Zaven A. Karian`;\ndescription \+ \"Two factro ANOVA (multiple obs. per cell)\";\n a := nops(D);\n \+ b := nops(D[1]);\n c := nops(D[1][1]);\n RowMean := array(1 .. \+ a);\n ColMean := array(1 .. b);\n CellMean := array(1 .. a, 1 .. b);\n GrandMean := 0;\n SSTO := 0;\n SSA := 0;\n SSB := 0 ;\n SSAB := 0;\n SSE := 0;\n for i to a do\n RowMean[i ] := 0;\n for j to b do for k to c do\n RowMean[ i] := RowMean[i] + D[i][j][k]\n od\n od;\n Ro wMean[i] := RowMean[i]/(b*c)\n od;\n for j to b do\n ColM ean[j] := 0;\n for i to a do for k to c do\n Col Mean[j] := ColMean[j] + D[i][j][k]\n od\n od;\n \+ ColMean[j] := ColMean[j]/(a*c)\n od;\n for i to a do for j to b do\n CellMean[i, j] := 0;\n for k to c do\n \+ CellMean[i, j] := CellMean[i, j] + D[i][j][k]\n \+ od;\n CellMean[i, j] := CellMean[i, j]/c\n od\n \+ od;\n for i to a do for j to b do for k to c do\n \+ GrandMean := GrandMean + D[i][j][k]\n od\n od\n o d;\n GrandMean := GrandMean/(a*b*c);\n for i to a do for j to b \+ do for k to c do\n SSTO := SSTO + (D[i][j][k] - GrandMe an)^2;\n SSE := SSE + (D[i][j][k] - CellMean[i, j])^2\n od\n od\n od;\n for i to a do SSA := SSA + (R owMean[i] - GrandMean)^2 od;\n for j to b do SSB := SSB + (ColMean[ j] - GrandMean)^2 od;\n SSA := b*c*SSA;\n SSB := a*c*SSB;\n S SAB := SSTO - SSA - SSB - SSE;\n MSA := SSA/(a - 1);\n MSB := SS B/(b - 1);\n MSAB := SSAB/((a - 1)*(b - 1));\n MSE := SSE/(a*b*( c - 1));\n PVR := 1 - FCDF(a - 1, a*b*(c - 1), MSA/MSE);\n PVC : = 1 - FCDF(b - 1, a*b*(c - 1), MSB/MSE);\n PVInt := 1 - FCDF((a - 1 )*(b - 1), a*b*(c - 1), MSAB/MSE)\n ;\n lprint(` S um of`);\n lprint(` Squares Degrees of Mean Sq`);\n \+ lprint(`Source (SS) Freedom (MS) F-R\\\n ati o p-value`);\n lprint(`----------------------------------------- -----\\\n ---------------`);\n printf(\"Row(A) %10.4f %5d %10.4f %8.4f %8.4f\",\n SSA, a - 1, MSA, MSA/MSE, PVR); \n lprint();\n printf(\"Col(B) %10.4f %5d %10.4f %8.4f \+ %8.4f\",\n SSB, b - 1, MSB, MSB/MSE, PVC);\n lprint();\n \+ printf(\"Int(AB) %10.4f %5d %10.4f %8.4f %8.4f\",\n SSA B, (a - 1)*(b - 1), MSAB, MSAB/MSE, PVInt);\n lprint();\n printf (\"Error %10.4f %5d %10.4f\", SSE, a*b*(c - 1),\n MSE);\n lprint();\n printf(\"Total %10.4f %5d\", SSTO, a*b*c - 1); \n lprint();\n lprint(`----------------------------------------- -----\\\n ---------------`)\n\nend proc;\n\n\nAnova2s := proc(D ::listlist)\nlocal a, b, tot, SSTO, SSA, SSB, SSE, RowMean, ColMean, i , j,\nGrandMean, MSA, MSB, MSE, PVR, PVC;\noption `July 5, 1993 -- Zav en A. Karian`;\ndescription \"Two factor ANOVA (one obs. per cell)\"; \n a := nops(D);\n b := nops(D[1]);\n tot := 0;\n SSTO := \+ 0;\n SSA := 0;\n SSB := 0;\n SSE := 0;\n RowMean := array( 1 .. a);\n ColMean := array(1 .. b);\n for j to b do ColMean[j] \+ := 0 od;\n for i to a do for j to b do tot := tot + D[i][j] od od; \n GrandMean := tot/(a*b);\n for i to a do RowMean[i] := Mean(D[ i]) od;\n for j to b do\n for i to a do ColMean[j] := ColMea n[j] + D[i][j] od;\n ColMean[j] := ColMean[j]/a\n od;\n f or i to a do for j to b do\n SSTO := SSTO + (D[i][j] - Gran dMean)^2\n od\n od;\n for i to a do SSA := SSA + (RowMean [i] - GrandMean)^2 od;\n for j to b do SSB := SSB + (ColMean[j] - G randMean)^2 od;\n SSA := b*SSA;\n SSB := a*SSB;\n SSE := SSTO - SSA - SSB;\n MSA := SSA/(a - 1);\n MSB := SSB/(b - 1);\n M SE := SSE/((a - 1)*(b - 1));\n PVR := 1 - FCDF(a - 1, (a - 1)*(b - \+ 1), MSA/MSE);\n PVC := 1 - FCDF(b - 1, (a - 1)*(b - 1), MSB/MSE);\n lprint(` Sum of`);\n lprint(` Squares Deg rees of Mean Sq`);\n lprint(`Source (SS) Freedom (MS ) F-R\\\n atio p-value`);\n lprint(`----------------- -----------------------------\\\n ---------------`);\n print f(\"Row(A) %10.4f %5d %10.4f %8.4f %8.4f\",\n SSA, a - 1, MSA, MSA/MSE, PVR);\n lprint();\n printf(\"Col(B) %10.4f \+ %5d %10.4f %8.4f %8.4f\",\n SSB, b - 1, MSB, MSB/MSE, PVC );\n lprint();\n printf(\"Error %10.4f %5d %10.4f\", SSE, \n (a - 1)*(b - 1), MSE);\n lprint();\n printf(\"Total \+ %10.4f %5d\", SSTO, a*b - 1);\n lprint();\n lprint(`---------- ------------------------------------\\\n ---------------`)\n\ne nd proc;\n\n\nAnovaBeta := proc()\nlocal n, k, xbar, xlist, ylist, alp hahat, betahat, sigma2hat,\nbetahat1, betahat2, ssto, ssr, sse, msr, m se, f, p;\noption `Mike, Andrew, and Christine were here.`;\ndescripti on \"Test beta = 0.\";\n if nargs < 1 or 2 < nargs then\n ER ROR(\"Must have list of lists, or (list,list)\")\n fi;\n n := no ps(args[1]);\n if nargs = 1 then\n xlist := [seq(args[1][k][ 1], k = 1 .. n)];\n ylist := [seq(args[1][k][2], k = 1 .. n)]\n else xlist := args[1]; ylist := args[2]\n fi;\n xbar := Mean (xlist);\n alphahat := Mean(ylist);\n k := 'k';\n betahat1 := sum(ylist[k]*(xlist[k] - xbar), k = 1 .. n);\n betahat2 := sum((xl ist[k] - xbar)^2, k = 1 .. n);\n betahat := betahat1/betahat2;\n \+ ssto := sum((ylist[k] - alphahat)^2, k = 1 .. n);\n ssr := sum(bet ahat^2*(xlist[k] - xbar)^2, k = 1 .. n);\n sse := ssto - ssr;\n \+ msr := ssr;\n mse := sse/(n - 2);\n f := msr/mse;\n p := 1 - \+ FCDF(1, n - 2, f);\n lprint(`-------------------------------------- --------\\\n -----------------`);\n lprint(` Su m of`);\n lprint(` Squares Degr of Mean Sq`);\n \+ lprint(` Source (SS) Freedom (MS) F\\\n -Ra tio p-value`);\n lprint(`---------------------------------------- ------\\\n -----------------`);\n printf(\"Regression %10.4f %5.0f %13.4f %10.4f %8.5f\",\n ssr, 1, msr, f, p);\n lprin t();\n printf(\" Error %13.4f %5.0f %13.4f\", sse, n - 2, mse);\n \+ lprint();\n printf(\" Total %13.4f %5.0f\", ssto, n - 1);\n \+ lprint();\n lprint(`---------------------------------------------- \\\n -----------------`)\n\nend proc;\n\n\nMedianTest := proc(X ::list(realcons), Y::list(realcons))\nlocal S, x, y, i, k;\noption\n`J oshua D. Levy ~ June 17, 1994, modified 6/1/96-zak`;\ndescription \" Median test\";\n if nargs = 2 then\n if type(nops(X) + nops( Y), even) then\n k := 1/2*nops(X) + 1/2*nops(Y)\n el se ERROR(\"sum of numbers of observations not even\")\n fi\n \+ elif type(args[3], posint) then k := args[3]\n else ERROR(\"bad ar gurment\")\n fi;\n S := sort([op(map(x, evalf(X))), op(map(y, ev alf(Y)))],\n (a, b) -> evalb(op(a) < op(b)));\n x := 1;\n \+ y := 0;\n S := eval(S);\n sum(S[i], i = 1 .. k)\n\nend proc;\n \n\nRunTest := proc(X::list(realcons), Y::\{realcons, list(realcons)\} )\nlocal S, x, y, i, runs;\noption\n`Joshua D. Levy ~ June 17, 1994, modified 6/1/96-zak`;\ndescription \"Run test. Number of runs in sort ed, combined l\\\nist of X and Y -or- number of runs above or below gi ven nu\\\nmber\";\n if type(Y, realcons) then\n S := map(una pply('sign'(x - evalf(Y)), x), X)\n else\n S := sort(\n \+ [op(map(x, evalf(X))), op(map(y, evalf(Y)))],\n (a, \+ b) -> evalb(op(a) < op(b)));\n x := 0;\n y := 1;\n \+ S := eval(S)\n fi;\n runs := 1;\n for i to nops(S) - 1 do\n if S[i] <> S[i + 1] then runs := runs + 1 fi\n od;\n run s\n\nend proc;\n\n\nSignTest := proc(C::boolean)\nlocal i, k, vars;\no ption\n`Joshua D. Levy ~ June 17, 1994, modified 6/1/96-zak`;\ndescr iption\n\"Number of elements in list(s) for which C is true\";\n if not type([args[2 .. nargs]], list(name = list)) then\n ERROR( \"all arguments beyond first must be of the f\\\n orm name=list \")\n fi;\n vars := select(type, indets(C), name);\n if vars \+ = \{\} then\n ERROR(\"test expression contains no variable\")\n fi;\n if vars <> indets([args[2 .. nargs]]) then ERROR(\"test\\ \n expression variable(s) must be defined in argumen\\\n \+ ts\")\n fi;\n subs(true = 1, false = 0, convert([seq(evalb(eva lf(subs(\n seq(lhs(args[i]) = rhs(args[i])[k], i = 2 .. nargs), \n C))), k = 1 .. nops(rhs(args[2])))], `+`))\n\nend proc;\n\n \nWilcoxon := proc(L::list(realcons), m::realcons)\nlocal S, n, w, x, \+ i, j, k;\noption `Joshua D. Levy ~ June 17, 1994`;\ndescription \"Wi lcoxon test for one median.\";\n S := sort([seq(evalf(x - m), x = L )],\n (a, b) -> evalb(abs(a) < abs(b)));\n for i while S[i] \+ = 0 do od;\n S := S[i .. nops(S)];\n n := nops(S);\n w := 0; \n for i to n do\n for j from i to n - 1 while abs(S[j]) = a bs(S[j + 1])\n do\n od;\n if j = i then w := w + \+ sign(S[i])*i\n else\n w :=\n w + 1/2* (i + j)*sum('sign(S[k])', k = i .. j)\n ;\n \+ i := j\n fi\n od;\n w\n\nend proc;\n\n\nWilcoxon2 := proc (X::list(realcons), Y::list(realcons))\nlocal S, n, x, y, w, i, j, k; \noption\n`Joshua D. Levy ~ June 17, 1994, modified 6/1/96-zak`;\nde scription \"Wilcoxon test for equality of two medians. Ret\\\nurns sum of ranks of Y in the sorted combination of X and \\\nY\";\n S := s ort([op(map(x, evalf(X))), op(map(y, evalf(Y)))],\n (a, b) -> e valb(op(a) < op(b)));\n n := nops(S);\n w := 0;\n for i to n \+ do\n for j from i to n - 1 while op(S[j]) = op(S[j + 1])\n \+ do\n od;\n if j = i then\n w := w + eval(s ubs(x = 0, y = 1, S[i]))*i\n else\n w := w + 1/2*(i \+ + j)*eval(\n subs(x = 0, y = 1, sum(S[k], k = i .. j))) ;\n i := j\n fi\n od;\n w\n\nend proc;\n\n\nCo ntingency := proc(L::list(list(nonnegint)))\nlocal h, k, n, npr, npc, \+ q, c, i, j;\noption `Joshua D. Levy ~ June 21, 1994`;\ndescription \+ \"Contingency table\";\n k := nops(L);\n h := nops(L[1]);\n n := sum('convert(L[i], `+`)', i = 1 .. k);\n npr := [seq(sum(L[i][j ], j = 1 .. h), i = 1 .. k)];\n npc := [seq(sum(L['i'][j], 'i' = 1 \+ .. k), j = 1 .. h)];\n q := 0;\n lprint();\n for i to k do\n \+ for j to h do printf(\"%7.0f \", L[i][j]) od;\n printf( \" %7.0f\\n \", npr[i]);\n for j to h do\n c := npr [i]*npc[j]/n;\n q := q + evalf((L[i][j] - c)^2/c);\n \+ printf(\"(%6.2f)\", c)\n od;\n lprint()\n od;\n \+ lprint();\n for j to h do printf(\"%7.0f \", npc[j]) od;\n pr intf(\" %7.0f\", n);\n lprint();\n lprint();\n printf(\"Chi-s quare statsistic (%.0f degrees of freedom\\\n ): %5g\\n\", (h - 1)*(k - 1), q);\n printf(\"p-value: %.6g\\n\",\n 1 - Chisqu areCDF((h - 1)*(k - 1), q));\n lprint()\n\nend proc;\n\n\nConvoluti on := proc(X1::list, X2::list)\nlocal dim1, dim2, l1, l2, j1, j2, i, L 1, L2, N1, N2, A, k,\nSUM, j, f1, f2, t;\noption `September 7, 1993 -- Zaven A. Karian`;\ndescription\n\"Computation of probabilities via th e convolution formula\";\n j1 := sum(X1[2*i], i = 1 .. 1/2*nops(X1) );\n j2 := sum(X2[2*i], i = 1 .. 1/2*nops(X2));\n if .1*10^(-6) \+ < abs(j1 - 1) or .1*10^(-6) < abs(j2 - 1)\n then ERROR(\"Probabilit ies must add up to 1\")\n fi;\n for i to 1/2*nops(X1) do\n \+ if X1[2*i - 1] < 0 then ERROR(\n \"Random variable values must be non-negative\")\n fi\n od;\n for i to 1/2*nops(X 2) do\n if X2[2*i - 1] < 0 then ERROR(\n \"Random va riable values must be non-negative\")\n fi\n od;\n dim1 : = nops(X1) + 2*X1[1];\n dim2 := nops(X2) + 2*X2[1];\n l1 := arra y(1 .. dim1);\n l2 := array(1 .. dim2);\n j1 := 1;\n j2 := 1; \n for i from 0 to X1[1] - 1 do\n l1[j1] := 0; l1[j1 + 1] := 0; j1 := j1 + 2\n od;\n for i to nops(X1) do l1[j1] := X1[i]; j 1 := j1 + 1 od;\n for i from 0 to X2[1] - 1 do\n l2[j2] := 0 ; l2[j2 + 1] := 0; j2 := j2 + 2\n od;\n for i to nops(X2) do l2[ j2] := X2[i]; j2 := j2 + 1 od;\n L1 := [seq(l1[i], i = 1 .. j1 - 1) ];\n L2 := [seq(l2[i], i = 1 .. j2 - 1)];\n N1 := 1/2*nops(L1); \n N2 := 1/2*nops(L2);\n A := array(1 .. 2*N1 + 2*N2 - 2);\n \+ for k from 0 to N1 + N2 - 2 do\n SUM := 0;\n for j from \+ 0 to k do\n if j < 0 or N1 - 1 < j then f1 := 0\n \+ else f1 := L1[2*j + 2]\n fi;\n if k - j < 0 or N2 - 1 < k - j then f2 := 0\n else f2 := L2[2*k - 2*j + 2] \n fi;\n SUM := SUM + f1*f2\n od;\n \+ A[2*k + 1] := k;\n A[2*k + 2] := SUM\n od;\n j := 0;\n \+ while A[j + 2] = 0 do j := j + 2 od;\n [seq(A[t], t = j + 1 .. 2 *N1 + 2*N2 - 2)]\n\nend proc;\n\n\nCraps := proc()\nlocal A, D, j, t, \+ point, roll, game;\noption `August 23, 1993 -- Zaven A. Karian`;\ndesc ription \"simulation of a single craps game\";\n A := array(1 .. 10 0);\n j := 2;\n D := Die(6, 2);\n point := D[1] + D[2];\n \+ A[j] := point;\n j := j + 1;\n if point = 7 or point = 11 then g ame := 1\n elif point = 2 or point = 3 or point = 12 then game := - 1\n else game := 0\n fi;\n while game = 0 do\n D := Di e(6, 2);\n roll := D[1] + D[2];\n A[j] := roll;\n \+ j := j + 1;\n if roll = point then game := 1\n elif rol l = 7 then game := -1\n fi\n od;\n if game = -1 then A[1] := 0 else A[1] := 1 fi;\n [seq(A[t], t = 1 .. j - 1)]\n\nend proc; \n\n\nFastCraps := proc(n::posint)\nlocal A, D, i, point, roll, game; \noption `August 23, 1993 -- Zaven A. Karian`;\ndescription \"simulati on of n craps games\";\n A := array(1 .. n);\n for i to n do\n \+ D := Die(6, 2);\n point := D[1] + D[2];\n if point = 7 or point = 11 then game := 1\n elif point = 2 or point = 3 or point = 12 then\n game := -1\n else game := 0\n \+ fi;\n while game = 0 do\n D := Die(6, 2);\n \+ roll := D[1] + D[2];\n if roll = point then game \+ := 1\n elif roll = 7 then game := -1\n fi\n \+ od;\n if game = -1 then A[i] := 0 else A[i] := 1 fi\n od; \n [seq(A[i], i = 1 .. n)]\n\nend proc;\n\n\nGraphRandWalk := proc( pn::numeric, ps::numeric, pe::numeric, steps::posint)\nlocal A, loc, r , j, L;\noption `July 5, 1993 -- Zaven A. Karian`;\ndescription \"Grap hic display of random walk\";\n A := array(1 .. 2*steps);\n loc \+ := array(1 .. 2);\n loc[1] := 0;\n loc[2] := 0;\n for j to st eps do\n r := rng();\n if r <= pn then loc[2] := loc[2] \+ + 1\n elif r <= pn + ps then loc[2] := loc[2] - 1\n elif r <= pn + ps + pe then loc[1] := loc[1] + 1\n else loc[1] := l oc[1] - 1\n fi;\n A[2*j - 1] := loc[1];\n A[2*j] \+ := loc[2]\n od;\n L := [seq([A[2*j - 1], A[2*j]], j = 1 .. steps )];\n plot(L)\n\nend proc;\n\n\nGraphRandWalk2 := proc(pn::numeric, pe::numeric, steps::posint)\nlocal A, loc, r1, r2, j, L;\noption `Jul y 5, 1993 -- Zaven A. Karian`;\ndescription \"Graphic display of rando m walk\";\n A := array(1 .. 2*steps);\n loc := array(1 .. 2);\n \+ loc[1] := 0;\n loc[2] := 0;\n for j to steps do\n r1 : = rng();\n r2 := rng();\n if r1 <= pn then loc[2] := loc [2] + 1\n else loc[2] := loc[2] - 1\n fi;\n if r2 <= pe then loc[1] := loc[1] + 1\n else loc[1] := loc[1] - 1\n \+ fi;\n A[2*j - 1] := loc[1];\n A[2*j] := loc[2]\n \+ od;\n L := [seq([A[2*j - 1], A[2*j]], j = 1 .. steps)];\n plo t(L)\n\nend proc;\n\n\nMarginalRelFreq := proc(A::listlist)\nlocal NA, rmin, rmax, cmin, cmax, i, r, c, RCOUNT, CCOUNT,\nrcount, ccount, RC, CC;\noption `July 5, 1993 -- Zaven A. Karian`;\ndescription\n\"Margin al relative frequencies of the list-of-lists, A\";\n NA := nops(A); \n rmin := A[1][2];\n rmax := A[1][2];\n cmin := A[1][1];\n \+ cmax := A[1][1];\n for i from 2 to NA do\n if A[i][2] < rm in then rmin := A[i][2] fi;\n if rmax < A[i][2] then rmax := A[ i][2] fi;\n if A[i][1] < cmin then cmin := A[i][1] fi;\n \+ if cmax < A[i][1] then cmax := A[i][1] fi\n od;\n if rmin < cmi n then cmin := rmin else rmin := cmin fi;\n if cmax < rmax then cma x := rmax else rmax := cmax fi;\n RCOUNT := array(rmin .. rmax);\n \+ CCOUNT := array(cmin .. cmax);\n for r from rmin to rmax do\n \+ rcount := 0;\n for i to NA do\n if A[i][2] = r \+ then rcount := rcount + 1 fi\n od;\n RCOUNT[r] := rcount \n od;\n for c from cmin to cmax do\n ccount := 0;\n \+ for i to NA do\n if A[i][1] = c then ccount := ccount + \+ 1 fi\n od;\n CCOUNT[c] := ccount\n od;\n RC := [se q(RCOUNT[r]/NA, r = rmin .. rmax)];\n CC := [seq(CCOUNT[c]/NA, c = \+ cmin .. cmax)];\n [CC, RC]\n\nend proc;\n\n\nRandWalk := proc(pn::n umeric, ps::numeric, pe::numeric,steps::posint, n::posint)\nlocal A, l oc, i, j, r;\noption `July 5, 1993 -- Zaven A. Karian`;\ndescription \+ \"Path of random walk\";\n A := array(1 .. n);\n loc := array(1 \+ .. 2);\n for i to n do\n loc[1] := 0;\n loc[2] := 0; \n for j to steps do\n r := rng();\n if r <= pn then loc[2] := loc[2] + 1\n elif r <= pn + ps then l oc[2] := loc[2] - 1\n elif r <= pn + ps + pe then loc[1] := loc[1] + 1\n else loc[1] := loc[1] - 1\n fi\n \+ od;\n A[i] := [loc[1], loc[2]]\n od;\n [seq(A[i], i = 1 .. n)]\n\nend proc;\n\n\nRandWalk2 := proc(pn::numeric, pe::numer ic, steps::posint, n::posint)\nlocal r1, r2, A, loc, j, i;\noption `Ju ly 5, 1993 -- Zaven A. Karian`;\ndescription \"Path of random walk\"; \n A := array(1 .. n);\n loc := array(1 .. 2);\n for i to n d o\n loc[1] := 0;\n loc[2] := 0;\n for j to steps \+ do\n r1 := rng();\n r2 := rng();\n if r1 <= pn then loc[2] := loc[2] + 1\n else loc[2] := loc[2] - 1\n fi;\n if r2 <= pe then loc[1] := loc[1] + 1\n else loc[1] := loc[1] - 1\n fi\n od; \n A[i] := [loc[1], loc[2]]\n od;\n [seq(A[i], i = 1 .. n )]\n\nend proc;\n\n\nConfIntAvLen := proc(LL::listlist)\nlocal n, Len, i;\noption `June 13, 1994 -- Zaven A. Karian`;\ndescription \"Average length of confidence intervals\";\n n := nops(LL);\n Len := 0; \n for i to n do Len := Len + LL[i][2] - LL[i][1] od;\n Len/n\n \nend proc;\n\n\nConfIntMean := proc()\nlocal LL, z, L, A, i, n, m, a, s, S, alpha, ConfLev, ptr;\noption `February 11, 1994 -- Zaven A. Kar ian`;\ndescription\n\"Compute end-points of confidence intervals for m u\";\n if not type(args[2], numeric) or args[2] < 0 or\n 100 < a rgs[2] then ERROR(\"Second argument must be a pe\\\n rcentage; \+ i.e. between 0 and 100\")\n fi;\n if nargs = 3 and args[3] <= 0 \+ then ERROR(\n \"Third argument, when present, must be positive \")\n fi;\n LL := args[1];\n ConfLev := args[2];\n if narg s = 3 then S := sqrt(args[3]) fi;\n n := nops(LL);\n m := nops(L L[1]);\n alpha := 1/2 + 1/200*ConfLev;\n if nargs = 3 or 24 < m \+ then z := NormalP(0, 1, alpha)\n else z := TP(m - 1, alpha)\n fi ;\n A := array(1 .. 2*n);\n ptr := 1;\n for i to n do\n \+ L := LL[i];\n a := Mean(L);\n if nargs = 2 then s := S tDev(L) else s := S fi;\n A[ptr] := evalf(a - z*s/sqrt(m));\n \+ A[ptr + 1] := evalf(a + z*s/sqrt(m));\n ptr := ptr + 2\n \+ od;\n [seq([A[2*i - 1], A[2*i]], i = 1 .. n)]\n\nend proc;\n\n\n ConfIntPlot := proc(LL::listlist, Value::realcons)\nlocal n, A, i;\nop tion `February 11, 1994 -- Zaven A. Karian`;\ndescription \"Graphic di splay of confidence intervals\";\n n := nops(LL);\n A := \{\};\n for i to n do\n A := A union \{[[LL[i][1], i], [LL[i][2], i ]]\}\n od;\n if nargs = 2 then A := A union \{[[Value, 1], [Valu e, n]]\}\n fi;\n plot(A, color = blue)\n\nend proc;\n\n\nConfInt Prop := proc(LL::listlist, ConfLev::realcons)\nlocal L, A, i, j, n, m, alpha, z, p, f, ptr;\noption `February 11, 1994 -- Zaven A. Karian`; \ndescription \"Compute confidence intervals for p\";\n if ConfLev \+ < 0 or 100 < ConfLev then ERROR(\"Second arg\\\n ument must be \+ a percentage; i.e. between 0 and 100\\\n \")\n fi;\n n := nops(LL);\n m := nops(LL[1]);\n alpha := 1/2 + 1/200*ConfLev;\n z := NormalP(0, 1, alpha);\n A := array(1 .. 2*n);\n ptr := \+ 1;\n for i to n do\n L := LL[i];\n p := sum(L[j], j = 1 .. m)/m;\n f := z*sqrt(p*(1 - p)/m);\n A[ptr] := eval f(p - f);\n A[ptr + 1] := evalf(p + f);\n ptr := ptr + 2 \n od;\n [seq([A[2*i - 1], A[2*i]], i = 1 .. n)]\n\nend proc;\n \n\nConfIntSuc := proc(L::listlist, Value::realcons)\nlocal n, count, \+ i, val, LL;\noption `June 13, 1994 -- Zaven A. Karian`;\ndescription\n \"number of confidence intervals that contain value\";\n val := eva lf(Value);\n LL := evalf(L);\n n := nops(LL);\n count := 0;\n for i to n do\n if val <= LL[i][2] and LL[i][1] <= val then \n count := count + 1\n fi\n od;\n count\n\nen d proc;\n\n\nConfIntVar := proc(LL::listlist, ConfLev::realcons)\nloca l L, A, i, n, m, ss, ptr, chi1, chi2, CHI1, CHI2;\noption `February 11 , 1994 -- Zaven A. Karian`;\ndescription \"Compute confidence interval s for sigma**2\";\n if ConfLev < 0 or 100 < ConfLev then ERROR(\"Se cond arg\\\n ument must be a percentage; i.e. between 0 and 100 \\\n \")\n fi;\n n := nops(LL);\n m := nops(LL[1]);\n \+ chi1 := 1/2 - 1/200*ConfLev;\n chi2 := 1 - chi1;\n CHI1 := (m - 1)/ChisquareP(m - 1, chi1);\n CHI2 := (m - 1)/ChisquareP(m - 1, \+ chi2);\n A := array(1 .. 2*n);\n ptr := 1;\n for i to n do\n \+ L := LL[i];\n ss := Variance(L);\n A[ptr] := ss*C HI2;\n A[ptr + 1] := ss*CHI1;\n ptr := ptr + 2\n od; \n [seq([A[2*i - 1], A[2*i]], i = 1 .. n)]\n\nend proc;\n\n\nCorrel ation := proc()\nlocal X, Y, N, SX, SY, XY, SumX, SumY, i, CXY;\noptio n `August 12, 1993 -- Zaven A. Karian`;\ndescription \"Mean correlatio n of the lists X and Y\";\n if nargs = 1 and not type(args[1], list list) then ERROR(\n \"When one argument is used, it must be a l ist of l\\\n ists\")\n fi;\n if nargs = 2 and\n not (t ype(args[1], list) and type(args[2], list)) then\n ERROR(\n \+ \"When two arguments are used, both must be lists\")\n fi;\n \+ if nargs = 2 and nops(args[1]) <> nops(args[2]) then\n ERROR( \"The two lists must be of the same length\")\n fi;\n if nargs = 1 then\n X := [seq(args[1][i][1], i = 1 .. nops(args[1]))];\n \+ Y := [seq(args[1][i][2], i = 1 .. nops(args[1]))]\n else X : = args[1]; Y := args[2]\n fi;\n N := nops(X);\n SX := StDev(X );\n SY := StDev(Y);\n XY := 0;\n SumX := 0;\n SumY := 0; \n for i to N do\n XY := XY + X[i]*Y[i];\n SumX := Su mX + X[i];\n SumY := SumY + Y[i]\n od;\n CXY := (N*XY - S umX*SumY)/(N*(N - 1));\n evalf(CXY/(SY*SX))\n\nend proc;\n\n\nLinRe g := proc()\nlocal X, Y, N, XY, SumX, SumY, i, CXY, slope, const, x;\n option `August 12, 1993 -- Zaven A. Karian`;\ndescription \"Slope and \+ intercept of regression line\";\n if nargs = 2 and not type(args[1] , listlist) then ERROR(\n \"When tow arguments are used the fir st must be a l\\\n ist of lists\")\n fi;\n if nargs = 3 a nd\n not (type(args[1], list) and type(args[2], list)) then\n \+ ERROR(\"When three arguments are used, the first tw\\\n o mus t be lists\")\n fi;\n if nargs = 3 and nops(args[1]) <> nops(arg s[2]) then\n ERROR(\"The two lists must be of the same length\" )\n fi;\n if nargs = 2 then\n X := [seq(args[1][i][1], i \+ = 1 .. nops(args[1]))];\n Y := [seq(args[1][i][2], i = 1 .. nop s(args[1]))];\n x := args[2]\n else X := args[1]; Y := args[ 2]; x := args[3]\n fi;\n N := nops(X);\n XY := 0;\n SumX : = 0;\n SumY := 0;\n for i to N do\n XY := XY + X[i]*Y[i]; \n SumX := SumX + X[i];\n SumY := SumY + Y[i]\n od;\n CXY := evalf((N*XY - SumX*SumY)/(N*(N - 1)));\n slope := CXY/Va riance(X);\n const := Mean(Y) - Mean(X)*slope;\n const + slope*x \n\nend proc;\n\n\nPlotPolyReg := proc()\nlocal X, Y, A, i, k, n, B, C , j, s, BI, eqn, x, pts, c, m;\noption `January 27, 1994 -- Zaven A. K arian`;\ndescription\n\"Produce scatter plot of X-Y pairs with poly. r eg. line\";\n if nargs = 2 and not type(args[1], listlist) then ERR OR(\n \"When two arguments are used, the first must be a \\\n \+ list of lists\")\n fi;\n if nargs = 3 and\n not (type(a rgs[1], list) and type(args[2], list)) then\n ERROR(\"When thre e arguments are used, the first tw\\\n o must be lists\")\n \+ fi;\n if nargs = 3 and nops(args[1]) <> nops(args[2]) then\n \+ ERROR(\"The two lists must have the same length\")\n fi;\n Digi ts := 50;\n if nargs = 2 then\n X := [seq(args[1][i][1], i = 1 .. nops(args[1]))];\n Y := [seq(args[1][i][2], i = 1 .. nops (args[1]))];\n m := args[2]\n else X := args[1]; Y := args[2 ]; m := args[3]\n fi;\n n := nops(X);\n B := array(1 .. m + 1 , 1 .. m + 1);\n C := array(1 .. m + 1, 1 .. 1);\n B[1, 1] := n; \n for i to m + 1 do for j from i to m + 1 do\n if j = 1 then j := j + 1 fi;\n s := 0;\n s := convert([s eq(X[k]^(i + j - 2), k = 1 .. n)],\n `+`);\n \+ B[i, j] := s;\n B[j, i] := s\n od\n od;\n C[1 , 1] := convert(Y, `+`);\n for i from 2 to m + 1 do\n s := 0 ;\n s := convert([seq(Y[k]*X[k]^(i - 1), k = 1 .. n)],\n \+ `+`);\n C[i, 1] := s\n od;\n BI := linalg[inverse](B );\n A := linalg[multiply](BI, C);\n eqn := 0;\n for i to m + 1 do eqn := eqn + A[i, 1]*x^(i - 1) od;\n pts := ScatPlot(X, Y);\n c := plot(eqn, x = Min(X) .. Max(X));\n with(plots, display);\n display(\{pts, c\})\n\nend proc;\n\n\nPolyReg := proc()\nlocal X, \+ Y, A, i, n, B, C, j, s, BI, eqn, x, k, m;\noption `January 27, 1994 -- Zaven A. Karian`;\ndescription\n\"Produce scatter plot of X-Y pairs w ith poly. reg. line\";\n if nargs = 3 and not type(args[1], listlis t) then ERROR(\n \"When three arguments are used, the first mus t be \\\n a list of lists\")\n fi;\n if nargs = 4 and\n \+ not (type(args[1], list) and type(args[2], list)) then\n ERRO R(\"When four arguments are used, the first two\\\n must be li sts\")\n fi;\n if nargs = 4 and nops(args[1]) <> nops(args[2]) t hen\n ERROR(\"The two lists must have the same length\")\n f i;\n Digits := 50;\n if nargs = 3 then\n X := [seq(args[1 ][i][1], i = 1 .. nops(args[1]))];\n Y := [seq(args[1][i][2], i = 1 .. nops(args[1]))];\n m := args[2];\n x := args[3] \n else\n X := args[1];\n Y := args[2];\n m := args[3];\n x := args[4]\n fi;\n n := nops(X);\n B := \+ array(1 .. m + 1, 1 .. m + 1);\n C := array(1 .. m + 1, 1 .. 1);\n \+ B[1, 1] := n;\n for i to m + 1 do for j from i to m + 1 do\n \+ if j = 1 then j := j + 1 fi;\n s := 0;\n \+ s := convert([seq(X[k]^(i + j - 2), k = 1 .. n)],\n `+ `);\n B[i, j] := s;\n B[j, i] := s\n od\n od;\n C[1, 1] := convert(Y, `+`);\n for i from 2 to m + 1 do \n s := 0;\n s := convert([seq(Y[k]*X[k]^(i - 1), k = 1 \+ .. n)],\n `+`);\n C[i, 1] := s\n od;\n BI := l inalg[inverse](B);\n A := linalg[multiply](BI, C);\n eqn := 0;\n for i to m + 1 do eqn := eqn + A[i, 1]*x^(i - 1) od;\n evalf(eq n, 10)\n\nend proc;\n\n\nRegAnal := proc()\nlocal n, k, xbar, xlist, y list, CL, c1, c2, c3, alphahat,\nbetahat, sigma2hat, betahat1, betahat 2, alphaconf, betaconf,\nsigma2conf;\noption `Mike, Andrew, and Christ ine were here.`;\ndescription \"Regression Analysis\";\n if nargs < = 1 or 3 < nargs then ERROR(\"Must have list \\\n of lists and \+ confidence level list, or (list,list,\\\n confidence list)\")\n fi;\n n := nops(args[1]);\n if nargs = 2 then\n xlist := [seq(args[1][k][1], k = 1 .. n)];\n ylist := [seq(args[1][k ][2], k = 1 .. n)];\n CL := args[2]\n else xlist := args[1]; ylist := args[2]; CL := args[3]\n fi;\n xbar := Mean(xlist);\n \+ alphahat := Mean(ylist);\n k := 'k';\n betahat1 := sum(ylist[ k]*(xlist[k] - xbar), k = 1 .. n);\n betahat2 := sum((xlist[k] - xb ar)^2, k = 1 .. n);\n betahat := betahat1/betahat2;\n sigma2hat \+ := sum(\n (ylist[k] - alphahat - betahat*(xlist[k] - xbar))^2, \n k = 1 .. n)/n;\n alphaconf[1] := alphahat -\n TP(n - 2, 1 - 1/2*CL[1])*sqrt(1.0*sigma2hat/(n - 2));\n alphaconf[2] := alphahat +\n TP(n - 2, 1 - 1/2*CL[1])*sqrt(1.0*sigma2hat/(n - \+ 2));\n betaconf[1] := betahat - TP(n - 2, 1 - 1/2*CL[2])*\n \+ sqrt(1.0*n*sigma2hat/((n - 2)*betahat2));\n betaconf[2] := betahat \+ + TP(n - 2, 1 - 1/2*CL[2])*\n sqrt(1.0*n*sigma2hat/((n - 2)*bet ahat2));\n sigma2conf[1] :=\n n*sigma2hat/ChisquareP(n - 2, \+ 1 - 1/2*CL[3]);\n sigma2conf[2] := n*sigma2hat/ChisquareP(n - 2, 1/ 2*CL[3])\n ;\n c1 := 1 - CL[1];\n c2 := 1 - CL[2];\n c 3 := 1 - CL[3];\n lprint(`----------------------------------------- -----\\\n ---------------`);\n lprint(` Poin t Confidence C\\\n onfidence`);\n lprint(`Paramet er Estimates Level \\\n Interval`);\n lp rint(`----------------------------------------------\\\n ------ ---------`);\n printf(\"alpha %13.4f %10.2f [%1.4f,%1.4 f]\\n\",\n alphahat, c1, alphaconf[1], alphaconf[2]);\n lpri nt();\n printf(\"beta %13.4f %10.2f [%1.4f,%1.4f]\\n\" ,\n betahat, c2, betaconf[1], betaconf[2]);\n lprint();\n \+ printf(\"sigma^2 %13.4f %10.2f [%1.4f,%1.4f]\\n\",\n \+ sigma2hat, c3, sigma2conf[1], sigma2conf[2]);\n lprint(`---------- ------------------------------------\\\n ---------------`)\n\ne nd proc;\n\n\nRegBand :=proc()\nlocal X, Y, N, A, Minx, Maxx, Miny, Ma xy, XR, YR, EX, EY, i,\nx, MX, alpha, beta, Sig2, SSQ, T, Option, Conf Lev, Term1,\nTerm2, upper, lower, p1, p2, p3, p4, L, Ymin, Ymax;\nopti on `August 21, 1993 -- Zaven A. Karian`;\ndescription\n\"Produce scatt er plot of X-Y pairs with the reg. line\";\n if nargs = 3 and not t ype(args[1], listlist) then ERROR(\n \"When one argument is use d, it must be a list of l\\\n ists\")\n fi;\n if nargs = \+ 4 and\n not (type(args[1], list) and type(args[2], list)) then\n \+ ERROR(\n \"When two arguments are used, both must be lists \")\n fi;\n if nargs = 4 and nops(args[1]) <> nops(args[2]) then \n ERROR(\"The two lists must have the same length\")\n fi; \n if nargs = 3 then\n X := [seq(args[1][i][1], i = 1 .. nop s(args[1]))];\n Y := [seq(args[1][i][2], i = 1 .. nops(args[1]) )];\n ConfLev := args[2];\n Option := args[3]\n else \n X := args[1];\n Y := args[2];\n ConfLev := arg s[3];\n Option := args[4]\n fi;\n N := nops(X);\n A := array(1 .. 2*N);\n MX := Mean(X);\n for i to N do A[2*i - 1] := X[i]; A[2*i] := Y[i] od;\n Minx := Min(X);\n Maxx := Max(X);\n \+ Miny := Min(Y);\n Maxy := Max(Y);\n EX := 1/20*Maxx - 1/20*Mi nx;\n EY := 1/20*Maxy - 1/20*Miny;\n XR := Minx - EX .. Maxx + E X;\n L := LinReg(X, Y, x);\n beta := op(1, op(2, L));\n alpha := op(1, L) + beta*MX;\n Sig2 := 0;\n for i to N do\n Si g2 := Sig2 + (Y[i] - alpha - beta*(X[i] - MX))^2\n od;\n Sig2 := Sig2/N;\n SSQ := (N - 1)*Variance(X);\n if N < 33 then T := TP( N - 2, .01*ConfLev)\n else T := NormalP(0, 1, .01*ConfLev)\n fi; \n Term1 := alpha + beta*(x - MX);\n if Option = \"C\" then Term 2 :=\n T*sqrt(N*Sig2*(1/N + (x - MX)^2/SSQ)/(N - 2))\n else \+ Term2 :=\n T*sqrt(N*Sig2*(1 + 1/N + (x - MX)^2/SSQ)/(N - 2))\n \+ fi;\n upper := Term1 + Term2;\n lower := Term1 - Term2;\n \+ Ymin := min(Miny - EY, subs(x = Minx, lower),\n subs(x = Maxx, \+ lower));\n Ymax := max(Maxy + EY, subs(x = Minx, upper),\n s ubs(x = Maxx, upper));\n YR := Ymin .. Ymax;\n p1 := plot([seq([ A[2*i - 1], A[2*i]], i = 1 .. N)], XR,\n YR, style = POINT, col or = black);\n p2 := plot(L, x = XR, YR, color = red);\n p3 := p lot(upper, x = XR, YR, color = blue);\n p4 := plot(lower, x = XR, Y R, color = blue);\n with(plots, display);\n display(\{p1, p2, p3 , p4\})\n\nend proc;\n\n\nResiduals := proc()\nlocal X, Y, N, Minx, Ma xx, XR, YR, p2, P, EX, EY, i, x, L,\nalpha, beta, Res, p, f, m;\noptio n `August 21, 1993 -- Zaven A. Karian`;\ndescription\n\"Produce scatte r plot of X-Y pairs with the reg. line\";\n if nargs = 2 and not ty pe(args[1], listlist) then ERROR(\n \"When two arguments are us ed, the first two must b\\\n e a list of lists\")\n fi;\n \+ if nargs = 3 and\n not (type(args[1], list) and type(args[2], list )) then\n ERROR(\"When three arguments are used, the first tw\\ \n o must be lists\")\n fi;\n if nargs = 3 and nops(args[ 1]) <> nops(args[2]) then\n ERROR(\"The two lists must have the same length\")\n fi;\n if nargs = 2 then\n X := [seq(arg s[1][i][1], i = 1 .. nops(args[1]))];\n Y := [seq(args[1][i][2] , i = 1 .. nops(args[1]))];\n m := args[2]\n else X := args[ 1]; Y := args[2]; m := args[3]\n fi;\n N := nops(X);\n if m = 1 then L := LinReg(X, Y, x)\n else L := PolyReg(X, Y, m, x)\n f i;\n f := student[makeproc](L, x);\n Res := [seq(Y[i] - f(X[i]), i = 1 .. N)];\n Minx := Min(X);\n Maxx := Max(X);\n EX := 1/ 20*Maxx - 1/20*Minx;\n XR := Minx - EX .. Maxx + EX;\n EY := .05 000000000*Max(Res) - .05000000000*Min(Res);\n YR := Min(Res) - EY . . Max(Res) + EY;\n P := \{\};\n for i to N do P := P union \{[X[ i], Res[i]]\} od;\n p2 := plot(P, XR, YR, style = POINT, color = re d);\n print(p2);\n lprint(` Independent`);\n lpri nt(` Variable Residual`);\n lprint(` \+ -----------------------`);\n for i to N do\n printf(\" %20.4 f %13.4f\", X[i], Res[i]); lprint()\n od\n\nend proc;\n\n\nScatPlot := proc()\nlocal X, Y, N, A, Minx, Maxx, Miny, Maxy, XR, YR, EX, EY, \+ i;\noption `January 27, 1994 -- Zaven A. Karian`;\ndescription \"Produ ce a scatter plot of X-Y pairs\";\n if nargs = 1 and not type(args[ 1], listlist) then ERROR(\n \"When one argument is used it must be a list of li\\\n sts\")\n fi;\n if nargs = 2 and\n \+ not (type(args[1], list) and type(args[2], list)) then\n ERROR (\n \"When two arguments are used, both must be lists\")\n f i;\n if nargs = 2 and nops(args[1]) <> nops(args[2]) then\n \+ ERROR(\"The two lists must have the same length\")\n fi;\n if na rgs = 1 then\n X := [seq(args[1][i][1], i = 1 .. nops(args[1])) ];\n Y := [seq(args[1][i][2], i = 1 .. nops(args[1]))]\n els e X := args[1]; Y := args[2]\n fi;\n N := nops(X);\n A := arr ay(1 .. 2*N);\n for i to N do A[2*i - 1] := X[i]; A[2*i] := Y[i] od ;\n Minx := Min(X);\n Maxx := Max(X);\n Miny := Min(Y);\n \+ Maxy := Max(Y);\n EX := 1/20*Maxx - 1/20*Minx;\n EY := 1/20*Maxy - 1/20*Miny;\n XR := Minx - EX .. Maxx + EX;\n YR := Miny - EY \+ .. Maxy + EY;\n plot([seq([A[2*i - 1], A[2*i]], i = 1 .. N)], XR, Y R,\n style = POINT)\n\nend proc;\n\n\nScatPlotLine, proc()\nloc al X, Y, N, A, Minx, Maxx, Miny, Maxy, XR, YR, EX, EY, i,\nx, p1, p2, \+ L;\noption `August 21, 1993 -- Zaven A. Karian`;\ndescription\n\"Produ ce scatter plot of X-Y pairs with the reg. line\";\n if nargs = 1 a nd not type(args[1], listlist) then ERROR(\n \"When one argumen t is used, it must be a list of l\\\n ists\")\n fi;\n if \+ nargs = 2 and\n not (type(args[1], list) and type(args[2], list)) t hen\n ERROR(\n \"When two arguments are used, both must \+ be lists\")\n fi;\n if nargs = 2 and nops(args[1]) <> nops(args[ 2]) then\n ERROR(\"The two lists must have the same length\")\n fi;\n if nargs = 1 then\n X := [seq(args[1][i][1], i = 1 .. nops(args[1]))];\n Y := [seq(args[1][i][2], i = 1 .. nops(a rgs[1]))]\n else X := args[1]; Y := args[2]\n fi;\n N := nops (X);\n A := array(1 .. 2*N);\n for i to N do A[2*i - 1] := X[i]; A[2*i] := Y[i] od;\n Minx := Min(X);\n Maxx := Max(X);\n Min y := Min(Y);\n Maxy := Max(Y);\n EX := 1/20*Maxx - 1/20*Minx;\n \+ EY := 1/20*Maxy - 1/20*Miny;\n XR := Minx - EX .. Maxx + EX;\n \+ YR := Miny - EY .. Maxy + EY;\n L := LinReg(X, Y, x);\n p1 := \+ plot([seq([A[2*i - 1], A[2*i]], i = 1 .. N)], XR,\n YR, style = POINT);\n p2 := plot(L, x = XR, YR);\n with(plots, display);\n \+ display(\{p1, p2\})\n\nend proc;\n\n\nBoxWhisker := proc()\nlocal N , i, A, D, m, M, Minn, Maxx, P25, P50, P75, eh;\noption `August 19, 19 93 -- Zaven A. Karian`;\ndescription \"Box-and whisker display\";\n \+ N := nargs;\n Minn := Min(args[1]);\n Maxx := Max(args[1]);\n \+ for i to N do\n if not type(args[i], list) then\n \+ ERROR(\"All arguments must be lists\")\n fi\n od;\n A := \+ \{\};\n for i to N do\n D := args[i];\n m := Min(D); \n M := Max(D);\n if m < Minn then Minn := m fi;\n \+ if Maxx < M then Maxx := M fi;\n P25 := Percentile(D, .25);\n P50 := Percentile(D, .50);\n P75 := Percentile(D, .75); \n A := A union \{[[P75, i - .4], [P75, i + .4]],\n \+ [[P25, i + .4], [P75, i + .4]],\n [[P25, i - .4], [P25, i + .4]],\n [[P50, i - .4], [P50, i + .4]],\n [[m, \+ i], [P25, i]], [[P75, i], [M, i]],\n [[P25, i - .4], [P75, \+ i - .4]]\}\n od;\n eh := 1/20*Maxx - 1/20*Minn;\n plot(A, Min n - eh .. Maxx + eh, 0 .. N + 1, color = blue)\n\nend proc;\n\n\nHisto gram := proc()\nlocal YYY, YY, Y, n, l, u, v, i, j, k, p, t, U, V, cou nt,\ntruemin, truemax, nint;\noption `July 7, 1993 -- Zaven A. Karian` ;\ndescription \"Graphic display of histogram\";\n YY := sort(evalf (args[1]));\n n := nops(YY);\n if nargs = 1 then\n truemi n := YY[1]; truemax := YY[n]; nint := 8\n else\n truemin := \+ lhs(args[2]);\n truemax := rhs(args[2]);\n nint := args[ 3]\n fi;\n count := 0;\n for i to n do\n if YY[i] < tr uemin or truemax < YY[i] then\n count := count + 1\n \+ fi\n od;\n if 0 < count then print(\"WARNING: There are\", coun t,\n \" data points out of plot range.\")\n fi;\n j := 1; \n YYY := array(1 .. n);\n u := array(1 .. nint + 1);\n v := \+ array(1 .. nint + 1);\n for i to n do\n if truemin <= YY[i] \+ and YY[i] <= truemax then\n YYY[j] := YY[i]; j := j + 1\n \+ fi\n od;\n n := j - 1;\n Y := [seq(YYY[k], k = 1 .. n)] ;\n p := array(1 .. 8*nint);\n l := (truemax - truemin)/nint;\n \+ u[1] := truemin;\n v[1] := 0;\n for i from 2 to nint + 1 do\n u[i] := u[i - 1] + l; v[i] := 0\n od;\n if u[nint + 1] < Y[n] then u[nint + 1] := Y[n] fi;\n i := 2;\n for j to n do\n \+ if Y[j] <= u[i] then v[i - 1] := v[i - 1] + 1\n else i := i + 1; j := j - 1\n fi\n od;\n t := 1;\n for i to nin t do\n U[t] := u[i];\n V[t] := 0;\n t := t + 1;\n U[t] := U[t - 1];\n V[t] := v[i];\n t := t + 1; \n U[t] := U[t - 1] + l;\n V[t] := v[i];\n t := t + 1;\n U[t] := U[t - 1];\n V[t] := 0;\n t := t + 1\n od;\n for i to t - 1 do\n p[2*i - 1] := U[i]; p[2*i] := V[i]/(l*n)\n od;\n plot([seq([p[2*i - 1], p[2*i]], i = 1 .. \+ t - 1)],\n color = red, linestyle = 3)\n\nend proc;\n\n\nOgive, proc(L::list, R::range, n::posint)\nlocal Min, Max, freq, A, X, LEP, \+ REP, j, inc, i, N;\noption `July 5, 1993 -- Zaven A. Karian`;\ndescrip tion \"Graphic display of ogive for L\";\n Min := lhs(R);\n Max \+ := rhs(R);\n freq := array(1 .. n);\n A := array(1 .. 2*n + 2); \n N := nops(L);\n X := sort(L);\n LEP := Min;\n inc := (M ax - Min)/n;\n REP := LEP + inc;\n for j to n do freq[j] := 0 od ;\n for j to n do\n for i to N do\n if X[i] < REP then freq[j] := freq[j] + 1 fi\n od;\n LEP := REP;\n \+ REP := LEP + inc\n od;\n A[1] := Min;\n A[2] := 0;\n \+ for j to n do\n A[2*j + 1] := A[2*j - 1] + inc;\n A[2*j \+ + 2] := freq[j]/N\n od;\n [seq([A[2*i - 1], A[2*i]], i = 1 .. n \+ + 1)];\n plot(%, color = blue, linestyle = 3)\n\nend proc;\n\n\nPlo tDiscCDF := proc(expr::algebraic, R::range)\nlocal A, MIN, p, i, var, \+ Min, Max, newp;\noption `August 8, 1993 -- Zaven A. Karian`;\ndescript ion \"Graphic display of distribution expr\";\n Min := lhs(R);\n \+ Max := rhs(R);\n A := \{\};\n MIN := Min;\n p := 0;\n var := op(select(type, indets(expr), name));\n for i to Max - Min + 1 \+ do\n newp := eval(subs(var = MIN, expr));\n if p < 1 and 0 <= newp then\n p := p + eval(subs(var = MIN, expr))\n \+ fi;\n A := A union \{[[MIN, p], [MIN + 1, p]]\};\n \+ MIN := MIN + 1\n od;\n plot(A, color = blue)\n\nend proc;\n\n\nP lotEmpCDF := proc(L::list, R::range)\nlocal i, A, N, X;\noption `Augus t 8, 1993 -- Zaven A. Karian`;\ndescription \"Graphic display of empir ical distribution L\";\n N := nops(L);\n X := sort(L);\n A := \{[[min(lhs(R), X[1]), 0], [X[1], 0]]\};\n for i to N - 1 do\n \+ if X[i] < X[i + 1] then\n A := A union \{[[X[i], i/N], \+ [X[i + 1], i/N]]\}\n fi\n od;\n A := A union \{[[X[N], 1] , [max(rhs(R), X[N]), 1]]\};\n plot(A, color = red, linestyle = 3) \n\nend proc;\n\n\nPlotEmpPDF := proc()\nlocal YYY, YY, Y, n, l, u, v, i, j, k, p, t, U, V, count,\ntruemin, truemax, nint;\noption `June 1, 1996 -- Zaven A. Karian`;\ndescription \"Graphic display of histogram \";\n YY := sort(evalf(args[1]));\n n := nops(YY);\n if nargs = 1 then\n truemin := YY[1]; truemax := YY[n]; nint := 8\n \+ else\n truemin := lhs(args[2]);\n truemax := rhs(args[2] );\n nint := args[3]\n fi;\n count := 0;\n for i to n \+ do\n if YY[i] < truemin or truemax < YY[i] then\n co unt := count + 1\n fi\n od;\n if 0 < count then print(\"W ARNING: There are\", count,\n \" data points out of plot range. \")\n fi;\n j := 1;\n YYY := array(1 .. n);\n u := array(1 .. nint + 1);\n v := array(1 .. nint + 1);\n for i to n do\n \+ if truemin <= YY[i] and YY[i] <= truemax then\n YYY[j] := YY[i]; j := j + 1\n fi\n od;\n n := j - 1;\n Y := \+ [seq(YYY[k], k = 1 .. n)];\n p := array(1 .. 2*nint + 4);\n l := (truemax - truemin)/nint;\n u[1] := truemin;\n v[1] := 0;\n \+ for i from 2 to nint + 1 do\n u[i] := u[i - 1] + l; v[i] := 0\n od;\n if u[nint + 1] < Y[n] then u[nint + 1] := Y[n] fi;\n i := 2;\n for j to n do\n if Y[j] <= u[i] then v[i - 1] := v[ i - 1] + 1\n else i := i + 1; j := j - 1\n fi\n od;\n U[1] := truemin - 1/2*l;\n V[1] := 0;\n t := 2;\n for i t o nint do\n U[t] := u[i] + 1/2*l; V[t] := v[i]; t := t + 1\n \+ od;\n U[nint + 2] := truemax + 1/2*l;\n V[nint + 2] := 0;\n \+ for i to nint + 2 do\n p[2*i - 1] := U[i]; p[2*i] := V[i]/(l*n) \n od;\n [seq([p[2*i - 1], p[2*i]], i = 1 .. nint + 2)];\n pl ot(%, color = red, linestyle = 3)\n\nend proc;\n\n\nPlotRunningAverage := proc(L::list)\nlocal N, SUM, A, i;\noption `August 23, 1993 -- Zav en A. Karian`;\ndescription \"Plotting the running average of a list\" ;\n N := nops(L);\n SUM := L[1];\n A := array(1 .. 2*N);\n \+ A[1] := 1;\n A[2] := L[1];\n for i from 2 to N do\n SUM \+ := SUM + L[i]; A[2*i - 1] := i; A[2*i] := SUM/i\n od;\n plot([se q([A[2*i - 1], A[2*i]], i = 1 .. N)], color = red)\n\nend proc;\n\n\nP robHist := proc()\nlocal n, L, lp, ap, A, MIN, dim, i, p, var, Min, Ma x;\noption `January 26, 1994 -- Zaven A. Karian`;\ndescription \"Graph ic display of probability histogram\";\n if nargs = 1 and not type( args[1], list) then ERROR(\"I\\\n f one argument is used, it mu st be a probability l\\\n ist\")\n fi;\n if nargs = 2 and not type(args[1], algebraic) then ERROR(\n \"If two arguments \+ are used, the first must be an e\\\n xpression\")\n fi;\n \+ if nargs = 2 and not type(args[2], range) then ERROR(\"\\\n If two arguments are used, the second must be a ra\\\n nge\")\n \+ fi;\n if nargs = 1 then\n n := nops(args[1]);\n L : = copy(args[1]);\n A := array(1 .. 4*n);\n lp := 1;\n \+ ap := 1;\n for i to 1/2*n do\n A[ap] := L[lp] - .5;\n A[ap + 1] := 0;\n A[ap + 2] := A[ap];\n \+ A[ap + 3] := L[lp + 1];\n A[ap + 4] := L[lp] + .5 ;\n A[ap + 5] := L[lp + 1];\n A[ap + 6] := A[ap \+ + 4];\n A[ap + 7] := 0;\n lp := lp + 2;\n \+ ap := ap + 8\n od;\n plot([seq([A[2*i - 1], A[2*i]] , i = 1 .. 2*n)],\n color = blue)\n else\n Min := lhs(args[2]);\n Max := rhs(args[2]);\n dim := 8*trunc(M ax - Min + 1/1000000000) + 8;\n A := array(1 .. dim);\n \+ MIN := Min;\n var := op(select(type, indets(args[1]), name));\n for i to Max - Min + 1 do\n p := evalf(subs(var = M IN, args[1]));\n A[8*i - 7] := MIN - .5;\n A[8*i - 6] := 0;\n A[8*i - 5] := MIN - .5;\n A[8*i - \+ 4] := p;\n A[8*i - 3] := MIN + .5;\n A[8*i - 2] \+ := p;\n A[8*i - 1] := MIN + .5;\n A[8*i] := 0;\n MIN := MIN + 1\n od;\n [seq([A[2*i - 1], A[2 *i]], i = 1 .. 4*Max - 4*Min + 4)\n ];\n plot(%, col or = blue)\n fi\n\nend proc;\n\n\nQQ := proc()\nlocal X, Y, N, M, X X, YY, MIN, P, NY, NEWY, NEWYL, i, r, f;\noption `August 12, 1993 -- Z aven A. Karian`;\ndescription \"Produce the Q-Q plot of the lists X an d Y\";\n if nargs = 1 and not type(args[1], listlist) then ERROR(\n \"When one argument is used, it must be a list of l\\\n \+ ists\")\n fi;\n if nargs = 2 and\n not (type(args[1], list) \+ and type(args[2], list)) then\n ERROR(\n \"When two argu ments are used, both must be lists\")\n fi;\n if nargs = 1 then \n X := [seq(args[1][i][1], i = 1 .. nops(args[1]))];\n \+ Y := [seq(args[1][i][2], i = 1 .. nops(args[1]))]\n else X := args[ 1]; Y := args[2]\n fi;\n N := nops(X);\n M := nops(Y);\n i f N <= M then XX := sort(X); YY := sort(Y)\n else XX := sort(Y); YY := sort(X)\n fi;\n MIN := min(N, M);\n NEWY := array(1 .. MI N);\n P := [seq(i/(MIN + 1), i = 1 .. MIN)];\n NY := nops(YY);\n for i to MIN do\n r := trunc(P[i]*(NY + 1));\n f := \+ P[i]*(NY + 1) - r;\n if f = 0 then NEWY[i] := YY[r]\n el se NEWY[i] := YY[r] + f*(YY[r + 1] - YY[r])\n fi\n od;\n \+ NEWYL := [seq(NEWY[i], i = 1 .. MIN)];\n ScatPlot(XX, NEWYL)\n\nend proc;\n\n\nQQFit := proc(L::list, expr::algebraic, R::range)\nlocal\n var, DataList, DistList, NQuantiles, Range, PlotRange, i;\n var := \+ op(select(type, indets(expr), name));\n if nops([var]) = 0 or 1 < n ops([var]) then ERROR(\"The \\\n second argument must have exac tly one varaible in \\\n it\")\n fi;\n if nargs = 2 then \n Range := -infinity .. infinity;\n PlotRange := Min(L) .. Max(L)\n elif nargs = 3 then Range := R; PlotRange := R\n el se ERROR(\n \"This procedure requires either 2 or 3 arguments. \")\n fi;\n NQuantiles := nops(L);\n DataList := sort(L);\n \+ DistList := [seq(\n fsolve(i/(NQuantiles + 1) = expr, var, Ra nge),\n i = 1 .. NQuantiles)];\n plot(\{\n plot(x, x \+ = PlotRange), ScatPlot(DistList, DataList)\})\n\nend proc;\n\n\nTimePl ot, proc()\nlocal N, B, n, i, j, t, S, A;\noption `August 19, 1993 -- \+ Zaven A. Karian`;\ndescription \"Time series plot of a list\";\n N \+ := nargs;\n n := nops(args[1]);\n S := \{\};\n for i to N do \n if not type(args[i], list) then\n ERROR(\"All arg uments must be lists\")\n fi;\n if n < nops(args[i]) the n n := nops(args[i]) fi\n od;\n B := array(1 .. 2*n);\n for i to N do\n A := args[i];\n for j to nops(A) do B[2*j - 1 ] := j; B[2*j] := A[j]\n od;\n S := S union\n \+ \{[seq([B[2*t - 1], B[2*t]], t = 1 .. nops(A))]\}\n od;\n plot( S)\n\nend proc;\n\n\nXbarChart := proc(X::list, SampleSize::posint, sb ar::numeric)\nlocal N, XX, S, Sbar, n, Xbarbar, A3, i, offset, UCL, LC L,\npts, P1, P2;\noption `January 24, 1994 -- Zaven A. Karian`;\ndescr iption \"Produce the control chart for the mean \";\n N := nops(X); \n if type(X, listlist) then\n XX := [seq(Mean(X[i]), i = 1 \+ .. N)];\n S := [seq(StDev(X[i]), i = 1 .. N)];\n Sbar := Mean(S);\n n := nops(X[1])\n else XX := copy(X); Sbar := sb ar; n := SampleSize\n fi;\n Xbarbar := Mean(XX);\n A3 := 3*sq rt(n - 1)*GAMMA(1/2*n - 1/2)/(\n sqrt(n)*sqrt(2)*GAMMA(1/2*n)); \n offset := A3*Sbar;\n UCL := [[0, Xbarbar + offset], [N + 1, X barbar + offset]]\n ;\n LCL := [[0, Xbarbar - offset], [N + \+ 1, Xbarbar - offset]]\n ;\n pts := [seq([i, XX[i]], i = 1 .. N)];\n P1 := plot(\{UCL, LCL\}, color = red);\n P2 := plot(pts, color = blue);\n with(plots, display);\n display(\{P1, P2\})\n \nend proc;\n\n\nBetaP := proc(b3::algebraic, b4::algebraic, p::algebr aic)\nlocal eqn, C, x, y;\noption `August 12, 1993 -- Zaven A. Karian` ;\ndescription \"Beta percentile\";\n if not (type(b3, numeric) and type(b4, numeric) and\n type(p, numeric)) then RETURN('procname(ar gs)')\n fi;\n if b3 <= -1 then ERROR(\"First argument must be >= 1\") fi;\n if b4 <= -1 then ERROR(\"Second argument must be >= 1\" ) fi\n ;\n if p < 0 or 1 < p then\n ERROR(\"Third argumen t must be between 0 and 1\")\n fi;\n C := GAMMA(b3 + b4 + 2)/(GA MMA(b3 + 1)*GAMMA(b4 + 1));\n eqn := p = int(C*x^b3*(1 - x)^b4, x = 0 .. y);\n fsolve(eqn, y, y = 0 .. 1)\n\nend proc;\n\n\nChisquareP := proc(r::algebraic, p::algebraic)\nlocal eqn, a, x, y;\noption `Aug ust 12, 1993 -- Zaven A. Karian`;\ndescription \"Chi square percentile \";\n if not (type(r, numeric) and type(p, numeric)) then\n \+ RETURN('procname(args)')\n fi;\n if p < 0 or 1 < p then\n \+ ERROR(\"Second argument must be between 0 and 1\")\n fi;\n a := 1/(GAMMA(1/2*r)*2^(1/2*r));\n eqn := p = int(a*x^(1/2*r - 1)*exp(- 1/2*x), x = 0 .. y);\n fsolve(eqn, y, y = 0 .. infinity)\n\nend pr oc;\n\n\nExponentialP := proc(t::algebraic, p::algebraic)\nlocal eqn, \+ x, y;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescription \"Exp onential percentile\";\n if not (type(t, numeric) and type(p, numer ic)) then\n RETURN('procname(args)')\n fi;\n if p < 0 or \+ 1 < p then\n ERROR(\"Second argument must be between 0 and 1\") \n fi;\n if t <= 0 then ERROR(\"First argument must be positive \")\n fi;\n eqn := p = int(exp(- x/t)/t, x = 0 .. y);\n fsolv e(eqn, y, y = 0 .. infinity)\n\nend proc;\n\n\nFP := proc(df1::algebra ic, df2::algebraic, p::algebraic)\nlocal eqn, num, den, x, y;\noption \+ `August 12, 1993 -- Zaven A. Karian`;\ndescription \"F percentile\";\n if not (type(df1, numeric) and type(df2, numeric) and\n type(p, numeric)) then RETURN('procname(args)')\n fi;\n if p < 0 or 1 < p then\n ERROR(\"Third argument must be between 0 and 1\")\n \+ fi;\n if df1 < 0 or df2 < 0 then\n ERROR(\"First two argum ents must be positive\")\n fi;\n num := x^(1/2*df1 - 1)*GAMMA(1/ 2*df1 + 1/2*df2)*\n (df1/df2)^(1/2*df1);\n den := GAMMA(1/2* df1)*GAMMA(1/2*df2)*\n (1 + df1*x/df2)^(1/2*df1 + 1/2*df2);\n \+ eqn := p = int(num/den, x = 0 .. y);\n fsolve(eqn, y, y = 0 .. in finity)\n\nend proc;\n\n\nGBP := proc(b1::algebraic, b2::algebraic, b3 ::algebraic,b4::algebraic, p::algebraic)\nlocal eqn, C, CC, x, y;\nopt ion `August 12, 1993 -- Zaven A. Karian`;\ndescription \"Generalized b eta percentile\";\n if not (type(b1, numeric) and type(b2, numeric) and\n type(b3, numeric) and type(b4, numeric) and\n type(p, num eric)) then RETURN('procname(args)')\n fi;\n if b3 <= -1 then ER ROR(\"Third argument must be >= 1\") fi;\n if b4 <= -1 then ERROR( \"Fourth argument must be >= 1\") fi\n ;\n if p < 0 or 1 < p the n\n ERROR(\"Fifth argument must be between 0 and 1\")\n fi; \n C := GAMMA(b3 + b4 + 2)/(GAMMA(b3 + 1)*GAMMA(b4 + 1));\n CC : = (1/b2)^(b3 + b4 + 1);\n eqn := p =\n int(C*CC*(x - b1)^b3* (b1 + b2 - x)^b4, x = b1 .. y);\n fsolve(eqn, y)\n\nend proc;\n\n\n GammaP := proc(a::algebraic, t::algebraic, p::algebraic)\nlocal eqn, b , x, y;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescription \"G amma percentile\";\n if not (type(a, numeric) and type(t, numeric) \+ and\n type(p, numeric)) then RETURN('procname(args)')\n fi;\n \+ if p < 0 or 1 < p then\n ERROR(\"Third argument must be betwee n 0 and 1\")\n fi;\n if a <= 0 then ERROR(\"First argument must \+ be positive\")\n fi;\n if t <= 0 then ERROR(\"Second argument mu st be positive\")\n fi;\n b := 1/(GAMMA(a)*t^a);\n eqn := p = int(b*x^(a - 1)*exp(- x/t), x = 0 .. y);\n fsolve(eqn, y, y = 0 .. infinity)\n\nend proc;\n\n\nNormalP := proc(mu::algebraic, var::algeb raic, p::algebraic)\nlocal eqn, t, X, x, guess;\noption `August 12, 19 93 -- Zaven A. Karian`;\ndescription \"Normal percentile\";\n if no t (type(mu, numeric) and type(var, numeric) and\n type(p, numeric)) then RETURN('procname(args)')\n fi;\n if p < 0 or 1 < p then\n \+ ERROR(\"Third argument must be between 0 and 1\")\n fi;\n \+ if var <= 0 then\n ERROR(\"Second argument must be positive\") \n fi;\n t := x/sqrt(2);\n guess := 5.063291139*p^.135 - 5.06 3291139*(1 - p)^.135;\n eqn := p = evalf(.5 + .5*erf(t));\n X := fsolve(eqn, x, guess - .5 .. guess + .5);\n sqrt(var)*X + mu\n\nen d proc;\n\n\nTP := proc(df::algebraic, p::algebraic)\nlocal eqn, C1, C 2, C3, x, y;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescriptio n \"Students t percentile\";\n if not (type(df, numeric) and type(p , numeric)) then\n RETURN('procname(args)')\n fi;\n if p \+ < 0 or 1 < p then\n ERROR(\"Second argument must be between 0 a nd 1\")\n fi;\n C1 := GAMMA(1/2*df + 1/2);\n C2 := GAMMA(1/2* df)*sqrt(df*Pi);\n C3 := (1 + x^2/df)^(1/2*df + 1/2);\n eqn := p = int(C1/(C2*C3), x = -infinity .. y);\n fsolve(eqn, y)\n\nend pro c;\n\n\nUniformP := proc(R::range, p::algebraic)\nlocal a, b, x;\nopti on `August 12, 1993 -- Zaven A. Karian`;\ndescription \"uniform percen tile\";\n a := lhs(R);\n b := rhs(R);\n if not (type(a, numer ic) and type(b, numeric) and\n type(p, numeric)) then RETURN('procn ame(args)')\n fi;\n if p < 0 or 1 < p then\n ERROR(\"Seco nd argument must be between 0 and 1\")\n fi;\n if b < a then ERR OR(\"First argument must be a range\") fi;\n fsolve((x - a)/(b - a) = p, x)\n\nend proc;\n\n\nBernoulliCDF := proc(p::algebraic, y::algeb raic)\nlocal i;\noption `August 8, 1993 -- Zaven A. Karian`;\ndescript ion \"Bernoulli CDF\";\n if type(y, numeric) and y < 0 then 0\n \+ elif type(y, numeric) and 1 < y then 1\n else evalf(sum(p^i*(1 - p) ^(1 - i), i = 0 .. y))\n fi\n\nend proc;\n\n\nBetaCDF := proc(b3::a lgebraic, b4::algebraic, y::algebraic)\nlocal C, x, A;\noption `August 8, 1993 -- Zaven A. Karian`;\ndescription \"Beta CDF\";\n if type( y, numeric) and y < 0 then A := 0\n elif type(y, numeric) and 1 < y then A := 1\n else\n C := GAMMA(b3 + b4 + 2)/(GAMMA(b3 + 1) *GAMMA(b4 + 1))\n ;\n A := evalf(int(C*x^b3*(1 - x)^ b4, x = 0 .. y))\n fi;\n A\n\nend proc;\n\n\nBinomialCDF := proc (N::algebraic, p::algebraic, y::algebraic)\nlocal i;\noption `August 8 , 1993 -- Zaven A. Karian`;\ndescription \"Binomial CDF\";\n if typ e(y, numeric) and y < 0 then 0\n elif type(y, numeric) and type(N, \+ numeric) and N < y then\n 1\n else evalf(\n sum(binom ial(N, i)*p^i*(1 - p)^(N - i), i = 0 .. y))\n fi\n\nend proc;\n\n\n ChisquareCDF := proc(r::algebraic, y::algebraic)\nlocal a, x;\noption \+ `August 8, 1993 -- Zaven A. Karian`;\ndescription \"Chi square CDF\"; \n if type(y, numeric) and y < 0 then 0\n else\n a := 1/( GAMMA(1/2*r)*2^(1/2*r));\n evalf(int(a*x^(1/2*r - 1)*exp(- 1/2* x), x = 0 .. y))\n fi\n\nend proc;\n\n\nExponentialCDF := proc(t::a lgebraic, y::algebraic)\nlocal x;\noption `August 8, 1993 -- Zaven A. \+ Karian`;\ndescription \"Exponential CDF\";\n if type(y, numeric) an d y < 0 then 0\n else evalf(int(exp(- x/t)/t, x = 0 .. y))\n fi \n\nend proc;\n\n\nFCDF := proc(df1::algebraic, df2::algebraic, y::alg ebraic)\nlocal num, den, x, A;\noption `August 8, 1993 -- Zaven A. Kar ian`;\ndescription \"F distribution CDF\";\n if type(y, numeric) an d y < 0 then A := 0\n else\n num := x^(1/2*df1 - 1)*GAMMA(1/ 2*df1 + 1/2*df2)*\n (df1/df2)^(1/2*df1);\n den := GA MMA(1/2*df1)*GAMMA(1/2*df2)*\n (1 + df1*x/df2)^(1/2*df1 + 1 /2*df2);\n A := evalf(int(num/den, x = 0 .. y))\n fi;\n A \n\nend proc;\n\n\nGBCDF := proc(b1::algebraic, b2::algebraic,b3::alge braic, b4::algebraic, y::algebraic)\nlocal C, CC, x, A;\noption `Augus t 8, 1993 -- Zaven A. Karian`;\ndescription \"Generalized beta CDF\"; \n if type(b1, numeric) and type(y, numeric) and y < b1 then\n \+ A := 0\n elif type(b2, numeric) and type(y, numeric) and b2 < y \n then A := 1\n else\n C := GAMMA(b3 + b4 + 2)/(GAMMA(b3 + 1)*GAMMA(b4 + 1))\n ;\n CC := (1/b2)^(b3 + b4 + 1 );\n A := evalf(int(C*CC*(x - b1)^b3*(b1 + b2 - x)^b4,\n \+ x = b1 .. y))\n fi;\n A\n\nend proc;\n\n\nGammaCDF := proc( a::algebraic, t::algebraic, y::algebraic)\nlocal b, x, A;\noption `Aug ust 8, 1993 -- Zaven A. Karian`;\ndescription \"Gamma CDF\";\n if t ype(y, numeric) and y < 0 then A := 0\n else\n b := 1/(GAMMA (a)*t^a);\n A := evalf(int(b*x^(a - 1)*exp(- x/t), x = 0 .. y)) \n fi;\n A\n\nend proc;\n\n\nGeometricCDF := proc(p::algebraic, \+ y::algebraic)\nlocal i;\noption `August 8, 1993 -- Zaven A. Karian`;\n description \"Geometric CDF\";\n if type(y, numeric) and y < 1 then 0\n else evalf(sum(p*(1 - p)^(i - 1), i = 1 .. y))\n fi\n\nend \+ proc;\n\n\nNegBinomialCDF := proc(r::algebraic, p::algebraic, y::algeb raic)\nlocal i;\noption `August 8, 1993 -- Zaven A. Karian`;\ndescript ion \"Negative binomial CDF\";\n if type(y, numeric) and type(r, nu meric) and y < r then 0\n else evalf(sum(\n binomial(i - 1, \+ r - 1)*p^r*(1 - p)^(i - r),\n i = r .. y))\n fi\n\nend proc; \n\n\nNormalCDF := proc(mu::algebraic, var::algebraic, y::algebraic)\n local t;\noption `July 7, 1993 -- Zaven A. Karian`;\ndescription \"Nor mal CDF\";\n t := (y - mu)/sqrt(2*var); evalf(.5 + .5*erf(t))\n\nen d proc;\n\n\nPoissonCDF := proc(l::algebraic, y::algebraic)\nlocal i; \noption `Augusr 8, 1993 -- Zaven A. Karian`;\ndescription \"Poisson C DF\";\n if type(y, numeric) and y < 0 then 0\n else evalf(sum(l^ i*exp(-l)/i!, i = 0 .. y))\n fi\n\nend proc;\n\n\nTCDF := proc(df:: algebraic, y::algebraic)\nlocal C1, C2, C3, x;\noption `July 7, 1993 - - Zaven A. Karian`;\ndescription \"Students t CDF\";\n C1 := GAMMA( 1/2*df + 1/2);\n C2 := GAMMA(1/2*df)*sqrt(df*Pi);\n C3 := (1 + x ^2/df)^(1/2*df + 1/2);\n evalf(int(C1/(C2*C3), x = -infinity .. y)) \n\nend proc;\n\n\nUniformCDF := proc(R::range, y::algebraic)\nlocal l , r;\noption `August 8, 1993 -- Zaven A. Karian`;\ndescription \"Unifo rm CDF\";\n l := lhs(R);\n r := rhs(R);\n if type(l, numeric) and type(y, numeric) and y < l then 0\n elif type(r, numeric) and \+ type(y, numeric) and r < y then\n 1\n else evalf((y - l)/(r \+ - l))\n fi\n\nend proc;\n\n\nBernoulliPDF := proc(p::algebraic, x:: algebraic)\noption `August 9, 1993 -- Zaven A. Karian`;\ndescription \+ \"Bernoulli density\";\n if type(p, numeric) and (p <= 0 or 1 <= p) then ERROR(\n \"First argument must be a variable or between 0 an\\\n d 1\")\n fi;\n if type(x, numeric) and\n (x < \+ 0 or 1 < x or not type(x, integer)) then RETURN(0)\n fi;\n p^x*( 1 - p)^(1 - x)\n\nend proc;\n\n\nBetaPDF, proc(b3::algebraic, b4::alge braic, x::algebraic)\nlocal C;\noption `August 11, 1993 -- Zaven A. Ka rian`;\ndescription \"Beta density\";\n if type(b3, numeric) and b3 <= 0 then ERROR(\n \"First argument must be a variable or posi tive\")\n fi;\n if type(b4, numeric) and b4 <= 0 then ERROR(\n \+ \"Second argument must be a variable or positive\")\n fi;\n \+ if type(x, numeric) and (x <= 0 or 1 <= x) then RETURN(0)\n fi;\n C := GAMMA(b3 + b4 + 2)/(GAMMA(b3 + 1)*GAMMA(b4 + 1));\n C*x^b3 *(1 - x)^b4\n\nend proc;\n\n\nBinomialPDF, proc(N::algebraic, p::algeb raic, x::algebraic)\noption `July 7, 1993 -- Zaven A. Karian`;\ndescri ption \"Binomial density\";\n if type(N, numeric) and not type(N, i nteger) then ERROR(\n \"First argument must be a variable or a \+ positive i\\\n nteger\")\n fi;\n if type(N, numeric) and \+ N <= 0 then ERROR(\"First argum\\\n ent must be a variable or a positive integer\")\n fi;\n if type(p, numeric) and (p < 0 or \+ 1 < p) then ERROR(\"S\\\n econd argument must be an expression \+ or between 0 \\\n and 1\")\n fi;\n if type(x, numeric) an d type(N, numeric) and\n (x < 0 or N < x) then RETURN(0)\n fi;\n if type(x, numeric) and (not type(x, integer) or x < 0)\n then \+ RETURN(0)\n fi;\n binomial(N, x)*p^x*(1 - p)^(N - x)\n\nend proc ;\n\n\nBivariateNormalPDF, proc(mux::algebraic,varx::algebraic, x::alg ebraic, \n muy::algebraic,vary::algebraic, y::algebraic, rho::algeb raic)\nlocal A, B, C, Q;\noption `August 11, 1993 -- Zaven A. Karian`; \ndescription \"Bivariate normal density\";\n if type(varx, numeric ) and varx <= 0 then ERROR(\n \"Second argument must be a varia ble or positive\")\n fi;\n if type(vary, numeric) and vary <= 0 \+ then ERROR(\n \"Fifth argument must be a variable or positive\" )\n fi;\n if type(rho, numeric) and (rho < -1 or 1 < rho) then\n ERROR(\"Seventh argument must be a variable or betw\\\n \+ een -1 and +1\")\n fi;\n A := 2*Pi*sqrt(varx*vary*(1 - rho^2)); \n B := (x - mux)^2/varx + (y - muy)^2/vary;\n C := -2*rho*(x - \+ mux)*(y - muy)/sqrt(varx*vary);\n Q := (B + C)/(1 - rho^2);\n ex p(-Q)/A\n\nend proc;\n\n\nCauchyPDF := proc(x::algebraic)\noption `Jul y 16, 1993 -- Zaven A. Karian`;\ndescription \"Cauchy density\";\n \+ 1/(Pi*(1 + x^2))\n\nend proc;\n\n\nChisquarePDF := proc(r::algebraic, \+ x::algebraic)\nlocal a;\noption `August 11, 1993 -- Zaven A. Karian`; \ndescription \"Chi square density\";\n if type(r, numeric) and r < = 0 then ERROR(\n \"First argument must be a variable or positi ve\")\n fi;\n if type(x, numeric) and x <= 0 then RETURN(0) fi; \n a := 1/(GAMMA(1/2*r)*2^(1/2*r));\n a*x^(1/2*r - 1)*exp(- 1/2* x)\n\nend proc;\n\n\nDoubleExponentialPDF := proc(t::algebraic, sig::a lgebraic, x::algebraic)\nlocal A;\noption `August 11, 1993 -- Zaven A. Karian`;\ndescription \"Double exponential density\";\n if type(si g, numeric) and sig <= 0 then ERROR(\n \"Second argument must b e a variable or positive\")\n fi;\n A := - abs(x - t)/sig;\n \+ 1/2*exp(A)/sig\n\nend proc;\n\n\nExponentialPDF := proc(t::algebraic, \+ x::algebraic)\noption `August 11, 1993 -- Zaven A. Karian`;\ndescripti on \"Exponentiall density\";\n if type(t, numeric) and t <= 0 then \+ ERROR(\n \"First argument must be a variable or positive\")\n \+ fi;\n if type(x, numeric) and x < 0 then RETURN(0) fi;\n exp(- x/t)/t\n\nend proc;\n\n\nFPDF := proc(df1::algebraic, df2::algebraic, x::algebraic)\nlocal num, den;\noption `August 11, 1993 -- Zaven A. K arian`;\ndescription \"F density\";\n if type(df1, numeric) and not type(df1, posint) then\n ERROR(\"First argument must be a vari able or a posi\\\n tive integer\")\n fi;\n if type(df2, n umeric) and not type(df2, posint) then\n ERROR(\"Second argumen t must be a variable or a pos\\\n itive integer\")\n fi;\n \+ if type(x, numeric) and x <= 0 then RETURN(0) fi;\n num := x^(1/2 *df1 - 1)*GAMMA(1/2*df1 + 1/2*df2)*\n (df1/df2)^(1/2*df1);\n \+ den := GAMMA(1/2*df1)*GAMMA(1/2*df2)*\n (1 + df1*x/df2)^(1/2*d f1 + 1/2*df2);\n num/den\n\nend proc;\n\n\nGBDPDF := proc(b1::algeb raic, b2::algebraic,b3::algebraic, b4::algebraic, x::algebraic)\nlocal C, CC;\noption `August 11, 1993 -- Zaven A. Karian`;\ndescription \"G eneralized beta density\";\n if type(b3, numeric) and b3 <= -1 then \n ERROR(\"Third argument must be a variable or > -1\")\n fi ;\n if type(b4, numeric) and b4 <= -1 then\n ERROR(\"Fourth \+ argument must be a variable or > -1\")\n fi;\n C := GAMMA(b3 + b 4 + 2)/(GAMMA(b3 + 1)*GAMMA(b4 + 1));\n CC := (1/b2)^(b3 + b4 + 1); \n C*CC*(x - b1)^b3*(b1 + b2 - x)^b4\n\nend proc;\n\n\nGammaPDF := \+ proc(a::algebraic, t::algebraic, x::algebraic)\nlocal b;\noption `Augu st 11, 1993 -- Zaven A. Karian`;\ndescription \"Gamma density\";\n \+ if type(a, numeric) and a <= 0 then ERROR(\n \"First argument m ust be a variable or positive\")\n fi;\n if type(t, numeric) and t <= 0 then ERROR(\n \"Second argument must be a variable or p ositive\")\n fi;\n if type(x, numeric) and x < 0 then RETURN(0) \+ fi;\n b := 1/(GAMMA(a)*t^a);\n b*x^(a - 1)*exp(- x/t)\n\nend pro c;\n\n\nGeometricPDF := proc(p::algebraic, x::algebraic)\noption `Augu st 11, 1993 -- Zaven A. Karian`;\ndescription \"Geometric density\";\n if type(p, numeric) and (p < 0 or 1 < p) then ERROR(\"F\\\n \+ irst argument must be a variable or between 0 and \\\n 1\")\n \+ fi;\n if type(x, numeric) and (not type(x, integer) or x < 1)\n \+ then RETURN(0)\n fi;\n p*(1 - p)^(x - 1)\n\nend proc;\n\n\nHy pergeometricPDF := proc(n1::algebraic, n2::algebraic, r::algebraic, x: :algebraic)\noption `January 27, 1994 -- Zaven A. Karian`;\ndescriptio n \"Hypergeometric density\";\n if type(x, numeric) and type(r, num eric) and r < x then\n ERROR(\"Fourth argument cannot be larger than the t\\\n hird argument\")\n fi;\n if type(x, numer ic) and type(n1, numeric) and n1 < x then\n ERROR(\"Fourth argu ment cannot be larger than the f\\\n irst argument\")\n fi; \n if type(x, numeric) and type(r, numeric) and\n type(n2, numer ic) and n2 < r - x then ERROR(\"2nd argum\\\n ent cannot be sma ller than the 3rd-4th arguments\")\n fi;\n binomial(n1, x)*binom ial(n2, r - x)/binomial(n1 + n2, r)\n\nend proc;\n\n\nLogisticPDF := p roc(x::algebraic)\noption `July 16, 1993 -- Zaven A. Karian`;\ndescrip tion \"Logistic density\";\n exp(-x)/(1 + exp(-x))^2\n\nend proc;\n \n\nNegBinomialPDF := proc(r::algebraic, p::algebraic, x::algebraic)\n option `July 7, 1993 -- Zaven A. Karian`;\ndescription \"Negative bino mial density\";\n if type(p, numeric) and (p < 0 or 1 < p) then ERR OR(\"S\\\n econd argument must be a variable or between 0 and\\ \n 1\")\n fi;\n if type(r, numeric) and not type(r, inte ger) then ERROR(\n \"First argument must be a variable or an in teger\")\n fi;\n if type(x, numeric) and not type(x, integer) th en\n RETURN(0)\n fi;\n if type(x, numeric) and type(r, nu meric) and x < r then\n RETURN(0)\n fi;\n binomial(x - 1, r - 1)*p^r*(1 - p)^(x - r)\n\nend proc;\n\n\nNormalPDF := proc(mu::al gebraic, var::algebraic, x::algebraic)\nlocal a;\noption `August 11, 1 993 -- Zaven A. Karian`;\ndescription \"Normal density\";\n if type (var, numeric) and var <= 0 then ERROR(\n \"Second argument mus t be a variable or positive\")\n fi;\n a := 1/sqrt(2*Pi*var);\n \+ a*exp(- 1/2*(x - mu)^2/var)\n\nend proc;\n\n\nPoissonPDF := proc(l: :algebraic, x::algebraic)\noption `August 11, 1993 -- Zaven A. Karian` ;\ndescription \"Poisson density\";\n if type(l, numeric) and l <= \+ 0 then ERROR(\n \"First argument must be a variable or positive \")\n fi;\n if type(x, numeric) and (not type(x, integer) or x < 0)\n then RETURN(0)\n fi;\n l^x*exp(-l)/x!\n\nend proc;\n\n \nTPDF := proc(df::algebraic, x::algebraic)\nlocal C1, C2, C3;\noption `August 11, 1993 -- Zaven A. Karian`;\ndescription \"Students t densi ty\";\n if type(df, numeric) and not type(df, posint) then ERROR(\n \"First argument must be a variable or a positive i\\\n \+ nteger\")\n fi;\n C1 := GAMMA(1/2*df + 1/2);\n C2 := GAMMA(1 /2*df)*sqrt(df*Pi);\n C3 := (1 + x^2/df)^(1/2*df + 1/2);\n C1/(C 2*C3)\n\nend proc;\n\n\nUniformPDF := proc(R::range, x::algebraic)\nlo cal a, b;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescription \+ \"Uniform density\";\n a := lhs(R);\n b := rhs(R);\n if type( x, numeric) and type(a, numeric) and x < a then\n RETURN(0)\n \+ fi;\n if type(x, numeric) and type(b, numeric) and b < x then\n \+ RETURN(0)\n fi;\n 1/(b - a)\n\nend proc;\n\n\nWeibullPDF : = proc(a::algebraic, b::algebraic, x::algebraic)\nlocal A;\noption `Au gust 11, 1993 -- Zaven A. Karian`;\ndescription \"Weibull density\";\n if type(a, numeric) and a <= 0 then ERROR(\n \"First argume nt must be a variable or positive\")\n fi;\n if type(b, numeric) and b <= 0 then ERROR(\n \"Second argument must be a variable \+ or positive\")\n fi;\n if type(x, numeric) and x < 0 then RETORN (0) fi;\n A := a*(x/b)^(a - 1)/b;\n A*exp(-(x/b)^a)\n\nend proc; \n\n\nExponentialSumS := proc(t::numeric, n::posint, m::posint)\nlocal X, Y, i, j;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescriptio n \"Generates sums from exponential\";\n if t <= 0 then ERROR(\"Fir st argument must be positive\")\n fi;\n Y := array(1 .. m);\n \+ for i to m do\n X := ExponentialS(t, n);\n Y[i] := sum( X[j], j = 1 .. n)\n od;\n [seq(Y[i], i = 1 .. m)]\n\nend proc;\n \n\nNormTransVarS := proc(a::numeric, b::numeric, n::posint, m::posint )\nlocal X, Y, i;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescr iption \"Generates transformed variances from normal\";\n if b <= 0 then ERROR(\"Second argument must be positive\")\n fi;\n Y := a rray(1 .. m);\n for i to m do\n X := NormalS(a, b, n);\n \+ Y[i] := ((n - 1)*Variance(X))/b\n od;\n [seq(Y[i], i = 1 .. \+ m)]\n\nend proc;\n\n\nNormalMeanS := proc(a::numeric, b::numeric, n::p osint, m::posint)\noption `May 21, 1994 -- Zaven A. Karian`;\ndescript ion \"Generates means from normal\";\n if b <= 0 then ERROR(\"Secon d argument must be positive\")\n fi;\n NormalS(a, b/n, m)\n\nend proc;\n\n\nNormalMedianS := proc(a::numeric, b::numeric, n::posint, m ::posint)\nlocal X, Y, i;\noption `August 12, 1993 -- Zaven A. Karian` ;\ndescription \"Generates medians from normal\";\n if b <= 0 then \+ ERROR(\"Second argument must be positive\")\n fi;\n Y := array(1 .. m);\n for i to m do\n X := NormalS(a, b, n); Y[i] := Per centile(X, .5)\n od;\n [seq(Y[i], i = 1 .. m)]\n\nend proc;\n\n \nNormalSumS := proc(a::numeric, b::numeric, n::posint, m::posint)\nop tion `May 21, 1994 -- Zaven A. Karian`;\ndescription \"Generates sums \+ from normal\";\n if b <= 0 then ERROR(\"Second argument must be pos itive\")\n fi;\n NormalS(a, n*b, m)\n\nend proc;\n\n\nNormalVari anceS := proc(a::numeric, b::numeric, n::posint, m::posint)\nlocal X, \+ Y, i;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescription \"Gen erates variances from normal\";\n if b <= 0 then ERROR(\"Second arg ument must be positive\")\n fi;\n Y := array(1 .. m);\n for i to m do X := NormalS(a, b, n); Y[i] := Variance(X)\n od;\n [seq (Y[i], i = 1 .. m)]\n\nend proc;\n\n\nUniformMeanS := proc(R::range, n ::posint, m::posint)\nlocal X, Y, i;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescription \"Generates means from continuous uniform\"; \n if not (type(lhs(R), numeric) and type(rhs(R), numeric))\n o r rhs(R) < lhs(R) then\n ERROR(\"First argument must be a numer ic range\")\n fi;\n Y := array(1 .. m);\n for i to m do X := \+ UniformS(R, n); Y[i] := Mean(X) od;\n [seq(Y[i], i = 1 .. m)]\n\nen d proc;\n\n\nUniformMedianS := proc(R::range, n::posint, m::posint)\nl ocal X, Y, i;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescripti on \"Generates medians from continuous uniform\";\n if not (type(lh s(R), numeric) and type(rhs(R), numeric))\n or rhs(R) < lhs(R) the n\n ERROR(\"First argument must be a numeric range\")\n fi; \n Y := array(1 .. m);\n for i to m do\n X := UniformS(R, n); Y[i] := Percentile(X, .5)\n od;\n [seq(Y[i], i = 1 .. m)]\n \nend proc;\n\n\nUniformSumS := proc(R::range, n::posint, m::posint)\n local X, Y, i;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescript ion \"Generates sums from continuous uniform\";\n if not (type(lhs( R), numeric) and type(rhs(R), numeric))\n or rhs(R) < lhs(R) then \n ERROR(\"First argument must be a numeric range\")\n fi;\n Y := array(1 .. m);\n for i to m do\n X := UniformS(R, n ); Y[i] := convert(X, `+`)\n od;\n [seq(Y[i], i = 1 .. m)]\n\nen d proc;\n\n \nBernoulliS := proc(p::numeric, n::posint)\nlocal P;\nopt ion `August 12, 1993 -- Zaven A. Karian`;\ndescription \"n random nos. from Bernoulli dist.\";\n if p < 0 or 1 < p then\n ERROR(\" First argument must be between 0 and 1\")\n fi;\n P := [0, 1 - p , 1, p];\n DiscreteS(P, n)\n\nend proc;\n\n\nBetaS := proc(n1::nume ric, n2::numeric, n::posint)\nlocal i, y, yy, A, X;\noption `August 12 , 1993 -- Zaven A. Karian`;\ndescription \"Generates sample from beta \";\n if n1 <= 0 then ERROR(\"First argument must be positive\")\n \+ fi;\n if n2 <= 0 then ERROR(\"Second argument must be positive\" )\n fi;\n X := array(1 .. n);\n A := FS(2*n2 + 2, 2*n1 + 2, n );\n for i to n do\n yy := ((n2 + 1)*A[i])/(n1 + 1);\n \+ y := 1/(yy + 1);\n X[i] := y\n od;\n [seq(X[i], i = 1 . . n)]\n\nend proc;\n\n\nBinomialS := proc(N::posint, p::numeric, n::po sint)\nlocal A, i, L;\noption `August 12, 1993 -- Zaven A. Karian`;\nd escription \"n random nos. from B(n,p)\";\n if p < 0 or 1 < p then \n ERROR(\"Second argument must be between 0 and 1\")\n fi; \n A := array(1 .. 2*N + 2);\n for i from 0 to N do\n A[2 *i + 1] := i;\n A[2*i + 2] := binomial(N, i)*p^i*(1 - p)^(N - i )\n od;\n L := [seq(A[i], i = 1 .. 2*N + 2)];\n DiscreteS(L, \+ n)\n\nend proc;\n\n\nBivariateNormalS := proc(mux::numeric, varx::nume ric,muy::numeric, \n vary::numeric, rho::numeric, n::posint)\nlocal i, u, v, a, b, cmua, cmuy, cvar, A;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescription \"Generates sample from bivariate normal\"; \n if varx <= 0 then\n ERROR(\"Second argument must be posit ive\")\n fi;\n if vary <= 0 then\n ERROR(\"Fourth argumen t must be positive\")\n fi;\n if rho < -1 or 1 < rho then\n \+ ERROR(\"Fifth argument must be between -1 and +1\")\n fi;\n c var := vary*(1 - rho^2);\n cmua := rho*sqrt(vary/varx);\n A := a rray(1 .. n, 1 .. 2);\n for i to n do\n u := rng();\n \+ v := rng();\n a := sqrt(-2*log(u))*sin(2*Pi*v);\n a := \+ evalf(sqrt(varx)*a + mux);\n A[i, 1] := a;\n cmuy := muy + cmua*(a - mux);\n b := sqrt(-2*log(u))*cos(2*Pi*v);\n \+ b := evalf(sqrt(cvar)*b + cmuy);\n A[i, 2] := b\n od;\n \+ [seq([A[i, 1], A[i, 2]], i = 1 .. n)]\n\nend proc;\n\n\nCauchyS := pro c(n::posint)\nlocal i;\noption `July 6, 1993 -- Zaven A. Karian`;\ndes cription \"Generates sample from Cauchy\";\n [seq(evalf(tan(Pi*rng( ) - 1/2*Pi)), i = 1 .. n)]\n\nend proc;\n\n\nChisquareS := proc(r::pos int, n::posint)\nlocal i, j, X, A, t;\noption `May 27, 1994 -- Zaven A . Karian`;\ndescription \"Generates sample from standard normal\";\n \+ A := array(1 .. n);\n X := N01S(r*n);\n j := 1;\n for i to \+ n do\n A[i] :=\n convert([seq(X[t]^2, t = j .. j + r - 1)], `+`);\n j := j + r\n od;\n [seq(A[i], i = 1 .. n) ]\n\nend proc;\n\n\nContinuousS := proc(expr::algebraic, R::range, n:: posint)\nlocal i, var, x, F, CDF, t, ll, ul, pts, Y, j, Rand, a, b, c, d;\noption `January 24, 1994 -- Zaven A. Karian`;\ndescription \"Gen \+ sample (size n) with pdf given by expr on \\\nthe interval R\";\n v ar := op(select(type, indets(expr), name));\n ll := lhs(args[2]);\n ul := rhs(args[2]);\n F := int(expr, var = ll .. t);\n CDF : = student[makeproc](F, t);\n pts := array(1 .. 72);\n Y := array (1 .. n);\n j := 3;\n for i from 2 to 35 do\n x := ll + 1 /35*(i - 1)*(ul - ll);\n pts[j] := CDF(x);\n pts[j + 1] \+ := x;\n j := j + 2\n od;\n pts[1] := 0;\n pts[2] := ll ;\n pts[71] := 1;\n pts[72] := ul;\n Rand := RNG(n);\n for i to n do\n j := 1;\n while pts[2*j - 1] < Rand[i] do j := j + 1 od;\n a := pts[2*j - 3];\n b := pts[2*j - 2]; \n c := pts[2*j - 1];\n d := pts[2*j];\n Y[i] := \+ b + ((Rand[i] - a)*(d - b))/(c - a)\n od;\n [seq(Y[i], i = 1 .. \+ n)]\n\nend proc;\n\n\nDie := proc(m::posint, n::posint)\nlocal die, i; \noption `July 5, 1993 -- Zaven A. Karian`;\ndescription \"n rolls of \+ an m-sided die\";\n die := rand(1 .. m); [seq(die(), i = 1 .. n)]\n \nend proc;\n\n\nDiscUniformS := proc(R::range, n::posint)\nlocal i, a , b;\noption `July 6, 1993 -- Zaven A. Karian`;\ndescription \"Generat es sample from discrete uniform\";\n a := lhs(R);\n b := rhs(R); \n [seq(a + trunc((b - a + 1)*rng()), i = 1 .. n)]\n\nend proc;\n\n \nDiscreteS := proc()\nlocal N, A, B, ll, ul, ptr, var, F, i, r, j, L, expr, n;\noption `July 5, 1993 -- Zaven A. Karian`;\ndescription \"n \+ rand. obs. from dist. described by L\";\n if nargs = 2 and not type (args[1], list) then ERROR(\"W\\\n hen two arguments are used, \+ the first must be a li\\\n st\")\n fi;\n if nargs = 2 and not type(args[2], posint) then ERROR(\n \"When two arguments a re used, the second must be a\\\n positive integer\")\n fi; \n if nargs = 3 and not type(args[1], algebraic) then ERROR(\n \+ \"When three arguments are used, the first must be \\\n an e xpression\")\n fi;\n if nargs = 3 and not type(args[2], range) t hen ERROR(\"\\\n When three arguments are used, the second must be \\\n a range\")\n fi;\n if nargs = 3 and not type(arg s[3], posint) then ERROR(\n \"When three arguments are used, th e third must be \\\n a positive integer\")\n fi;\n if nar gs = 3 then\n expr := args[1];\n n := args[3];\n \+ var := op(select(type, indets(expr), name));\n ll := lhs(args[2 ]);\n ul := rhs(args[2]);\n B := array(1 .. 2*ul - 2*ll \+ + 2);\n ptr := 1;\n for i from ll to ul do\n \+ B[ptr] := i;\n B[ptr + 1] := subs(var = i, expr);\n \+ ptr := ptr + 2\n od;\n L := convert(B, list)\n el se L := args[1]; n := args[2]\n fi;\n N := 1/2*nops(L);\n A : = array(1 .. n);\n F := array(1 .. N);\n F[1] := L[2];\n for \+ i from 2 to N do F[i] := F[i - 1] + L[2*i] od;\n for i to n do\n \+ r := rng();\n if r < F[1] then A[i] := L[1] fi;\n i f F[N] < r then A[i] := L[2*N - 1] fi;\n for j from 2 to N do\n if r < F[j] and F[j - 1] <= r then\n A[i] : = L[2*j - 1]\n fi\n od\n od;\n [seq(A[i], i = \+ 1 .. n)]\n\nend proc;\n\n\nExponentialS := proc(t::numeric, n::posint) \nlocal i;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescription \+ \"Generates sample from exponential\";\n if t <= 0 then ERROR(\"Fir st argument must be positive\")\n fi;\n [seq(evalf(- t*log(1 - r ng())), i = 1 .. n)]\n\nend proc;\n\n\nFS := proc(n1::posint, n2::posi nt, n::posint)\nlocal i, A, B;\noption `August 12, 1993 -- Zaven A. Ka rian`;\ndescription \"Generates sample from F\";\n A := ChisquareS( n1, n);\n B := ChisquareS(n2, n);\n [seq(A[i]*n2/(n1*B[i]), i = \+ 1 .. n)]\n\nend proc;\n\n\nGammaS := proc(a::numeric, t::numeric, n::p osint)\nlocal WW, V, X, ay, Zee, d, q, theta, k, j, Y, U1, U2, b, P;\n if a <= 0 then ERROR(\"First argument must be positive\")\n fi; \n if t <= 0 then ERROR(\"Second argument must be positive\")\n \+ fi;\n X := array(1 .. n);\n k := 1;\n while k <= n do\n \+ if a < 1 then\n b := evalf((exp(1) + a)/exp(1));\n \+ U1 := RNG(1)[1];\n P := b*U1;\n if 1 < P th en\n Y := evalf(-ln((b - P)/a));\n U2 := RNG(1)[1];\n if U2 <= evalf(Y^(a - 1)) then\n \+ X[k] := t*Y; k := k + 1\n fi\n el se\n Y := P^(1/a);\n U2 := RNG(1)[1];\n \+ if U2 <= evalf(exp(-Y)) then\n X[k] \+ := t*Y; k := k + 1\n fi\n fi\n else\n ay := evalf(1/sqrt(2*a - 1));\n b := evalf(a - \+ ln(4));\n q := evalf(a + 1/ay);\n theta := 4.5; \n d := 1 + ln(theta);\n U1 := RNG(1)[1];\n \+ U2 := RNG(1)[1];\n V := evalf(ay*ln(U1/(1 - U1)));\n Y := evalf(a*exp(V));\n Zee := U1^2*U2;\n \+ WW := b + q*V - Y;\n if 0 <= evalf(WW + d - theta*Zee ) then\n X[k] := t*Y; k := k + 1\n else\n \+ if evalf(ln(Zee)) <= WW then\n X[k] := t*Y; k := k + 1\n fi\n fi\n fi\n \+ od;\n [seq(X[j], j = 1 .. n)]\n\nend proc;\n\n\nGeometricS := proc( p::numeric, n::posint)\nlocal A, SumProb, i, Prob, UpperLimit, L;\nopt ion `August 12, 1993 -- Zaven A. Karian`;\ndescription \"n random nos. from geometric dist.\";\n if p < 0 or 1 < p then\n ERROR(\" First argument must be between 0 and 1\")\n fi;\n A := array(1 . . 500);\n SumProb := 0;\n i := 1;\n if p = 1 then A[1] := 1; \+ A[2] := 1; UpperLimit := 2\n else\n while SumProb < .999 do \n Prob := p*(1 - p)^(i - 1);\n SumProb := SumPr ob + Prob;\n A[2*i - 1] := i;\n A[2*i] := Prob; \n i := i + 1\n od;\n A[2*i - 1] := i;\n \+ A[2*i] := 1 - SumProb;\n UpperLimit := 2*i\n fi;\n L : = [seq(A[i], i = 1 .. UpperLimit)];\n DiscreteS(L, n)\n\nend proc; \n\n\nHypergeometricS :=proc(n1::posint, n2::posint, r::posint, n::pos int)\nlocal lb, ub, A, ptr, x, L, i;\noption `January 27, 1994 -- Zave n A. Karian`;\ndescription \"Hypergeometric density\";\n lb := max( 0, r - n2);\n ub := min(r, n1);\n A := array(1 .. 2*ub - 2*lb + \+ 2);\n ptr := 1;\n for x from lb to ub do\n A[ptr] := x;\n A[ptr + 1] := binomial(n1, x)*binomial(n2, r - x)/\n \+ binomial(n1 + n2, r);\n ptr := ptr + 2\n od;\n L := [se q(A[i], i = 1 .. 2*ub - 2*lb + 2)];\n DiscreteS(L, n)\n\nend proc; \n\n\nLogisticS := proc(n::posint)\nlocal a, r, i, X;\noption `July 6, 1993 -- Zaven A. Karian`;\ndescription \"Generates sample from logist ic\";\n X := array(1 .. n);\n for i to n do\n r := rng(); a := evalf(log(r/(1 - r))); X[i] := a\n od;\n [seq(X[i], i = 1 \+ .. n)]\n\nend proc;\n\n\nN01S := proc(n::posint)\nlocal j, i, u, v, X, P2;\noption `July 6, 1993 -- Zaven A. Karian`;\ndescription \"Generat es sample from standard normal\";\n X := array(1 .. n + 1);\n j \+ := 1;\n P2 := evalf(2*Pi);\n for i to 1/2*n + 1/2 do\n u \+ := rng();\n v := rng();\n X[j] := evalf(sqrt(-2*log(u))* sin(P2*v));\n j := j + 1;\n X[j] := evalf(sqrt(-2*log(u) )*cos(P2*v));\n j := j + 1\n od;\n [seq(X[i], i = 1 .. n) ]\n\nend proc;\n\n\nNegBinomialS := proc(r::posint, p::numeric, n::pos int)\nlocal A, SumProb, i, Prob, UpperLimit, L;\noption `August 12, 19 93 -- Zaven A. Karian`;\ndescription \"n random nos. from negative bin omial dist.\";\n if p < 0 or 1 < p then\n ERROR(\"Second arg ument must be between 0 and 1\")\n fi;\n A := array(1 .. 100);\n SumProb := 0;\n i := r;\n while SumProb < .999 do\n P rob :=\n evalf(binomial(i - 1, r - 1)*p^r*(1 - p)^(i - r)) \n ;\n SumProb := SumProb + Prob;\n A[2*i - 2 *r + 1] := i;\n A[2*i - 2*r + 2] := Prob;\n i := i + 1\n od;\n A[2*i - 2*r + 1] := i;\n A[2*i - 2*r + 2] := 1 - SumPr ob;\n UpperLimit := 2*i - 2*r + 2;\n L := [seq(A[i], i = 1 .. Up perLimit)];\n DiscreteS(L, n)\n\nend proc;\n\n\nNormalS := proc(mu: :numeric, var::numeric, n::posint)\nlocal j, i, P2, u, v, a, b, X, sqr t_var;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescription \"Ge nerates sample from normal\";\n if var <= 0 then\n ERROR(\"S econd argument must be positive\")\n fi;\n X := array(1 .. n + 1 );\n j := 1;\n P2 := evalf(2*Pi);\n sqrt_var := sqrt(var);\n \+ for i to 1/2*n + 1/2 do\n u := rng();\n v := rng();\n a := sqrt(-2*log(u))*sin(P2*v);\n X[j] := evalf(sqrt_va r*a + mu);\n j := j + 1;\n b := sqrt(-2*log(u))*cos(P2*v );\n X[j] := evalf(sqrt_var*b + mu);\n j := j + 1\n o d;\n [seq(X[i], i = 1 .. n)]\n\nend proc;\n\n\nPoissonS := proc(l:: numeric, n::posint)\nlocal A, SumProb, i, Prob, UpperLimit, L;\noption `July 5, 1993 -- Zaven A. Karian`;\ndescription \"n random nos. from \+ Poisson dist.\";\n if l <= 0 then ERROR(\"First argument must be po sitive\")\n fi;\n A := array(1 .. 100);\n SumProb := 0;\n \+ i := 0;\n while SumProb < .999 do\n Prob := evalf(l^i*exp(-l )/i!);\n SumProb := SumProb + Prob;\n A[2*i + 1] := i;\n A[2*i + 2] := Prob;\n i := i + 1\n od;\n A[2*i + \+ 1] := i;\n A[2*i + 2] := 1 - SumProb;\n UpperLimit := 2*i + 2;\n L := [seq(A[i], i = 1 .. UpperLimit)];\n DiscreteS(L, n)\n\nend proc;\n\n\nRNG := proc(n::posint)\nlocal i;\noption `July 5, 1993 -- \+ Zaven A. Karian`;\ndescription \"n random numbers from [0,1]\";\n [ seq(Float(rand(), -12), i = 1 .. n)]\n\nend proc;\n\n\nTS := proc(nu:: posint, n::posint)\nlocal i, A, B;\noption `August 12, 1993 -- Zaven A . Karian`;\ndescription \"Generates sample from Student's t\";\n A \+ := N01S(n);\n B := ChisquareS(nu, n);\n [seq(A[i]/sqrt(B[i]/nu), i = 1 .. n)]\n\nend proc;\n\n\nUniformS := proc(R::range, n::posint) \nlocal a, b, i;\noption `August 12, 1993 -- Zaven A. Karian`;\ndescri ption \"Generates sample from continuous uniform\";\n a := lhs(R); \n b := rhs(R);\n if not (type(a, numeric) and type(b, numeric)) or b < a\n then ERROR(\n \"First argument must be a proper \+ numeric range\")\n fi;\n [seq(evalf((b - a)*rng() + a), i = 1 .. n)]\n\nend proc;\n\n\nMakeRandom := proc()\nglobal _seed;\noption `Ju ly 5, 1993 -- Zaven A. Karian`;\ndescription \"n random numbers from [ 0,1)\";\n _seed := round(123456789*time())\n\nend proc;\n\n\nrng := () -> Float(rand(), -12):\n\n\nClassFreq := proc(L::list, R::range, n ::posint)\nlocal N, i, A, SL, bdry, delta, ptr;\noption `January 26, 1 994 -- Zaven A. Karian`;\ndescription\n\"Class Frequencies of L on the range R in n intervals\";\n if nargs = 2 then N := 8 else N := n f i;\n A := array(1 .. N);\n SL := sort(L);\n delta := (rhs(R) \+ - lhs(R))/N;\n ptr := 1;\n for i to N do\n A[i] := 0;\n \+ bdry := i*delta + lhs(R);\n while ptr <= nops(SL) and SL[ ptr] <= bdry do\n A[i] := A[i] + 1; ptr := ptr + 1\n \+ od\n od;\n [seq(A[i], i = 1 .. N)]\n\nend proc;\n\n\nFreq := pr oc(L::list, R::range)\nlocal count, N, i, Min, Max;\noption `August 12 , 1993 -- Zaven A. Karian`;\ndescription \"Frequencies of L on the ran ge R\";\n Min := lhs(R);\n Max := rhs(R);\n if not (type(Min, numeric) and type(Max, numeric)) or\n Max < Min then ERROR(\n \+ \"Second argument must be a propoer numerical range\")\n fi;\n \+ count := array(Min .. Max);\n N := nops(L);\n for i from Min t o Max do count[i] := 0 od;\n for i to N do\n if Min <= L[i] \+ and L[i] <= Max then\n count[L[i]] := count[L[i]] + 1\n \+ fi\n od;\n [seq(count[i], i = Min .. Max)]\n\nend proc;\n\n \nKurtosis := proc(L::list)\nlocal sf, M, N, V, i, X;\noption `August \+ 19, 1993 -- Zaven A. Karian`;\ndescription \"Kurtosis or 3rd standardi zed moment\";\n if type(args[1], listlist) then\n X := [seq( L[i][1] $ L[i][2], i = 1 .. nops(L))]\n else X := L\n fi;\n s f := 0;\n M := Mean(X);\n N := nops(X);\n V := ((N - 1)*Varia nce(X))/N;\n for i to N do sf := sf + (X[i] - M)^4 od;\n sf/(N*V ^2)\n\nend proc;\n\n\nLocate := proc(X::list, item::algebraic)\nlocal \+ n, A, ptr, i, j;\noption `January 24, 1994 -- Zaven A. Karian`;\ndescr iption \"Find the indeces in the list X of the item\";\n n := nops( X);\n A := array(1 .. n);\n ptr := 1;\n for i to n do\n \+ if X[i] = item then A[ptr] := i; ptr := ptr + 1 fi\n od;\n [se q(A[j], j = 1 .. ptr - 1)]\n\nend proc;\n\n\nMax := proc(L::list)\nopt ion `July 5, 1993 -- Zaven A. Karian`;\ndescription \"Maximum value of the list L\";\n max(op(L))\n\nend proc;\n\n\nMean := proc(L::list) \nlocal S, N, i;\noption `July 5, 1993 -- Zaven A. Karian`;\ndescripti on \"Arithmetic mean of the list L\";\n if type(args[1], listlist) \+ then\n S := 0;\n N := 0;\n for i to nops(L) do\n \+ S := S + L[i][1]*L[i][2]; N := N + L[i][2]\n od\n \+ else S := convert(L, `+`); N := nops(L)\n fi;\n S/N\n\nend proc ;\n\nMedian := proc(L::list)\noption `May 26, 1994 -- Zaven A. Karian` ;\ndescription \"The 1n of the list L\";\n Percentile(L, .5)\n\nend proc;\n\n\nMin := proc(L::list)\noption `July 5, 1993 -- Zaven A. Kar ian`;\ndescription \"Minimum value of the list L\";\n min(op(L))\n \nend proc;\n\n\nPercentile := proc(L::list(numeric), p::numeric)\nloc al n, LL, F, r, ab;\noption `August 12, 1993 -- Zaven A. Karian`;\ndes cription \"The 100*p-th percentile of the list L\";\n if p < 0 or 1 < p then\n ERROR(\"Second argument must be between 0 and 1\") \n fi;\n n := nops(L);\n if p < 1/(n + 1) or n/(n + 1) < p th en\n ERROR(\"Percentile does not exist\")\n fi;\n LL := s ort(L);\n F := convert((n + 1)*p, fraction);\n r := trunc(F);\n \+ ab := F - r;\n LL[r] + ab*(LL[r + 1] - LL[r])\n\nend proc;\n\n\n Range := proc(L::list)\noption `July 5, 1993 -- Zaven A. Karian`;\ndes cription \"Range of the list L\";\n Max(L) - Min(L)\n\nend proc;\n \n\nRunningSum := proc(L::list)\nlocal N, A, i;\noption `August 23, 19 93 -- Zaven A. Karian`;\ndescription \"Running sum of a list\";\n N := nops(L);\n A := array(1 .. N);\n A[1] := L[1];\n for i fr om 2 to N do A[i] := A[i - 1] + L[i] od;\n [seq(A[i], i = 1 .. N)] \n\nend proc;\n\n\nSkewness := proc(L::list)\nlocal sc, M, N, V, i, X; \noption `August 19, 1993 -- Zaven A. Karian`;\ndescription \"Skewness or 3rd standardized moment\";\n if type(args[1], listlist) then\n \+ X := [seq(L[i][1] $ L[i][2], i = 1 .. nops(L))]\n else X := \+ L\n fi;\n sc := 0;\n M := Mean(X);\n N := nops(X);\n V \+ := ((N - 1)*Variance(X))/N;\n for i to N do sc := sc + (X[i] - M)^3 od;\n sc/(N*V^(3/2))\n\nend proc;\n\n\nStDev := proc(L::list)\nopt ion `July 5, 1993 -- Zaven A. Karian`;\ndescription \"Standard deviati on of the list L\";\n sqrt(Variance(L))\n\nend proc;\n\n\nStemLeaf \+ := proc(L::list(numeric))\nlocal n, stems, digits, Ls, range, logscale , x, fracwidth,\nminwidth, maxwidth, M, Label, labels, D, i, width, st emnum,\nstart, format, leaf, stem, dispstem, tenstem;\noption `Joshua \+ D. Levy ~ June 10, 1994`;\ndescription \"Ordered stem-and-leaf displ ay\";\n M := [1, 2, 5, 10];\n Label :=\n '[[\" \"], [\"* \", t, f, s, \"+\"], [\"*\", \"+\"], [\" \"]]';\n n := nops(L);\n \+ Ls := sort(evalf(L));\n range := Ls[n] - Ls[1];\n if 1 < nargs and type(args[2], posint) then\n digits := args[2]\n else d igits := min(length(op(1, range)), 2)\n fi;\n if 2 < nargs and t ype(args[3], posint) then\n stems := args[3]\n else stems := round(evalf(sqrt(n)))\n fi;\n logscale := digits - ilog10(range /stems) - 1;\n Ls := [seq(round(10^logscale*x), x = Ls)];\n frac width := (Ls[n] - Ls[1])/stems;\n minwidth := 10^ilog10(fracwidth); \n maxwidth := 10*minwidth;\n D := seq(evalf(abs(fracwidth - M[i ]*minwidth)),\n i = 1 .. 4);\n member(min(D), [D], 'i');\n \+ width := M[i]*minwidth;\n labels := nops(Label[i]);\n stemnum \+ := iquo(Ls[1], width);\n start := width*stemnum;\n leaf := 1;\n \+ format := cat(\" %0\", length(minwidth), \".0f\");\n printf(\"\\ n Stem Leaf\");\n for stem from start by width while stem <= Ls [n] do\n dispstem := iquo(stem, maxwidth, 'r');\n tenste m := maxwidth*dispstem;\n printf(\"\\n%6.0f%c |\", dispstem,\n \+ Label[i][1 + irem(stemnum, labels)]);\n for leaf fro m leaf while\n leaf <= n and Ls[leaf] < stem + width do\n \+ printf(format, Ls[leaf] - tenstem)\n od;\n stemnum := stemnum + 1\n od;\n if logscale <> 0 then printf(\"\\n \+ \\\n %.0f\\n(Multiply numbers by 10 .)\", -log scale)\n fi;\n lprint()\n\nend proc;\n\n\nVariance := proc(L::li st)\nlocal ssq, s, N, i, X;\noption `July 5, 1993 -- Zaven A. Karian`; \ndescription \"Variance of the list L\";\n if type(args[1], listli st) then\n X := [seq(L[i][1] $ L[i][2], i = 1 .. nops(L))]\n \+ else X := L\n fi;\n N := nops(X);\n s := convert(X, `+`);\n \+ ssq := convert([seq(X[i]^2, i = 1 .. N)], `+`);\n (ssq - s^2/N)/ (N - 1)\n\nend proc;\n\nend module:\n" }}}{EXCHG {PARA 4 "" 0 "" {TEXT -1 25 "Installation Instructions" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 12 "To save the " }{TEXT 262 4 "stat" }{TEXT -1 65 " package \+ into a library, remove the sharp character (#) from the " }{TEXT 259 7 "libname" }{TEXT -1 22 " assignment below and " }{TEXT 258 0 "" } {TEXT -1 124 "replace the example directory name \"C:/mylib/statpackag e\" with the name of the directory where you want the package saved. \+ " }{TEXT 264 45 "Make sure that this directory exists first! " } {TEXT -1 38 "Then remove the # characters from the " }{TEXT 260 5 "mar ch" }{TEXT -1 5 " and " }{TEXT 261 8 "savelib " }{TEXT -1 69 "commands . Then execute this worksheet by choosing from the menu bar " }} {PARA 0 "" 0 "" {TEXT 263 24 "Edit->Execute->Worksheet" }{TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "#libname:=\"C:/mylib/statpackage\",libname; " }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "#march('create',libname[1],1 00);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "#savelib('stat');" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 79 "Once you've saved the package, you can use it in any Maple worksheet \+ by typing:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 42 "#libname:=\"C:/mylib/statpackage\",libname; " }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "#with(stat);" }}}}{MARK "9 1 0" 79 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }