home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / STATISTI.PAK / DISCRETE.M < prev    next >
Encoding:
Text File  |  1992-07-29  |  26.3 KB  |  794 lines
  1.  
  2. (*:Version: Mathematica 2.0 *)
  3.  
  4. (*:Name: Statistics`DiscreteDistributions *)
  5.  
  6. (*:Title: Discrete Statistical Distributions *)
  7.  
  8. (*:Author:
  9.   David Withoff (Wolfram Research), April, 1990.
  10.   Revised March, 1991.
  11. *)
  12.  
  13. (*:Legal:
  14.   Copyright (c) 1990, 1991, Wolfram Research, Inc.
  15. *)
  16.  
  17. (*:Reference: Usage messages only. *)
  18.  
  19. (*:Summary:
  20.   Properties and functionals of discrete probability
  21.   distributions used in statistical computations.
  22. *)
  23.  
  24. (*:Keywords: discrete distribution *)
  25.  
  26. (*:Requirements: No special system requirements. *)
  27.  
  28. (*:Warning:
  29.   This package extends the definition of several descriptive statistics
  30.   functions and the definition of Random.  If the original usage messages
  31.   are reloaded, this change will not be reflected in the usage message,
  32.   although the extended functionality will remain.
  33. *)
  34.  
  35. (*:Sources:
  36.   Norman L. Johnson and Samuel Kotz, Discrete Distributions, 1969.
  37.   Stephen Kokoska and Christopher Nevison, Statistical Tables and
  38.     Formulae, 1989.
  39. *)
  40.  
  41. BeginPackage["Statistics`DiscreteDistributions`",
  42.              "Statistics`DescriptiveStatistics`",
  43.              "Statistics`Common`DistributionsCommon`"]
  44.  
  45.      (* Distributions *)
  46.  
  47. BernoulliDistribution::usage =
  48. "BernoulliDistribution[p] represents the Bernoulli distribution with mean p."
  49.  
  50. BinomialDistribution::usage =
  51. "BinomialDistribution[n, p] represents the Binomial distribution
  52. for n trials with probability p."
  53.  
  54. DiscreteUniformDistribution::usage =
  55. "DiscreteUniformDistribution[n] represents the discrete Uniform
  56. distribution with n possible outcomes."
  57.  
  58. GeometricDistribution::usage =
  59. "GeometricDistribution[p] represents the Geometric distribution
  60. for probability p."
  61.  
  62. HypergeometricDistribution::usage =
  63. "HypergeometricDistribution[n, nsucc, ntot] represents the
  64. Hypergeometric distribution for samples of size n drawn without
  65. replacement from a population with nsucc successes and total size
  66. ntot."
  67.  
  68. LogSeriesDistribution::usage =
  69. "LogSeriesDistribution[theta] represents the logarithmic series
  70. distribution with parameter theta."
  71.  
  72. NegativeBinomialDistribution::usage =
  73. "NegativeBinomialDistribution[n, p] represents the Negative Binomial
  74. distribution for failure count n and probability p."
  75.  
  76. PoissonDistribution::usage =
  77. "PoissonDistribution[mu] represents the Poisson distribution with mean mu."
  78.  
  79.     (*   Functions  *)
  80.  
  81. (* Introduce usage messages for new functions. *)
  82.  
  83. Domain::usage =
  84. "Domain[distribution] gives the domain of the specified distribution."
  85.  
  86. PDF::usage =
  87. "PDF[distribution, x] gives the probability density function
  88. of the specified statistical distribution evaluated at x."
  89.  
  90. CDF::usage =
  91. "CDF[distribution, x] gives cumulative distribution function of distribution
  92. at the point x, defined as the integral of the probability density of the
  93. the distribution from the lowest value in the domain to x."
  94.  
  95. CharacteristicFunction::usage =
  96. "CharacteristicFunction[distribution, t] gives the characteristic
  97. function of the specified statistical distribution as a function of
  98. the variable t."
  99.  
  100. (*
  101.   Extend usage messages for existing functions.  Look for the
  102.   indicated phrase to determine if this has already been done.
  103. *)
  104.  
  105. If[StringPosition[Mean::usage, "specified statistical distribution"] === {},
  106.  
  107. Mean::usage = Mean::usage <> " " <>
  108. "Mean[distribution] gives the mean of the specified statistical distribution.";
  109.  
  110. Variance::usage = Variance::usage <> " " <>
  111. "Variance[distribution] gives the variance of the specified
  112. statistical distribution.";
  113.  
  114. StandardDeviation::usage = StandardDeviation::usage <> " " <>
  115. "StandardDeviation[distribution] gives the standard deviation of the
  116. specified statistical distribution.";
  117.  
  118. Skewness::usage = Skewness::usage <> " " <>
  119. "Skewness[distribution] gives the coefficient of skewness of the
  120. specified statistical distribution.";
  121.  
  122. Kurtosis::usage = Kurtosis::usage <> " " <>
  123. "Kurtosis[distribution] gives the coefficient of kurtosis of the
  124. specified statistical distribution.";
  125.  
  126. KurtosisExcess::usage = KurtosisExcess::usage <> " " <>
  127. "KurtosisExcess[distribution] gives the kurtosis excess for the
  128. specified statistical distribution.";
  129.  
  130. Quantile::usage =
  131. "Quantile[distribution, q] gives the q-th quantile of the specified
  132. statistical distribution."
  133. ]
  134.  
  135. (*
  136.   Extend the usage message of Random.  Look for the indicated phrase
  137.   to determine if this has already been done.
  138. *)
  139.  
  140. If[StringPosition[Random::usage, "specified statistical distribution"] === {},
  141.  
  142. Random::usage = Random::usage <> " " <>
  143. "Random[distribution] gives a random number with the specified
  144. statistical distribution."
  145. ]
  146.  
  147. Begin["`Private`"]
  148.  
  149. (* Bernoulli Distribution *)
  150.  
  151. BernoulliQ[p_] := If[0 <= N[p] <= 1, True,
  152.     Message[BernoulliDistribution::bern, p]; False, True]
  153.  
  154. BernoulliDistribution::bern =
  155. "Bernoulli distribution parameter `1` is expected to be a probability
  156. between 0 and 1."
  157.  
  158. BernoulliDistribution/: Domain[BernoulliDistribution[___]]:={0, 1}
  159.  
  160. BernoulliDistribution/: PDF[BernoulliDistribution[p_], x_] :=
  161.     Which[x == 0, 1-p, x == 1, p, NumberQ[x], 0] /; BernoulliQ[p]
  162.  
  163. BernoulliDistribution/:
  164.         CDF[BernoulliDistribution[p_], x_] :=
  165.     Which[Negative[x], 0, x >= 1, 1, True,  1-p]  /; BernoulliQ[p]
  166.  
  167. BernoulliDistribution/:
  168.     Mean[BernoulliDistribution[p_]] := p /; BernoulliQ[p]
  169.  
  170. BernoulliDistribution/:
  171.     Variance[BernoulliDistribution[p_]] := p (1-p) /; BernoulliQ[p]
  172.  
  173. BernoulliDistribution/: StandardDeviation[BernoulliDistribution[p_]] :=
  174.     Sqrt[p (1-p)] /; BernoulliQ[p]
  175.  
  176. BernoulliDistribution/: Skewness[BernoulliDistribution[p_]] :=
  177.     (1 - 2 p)/Sqrt[p (1 - p)] /; BernoulliQ[p]
  178.  
  179. BernoulliDistribution/: Kurtosis[BernoulliDistribution[p_]] :=
  180.     3 + (1 - 6 p (1-p))/(p (1-p)) /; BernoulliQ[p]
  181.  
  182. BernoulliDistribution/: KurtosisExcess[BernoulliDistribution[p_]] :=
  183.     (1 - 6 p (1-p))/(p (1-p)) /; BernoulliQ[p]
  184.  
  185. BernoulliDistribution/: CharacteristicFunction[
  186.         BernoulliDistribution[p_], t_] := 1 - p + p Exp[I t] /; BernoulliQ[p]
  187.  
  188. BernoulliDistribution/: Quantile[BernoulliDistribution[p_], q_] :=
  189.     Which[q <= 1-p, 0, q > 1-p, 1] /; QuantileQ[q] && BernoulliQ[p]
  190.  
  191. BernoulliDistribution/: Random[BernoulliDistribution[p_]] :=
  192.         If[Random[] > p, 0, 1] /; BernoulliQ[p]
  193.  
  194. (* Binomial Distribution *)
  195.  
  196. BinomialPQ[n_, p_] := And[
  197.     If[N[n] > 0, If[IntegerQ[n], True,
  198.         Message[BinomialDistribution::trials, n]; False],
  199.         Message[BinomialDistribution::trials, n]; False, True],
  200.     If[0 <= N[p] <= 1, True,
  201.         Message[BinomialDistribution::probability, p]; False, True]
  202. ]
  203.  
  204. BinomialDistribution::trials =
  205. "The parameter `1` is expected to be a positive integer."
  206.         
  207. BinomialDistribution::probability =
  208. "The parameter `1` is expected to be a probability between 0 and 1."
  209.  
  210. BinomialDistribution/: Domain[BinomialDistribution[n_, p_]]:= Range[n+1]-1 /;
  211.                 IntegerQ[n] && 1<=n<=2 && BinomialPQ[n,p]
  212.                     
  213. BinomialDistribution/: Domain[BinomialDistribution[n_, p_]]:=
  214.     {0,1," ..." ,n}        /; BinomialPQ[n,p]
  215.  
  216. BinomialDistribution/: PDF[BinomialDistribution[n_, p_], x_] :=
  217.     If[0 <= x <= n,
  218.         If[IntegerQ[x], Binomial[n, x] p^x (1-p)^(n-x), 0], 0] /;
  219.         BinomialPQ[n, p]
  220.  
  221. BinomialDistribution/: CDF[BinomialDistribution[n_, p_], x_] :=
  222.     Which[
  223.         Negative[x], 0,
  224.         x >= n, 1,
  225.         True, BetaRegularized[1 - p, n-Floor[x], Floor[x]+1]
  226.     ] /; BinomialPQ[n, p]
  227.          
  228. BinomialDistribution/: Mean[BinomialDistribution[n_, p_]] := n p /;
  229.         BinomialPQ[n, p]
  230.  
  231. BinomialDistribution/: Variance[BinomialDistribution[n_, p_]] :=
  232.     n p (1-p) /; BinomialPQ[n, p]
  233.  
  234. BinomialDistribution/:
  235.         StandardDeviation[BinomialDistribution[n_, p_]] :=
  236.     Sqrt[n p (1-p)] /; BinomialPQ[n, p]
  237.  
  238. BinomialDistribution/: Skewness[BinomialDistribution[n_, p_]] :=
  239.     (1 - 2 p)/Sqrt[n p (1-p)] /; BinomialPQ[n, p]
  240.  
  241. BinomialDistribution/: Kurtosis[BinomialDistribution[n_, p_]] :=
  242.     3 + (1 - 6 p (1-p))/(n p (1-p)) /; BinomialPQ[n, p]
  243.  
  244. BinomialDistribution/: KurtosisExcess[BinomialDistribution[n_, p_]] :=
  245.     (1 - 6 p (1-p))/(n p (1-p)) /; BinomialPQ[n, p]
  246.  
  247. BinomialDistribution/: CharacteristicFunction[
  248.         BinomialDistribution[n_, p_], t_] :=
  249.     (1 - p + p Exp[I t])^n /; BinomialPQ[n, p]
  250.  
  251. BinomialDistribution/:
  252.         Quantile[BinomialDistribution[n_, p_], q_] :=
  253.     With[{result = iBinomialQuantile[n, p, q]},
  254.         result /; result =!= Fail] /; 
  255.     NumberQ[n] && NumberQ[N[p]] && BinomialPQ[n, p] && QuantileQ[q]
  256.  
  257. iBinomialQuantile[n_, p_, q_] := Module[{low, high, mid},
  258.   If[q==1, Return[n]];
  259.   If[N[(1-p)^n] < q,
  260.     low = 0; high = n;
  261.     While[high - low > 1,
  262.         mid = Floor[(high+low)/2];
  263.         If[N[CDF[BinomialDistribution[n, p], mid]] < q,
  264.             low = mid,
  265.             high = mid,
  266.             high = Fail; Break[]
  267.         ]
  268.     ]; high,
  269.     0,
  270.     Fail
  271.   ]
  272. ]
  273.  
  274. BinomialDistribution/: Random[BinomialDistribution[n_, p_]] :=
  275.     Quantile[BinomialDistribution[n, p], Random[]]
  276.  
  277. (* Discrete Uniform Distribution *)
  278.  
  279. DiscreteUniformQ[n_] := If[N[n] > 0, If[IntegerQ[n], True,
  280.         Message[DiscreteUniformDistribution::count, n]; False],
  281.         Message[DiscreteUniformDistribution::count, n]; False, True]
  282.  
  283. DiscreteUniformDistribution::count =
  284. "The parameter `1` is expected to be a positive integer."
  285.  
  286. DiscreteUniformDistribution/: Domain[DiscreteUniformDistribution[n_]]:=
  287.     Range[n] /; IntegerQ[n] && 1<=n<=3 && DiscreteUniformQ[n]
  288.  
  289. DiscreteUniformDistribution/:
  290.     Domain[DiscreteUniformDistribution[n_]]:={1,2, "..." ,n} /;
  291.                         DiscreteUniformQ[n]
  292.  
  293. DiscreteUniformDistribution/:
  294.         PDF[DiscreteUniformDistribution[n_], x_] :=
  295.     If[1 <= x <= n, If[IntegerQ[x], 1/n, 0], 0] /; DiscreteUniformQ[n]
  296.  
  297. DiscreteUniformDistribution/: CDF[DiscreteUniformDistribution[n_], x_] :=
  298.     Which[
  299.         Negative[x], 0,
  300.         x >= n, 1,
  301.         True, Max[0, Min[1, Floor[x]/n]]
  302.     ] /;    DiscreteUniformQ[n]
  303.  
  304. DiscreteUniformDistribution/:
  305.         Mean[DiscreteUniformDistribution[n_]] :=
  306.     (n + 1) / 2 /; DiscreteUniformQ[n]
  307.  
  308. DiscreteUniformDistribution/:
  309.         Variance[DiscreteUniformDistribution[n_]] :=
  310.     (n^2 - 1)/12 /; DiscreteUniformQ[n]
  311.  
  312. DiscreteUniformDistribution/:
  313.         StandardDeviation[DiscreteUniformDistribution[n_]] :=
  314.     Sqrt[(n^2-1)/12] /; DiscreteUniformQ[n]
  315.  
  316. DiscreteUniformDistribution/: Skewness[
  317.         DiscreteUniformDistribution[n_]] := 0 /; DiscreteUniformQ[n]
  318.  
  319. DiscreteUniformDistribution/:
  320.         Kurtosis[DiscreteUniformDistribution[n_]] :=
  321.     3/5 (3 - 4/(n^2-1)) /; DiscreteUniformQ[n]
  322.  
  323. DiscreteUniformDistribution/: KurtosisExcess[
  324.         DiscreteUniformDistribution[n_]] :=
  325.     -(12/(n^2-1) + 6)/5 /; DiscreteUniformQ[n]
  326.  
  327. DiscreteUniformDistribution/: CharacteristicFunction[
  328.         DiscreteUniformDistribution[n_], t_] :=
  329.     Sum[Exp[I t k],{k,1,n}]/n    /; DiscreteUniformQ[n]
  330.  
  331. DiscreteUniformDistribution/: Quantile[
  332.         DiscreteUniformDistribution[n_], q_] :=
  333.     Max[Ceiling[N[n q]], 1] /; DiscreteUniformQ[n] && QuantileQ[q]
  334.  
  335. DiscreteUniformDistribution/: Random[DiscreteUniformDistribution[n_]] :=
  336.     Random[Integer, {1, n}] /; DiscreteUniformQ[n]
  337.  
  338. (* Geometric Distribution *)
  339.  
  340. GeometricQ[p_] := If[0 <= N[p] <= 1, True,
  341.         Message[GeometricDistribution::probability, p]; False, True]
  342.  
  343. GeometricDistribution::probability =
  344. "The parameter `1` is expected to be a probability between 0 and 1."
  345.  
  346. GeometricDistribution/:
  347.     Domain[GeometricDistribution[___]]:= {0,1,2, "..."}
  348.  
  349. GeometricDistribution/:
  350.         PDF[GeometricDistribution[p_], x_] :=
  351.     If[!Negative[x], If[IntegerQ[x], p (1-p)^x, 0], 0] /; GeometricQ[p]
  352.  
  353. GeometricDistribution/: 
  354.         CDF[GeometricDistribution[p_], x_] :=
  355.     Which[
  356.         Negative[x], 0,
  357.         x == Infinity, 1,
  358.         True, 1 - (1-p)^Floor[x + 1]
  359.     ] /;    GeometricQ[p]
  360.  
  361. GeometricDistribution/:
  362.     Mean[GeometricDistribution[p_]] := 1/p - 1 /; GeometricQ[p]
  363.  
  364. GeometricDistribution/: Variance[GeometricDistribution[p_]] :=
  365.     (1 - p) / p^2 /; GeometricQ[p]
  366.  
  367. GeometricDistribution/: StandardDeviation[GeometricDistribution[p_]] :=
  368.     Sqrt[1 - p] / p /; GeometricQ[p]
  369.  
  370. GeometricDistribution/: Skewness[GeometricDistribution[p_]] :=
  371.     (2 - p)/Sqrt[1 - p] /; GeometricQ[p]
  372.  
  373. GeometricDistribution/: Kurtosis[GeometricDistribution[p_]] :=
  374.     (p^2 + 6 - 6 p) / (1 - p) + 3 /; GeometricQ[p]
  375.  
  376. GeometricDistribution/: KurtosisExcess[GeometricDistribution[p_]] :=
  377.     (p^2 + 6 - 6 p) / (1 - p) /; GeometricQ[p]
  378.  
  379. GeometricDistribution/: CharacteristicFunction[
  380.         GeometricDistribution[p_], t_] :=
  381.     p / (1 - (1-p) Exp[I t]) /; GeometricQ[p]
  382.  
  383. GeometricDistribution/: Quantile[GeometricDistribution[p_], q_] :=
  384.     With[{result = iGeometricQuantile[p, q]},
  385.         result /; result =!= Fail
  386.     ] /; GeometricQ[p] && QuantileQ[q]
  387.  
  388. iGeometricQuantile[p_, q_] :=
  389.    (If[q == 1, Return[Infinity]];
  390.     If[N[p] < q,
  391.     Ceiling[N[Log[1-q]/Log[1-p]]] - 1,
  392.     0,
  393.     Fail
  394.     ])
  395.  
  396. GeometricDistribution/: Random[GeometricDistribution[p_]] :=
  397.     Quantile[GeometricDistribution[p], Random[]]
  398.  
  399. (* Hypergeometric Distribution *)
  400.  
  401. HypergeometricPQ[n_, nsucc_, ntot_] := And[
  402. If[N[n] >= 0,
  403.     If[IntegerQ[n], True,
  404.         Message[HypergeometricDistribution::posint, n]; False],
  405.     Message[HypergeometricDistribution::posint, n]; False,
  406.     True],
  407. If[N[nsucc] >= 0,
  408.     If[IntegerQ[nsucc], True,
  409.         Message[HypergeometricDistribution::posint, nsucc]; False],
  410.     Message[HypergeometricDistribution::posint, nsucc]; False,
  411.     True],
  412. If[N[ntot] >= 0,
  413.     If[IntegerQ[ntot], True,
  414.         Message[HypergeometricDistribution::posint, ntot]; False],
  415.     Message[HypergeometricDistribution::posint, ntot]; False,
  416.     True],
  417. If[N[n] <= N[ntot], True,
  418.     Message[HypergeometricDistribution::ntoobig, n, ntot]; False, True],
  419. If[N[nsucc] <= N[ntot], True,
  420.     Message[HypergeometricDistribution::nsucctoobig, nsucc, ntot]; False, True]
  421. ]
  422.  
  423. HypergeometricDistribution::posint =
  424. "The parameter `1` is expected to be a positive integer."
  425.  
  426. HypergeometricDistribution::ntoobig =
  427. "The number of trials `1` without replacement must be less than
  428. the total number `2`."
  429.  
  430. HypergeometricDistribution::nsucctoobig =
  431. "The number of possible successes `1` must be less than the total
  432. number `2`."
  433.  
  434. HypergeometricDistribution/:
  435.     Domain[HypergeometricDistribution[n_, nsucc_, ntot_]]:=
  436.         {Max[0,n-ntot+nsucc], "...", Min[n,nsucc]} /;
  437.             HypergeometricPQ[n,nsucc,ntot]
  438.  
  439. HypergeometricDistribution/: PDF[
  440.         HypergeometricDistribution[n_, nsucc_, ntot_], x_] :=
  441.     If[Max[0,n-ntot+nsucc] <= x <= Min[n,nsucc],
  442.         If[IntegerQ[x],
  443.             Binomial[nsucc, x] Binomial[ntot-nsucc, n-x]/
  444.                 Binomial[ntot, n], 0],
  445.         0] /; HypergeometricPQ[n, nsucc, ntot]
  446.  
  447. HypergeometricDistribution/:
  448.         CDF[HypergeometricDistribution[n_, nsucc_, ntot_], x_] :=
  449.     Which[
  450.     x < Max[0,n-ntot+nsucc], 0,
  451.     x >= Min[n,nsucc], 1,
  452.     True, Sum[PDF[HypergeometricDistribution[n,nsucc,ntot],k],
  453.             {k,Max[0,n-ntot+nsucc],Min[x,n,nsucc]}]
  454.     ] /; HypergeometricPQ[n, nsucc, ntot]
  455.  
  456. HypergeometricDistribution/:
  457.         Mean[HypergeometricDistribution[n_, nsucc_, ntot_]] :=
  458.     n nsucc / ntot /; HypergeometricPQ[n, nsucc, ntot]
  459.  
  460. HypergeometricDistribution/:
  461.         Variance[HypergeometricDistribution[n_, nsucc_, ntot_]] :=
  462.     (ntot-n)/(ntot-1) n nsucc/ntot (1-nsucc/ntot) /;
  463.         HypergeometricPQ[n, nsucc, ntot]
  464.  
  465. HypergeometricDistribution/: StandardDeviation[
  466.         HypergeometricDistribution[n_, nsucc_, ntot_]] :=
  467.     Sqrt[(ntot-n)/(ntot-1) n nsucc (ntot-nsucc)]/ntot /;
  468.         HypergeometricPQ[n, nsucc, ntot]
  469.  
  470. HypergeometricDistribution/: Skewness[
  471.         HypergeometricDistribution[n_, nsucc_, ntot_]] :=
  472.     (ntot-2 nsucc) (ntot-2 n) Sqrt[ntot-1]/
  473.         ((ntot-2) Sqrt[n nsucc (ntot-nsucc) (ntot-n)]) /;
  474.             HypergeometricPQ[n, nsucc, ntot]
  475.  
  476. HypergeometricDistribution/: Kurtosis[
  477.         HypergeometricDistribution[n_, nsucc_, ntot_]] :=
  478.     ntot^2 (ntot-1)/((ntot-2)(ntot-3) n nsucc (ntot-nsucc)(ntot-n))*
  479.       (ntot (ntot+1) - 6 n (ntot-n) + 3 nsucc (ntot-nsucc)/ntot^2 *
  480.         (ntot^2 (n-2) - ntot n^2 + 6 n (ntot-n))) /;
  481.             HypergeometricPQ[n, nsucc, ntot]
  482.  
  483. HypergeometricDistribution/: KurtosisExcess[
  484.         HypergeometricDistribution[n_, nsucc_, ntot_]] :=
  485.     Kurtosis[HypergeometricDistribution[n, nsucc, ntot]] - 3 /;
  486.         HypergeometricPQ[n, nsucc, ntot]
  487.  
  488. HypergeometricDistribution/: CharacteristicFunction[
  489.         HypergeometricDistribution[n_, nsucc_, ntot_], t_] :=
  490.     ((ntot-n)! (ntot-nsucc)!/ntot!) *
  491.     Hypergeometric2F1[-n, -nsucc, ntot-nsucc-n+1, E^(I t)] /;
  492.         HypergeometricPQ[n, nsucc, ntot]
  493.  
  494. HypergeometricDistribution/: Quantile[
  495.         HypergeometricDistribution[n_, nsucc_, ntot_], q_] :=
  496.     With[{result = iHypergeometricQuantile[n, nsucc, ntot, q]},
  497.         result /; result =!= Fail
  498.     ] /; HypergeometricPQ[n, nsucc, ntot] && QuantileQ[q]
  499.  
  500. iHypergeometricQuantile[n_, nsucc_, ntot_, q_] :=
  501. Module[{kk = 0, sum, count = Max[0,n-ntot+nsucc]},
  502.     If[q == 1, Return[Min[nsucc, n]]];
  503.     sum = PDF[HypergeometricDistribution[n, nsucc, ntot], count];
  504.     If[NumberQ[sum],
  505.         While[sum < q,
  506.             count++;
  507.             sum += Binomial[nsucc, count] *
  508.                     Binomial[ntot-nsucc, n-count]/
  509.                         Binomial[ntot, n]
  510.         ],
  511.     (* else not a number *)
  512.         count = Fail
  513.     ];
  514.     count
  515. ]
  516.  
  517. HypergeometricDistribution/:
  518.         Random[HypergeometricDistribution[n_, nsucc_, ntot_]] :=
  519.     Quantile[HypergeometricDistribution[n,nsucc,ntot], Random[]] 
  520.  
  521. (* Logarithmic Series Distribution *)
  522.  
  523. LogSeriesQ[theta_] :=
  524.     If[0 <= N[theta] < 1, True,
  525.         Message[LogSeriesDistribution::outofrange, theta]; False,
  526.         True]
  527.  
  528. LogSeriesDistribution::outofrange =
  529. "Parameter `1` is expected to be less than 1 and non-negative."
  530.  
  531. LogSeriesDistribution/: Domain[LogSeriesDistribution[___]]:= {1,2,3, "..."}
  532.  
  533. LogSeriesDistribution/: PDF[LogSeriesDistribution[theta_], x_] :=
  534.     If[Positive[x], If[IntegerQ[x], theta^x / (-x Log[1-theta]), 0], 0] /;
  535.                             LogSeriesQ[theta]
  536.  
  537. LogSeriesDistribution/:
  538.         CDF[LogSeriesDistribution[theta_], x_] :=
  539.     Which[
  540.       x < 1, 0,
  541.       x == Infinity, 1,
  542.       True,
  543.        Sum[theta^k / k, {k, 1, x}] / (-Log[1-theta])
  544.     ]    /; LogSeriesQ[theta]
  545.  
  546. LogSeriesDistribution/: Mean[LogSeriesDistribution[theta_]] :=
  547.     theta/(1 - theta) / (-Log[1-theta]) /; LogSeriesQ[theta]
  548.  
  549. LogSeriesDistribution/: Variance[LogSeriesDistribution[theta_]] :=
  550.   Module[{alpha=-1/Log[1-theta]},
  551.     alpha theta (1-alpha theta)/(1-theta)^2
  552.   ]    /; LogSeriesQ[theta]
  553.  
  554. LogSeriesDistribution/: StandardDeviation[LogSeriesDistribution[theta_]] :=
  555.   Module[{alpha=-1/Log[1-theta]},
  556.     Sqrt[alpha theta (1-alpha theta)]/(1-theta)
  557.   ]    /; LogSeriesQ[theta]
  558.  
  559. LogSeriesDistribution/: Skewness[LogSeriesDistribution[theta_]] :=
  560.   Module[{alpha=-1/Log[1-theta]},
  561.     (1 + theta - 3 alpha theta + 2 alpha^2 theta^2)/(1-alpha theta)/
  562.     Sqrt[alpha theta(1-alpha theta)]
  563.   ]    /; LogSeriesQ[theta]
  564.  
  565. LogSeriesDistribution/: Kurtosis[LogSeriesDistribution[theta_]] :=
  566.   Module[{alpha=-1/Log[1-theta]},
  567.     (1 + 4 theta + theta^2 - 4 alpha theta (1+theta) + 6 alpha^2 theta^2 -
  568.     3 alpha^3 theta^3)/(alpha theta)/(1-alpha theta)^2
  569.   ]     /; LogSeriesQ[theta]
  570.  
  571. LogSeriesDistribution/: KurtosisExcess[LogSeriesDistribution[theta_]] :=
  572.     Kurtosis[LogSeriesDistribution[theta]] - 3 /; LogSeriesQ[theta]
  573.  
  574. LogSeriesDistribution/:
  575.         CharacteristicFunction[LogSeriesDistribution[theta_], t_] :=
  576.     Log[1 - theta E^(t I)] / Log[1-theta]    /; LogSeriesQ[theta]
  577.  
  578. LogSeriesDistribution/:
  579.         Quantile[LogSeriesDistribution[theta_], q_] :=
  580.     With[{result = iLogSeriesQuantile[theta, q]},
  581.         result /; result =!= Fail] /;
  582.              LogSeriesQ[theta] && QuantileQ[q]
  583.  
  584. iLogSeriesQuantile[theta_, q_] := Module[{high, low, mid},
  585.     If[q == 1, Return[Infinity]];
  586.     If[N[theta/(-Log[1-theta])] < q,
  587.         high = 2;
  588.         While[N[CDF[LogSeriesDistribution[theta], high]] < q, high *= 2];
  589.         low = high/2;
  590.         While[high - low > 1,
  591.             mid = (high+low)/2;
  592.             If[N[CDF[LogSeriesDistribution[theta], mid]] < q,
  593.                 low = mid,
  594.                 high = mid,
  595.                 low = high
  596.             ]
  597.         ]; high,
  598.     (* else *)
  599.         1,
  600.     (* indeterminate *)
  601.         Fail
  602.     ]
  603. ]
  604.  
  605. LogSeriesDistribution/: Random[LogSeriesDistribution[theta_]] :=
  606.        Quantile[LogSeriesDistribution[theta],Random[]] /; LogSeriesQ[theta]
  607.  
  608. (* Negative Binomial Distribution *)
  609.  
  610. NegativeBinomialPQ[n_, p_] := And[
  611.     If[N[n] > 0,
  612.         If[IntegerQ[n], True,
  613.             Message[NegativeBinomialDistribution::trials, n]; False],
  614.         Message[NegativeBinomialDistribution::trials, n]; False, True],
  615.     If[0 < N[p] <= 1, True,
  616.         Message[NegativeBinomialDistribution::probability, p]; False, True]
  617. ]
  618.  
  619. NegativeBinomialDistribution::trials =
  620. "The parameter `1` is expected to be a positive integer."
  621.         
  622. NegativeBinomialDistribution::probability =
  623. "The parameter `1` is expected to be a non-zero probability
  624. between 0 and 1."
  625.  
  626. NegativeBinomialDistribution/:
  627.     Domain[NegativeBinomialDistribution[___]]:={0,1,2, "..."}
  628.  
  629. NegativeBinomialDistribution/:
  630.         PDF[NegativeBinomialDistribution[n_, p_], x_] :=
  631.     If[!Negative[x],
  632.         If[IntegerQ[x], Binomial[x+n-1, n-1] p^n (1-p)^x, 0],
  633.         0] /; NegativeBinomialPQ[n, p]
  634.  
  635. NegativeBinomialDistribution/:
  636.         CDF[NegativeBinomialDistribution[n_, p_], x_] :=
  637.     Which[
  638.       Negative[x], 0,
  639.       x == Infinity, 1,
  640.       True,
  641.     BetaRegularized[p, n, Floor[x] + 1]
  642.     ]    /; NegativeBinomialPQ[n, p]
  643.  
  644. NegativeBinomialDistribution/:
  645.         Mean[NegativeBinomialDistribution[n_, p_]] :=
  646.     n (1-p)/p /; NegativeBinomialPQ[n, p]
  647.  
  648. NegativeBinomialDistribution/: Variance[
  649.                 NegativeBinomialDistribution[n_, p_]] :=
  650.         n (1-p) / p^2 /; NegativeBinomialPQ[n, p]
  651.  
  652. NegativeBinomialDistribution/: StandardDeviation[
  653.                 NegativeBinomialDistribution[n_, p_]] :=
  654.         Sqrt[n (1-p)] / p /; NegativeBinomialPQ[n, p]
  655.  
  656. NegativeBinomialDistribution/: Skewness[
  657.                 NegativeBinomialDistribution[n_, p_]] :=
  658.         (2 - p) / Sqrt[n (1-p)] /; NegativeBinomialPQ[n, p]
  659.  
  660. NegativeBinomialDistribution/:
  661.         Kurtosis[NegativeBinomialDistribution[n_, p_]] :=
  662.     3 + (p^2 + 6 - 6 p) / (n (1-p)) /; NegativeBinomialPQ[n, p]
  663.  
  664. NegativeBinomialDistribution/:
  665.         KurtosisExcess[NegativeBinomialDistribution[n_, p_]] :=
  666.     (p^2 + 6 - 6 p) / (n (1-p)) /; NegativeBinomialPQ[r, p]
  667.  
  668. NegativeBinomialDistribution/: CharacteristicFunction[
  669.         NegativeBinomialDistribution[n_, p_], t_] :=
  670.     (p/(1 - (1-p) Exp[I t]))^n /; NegativeBinomialPQ[n, p]
  671.  
  672. NegativeBinomialDistribution/:
  673.         Quantile[NegativeBinomialDistribution[n_, p_], q_] :=
  674.     With[{result = iNegativeBinomialQuantile[n, p, q]},
  675.         result /; result =!= Fail] /;
  676.             NegativeBinomialPQ[n, p] && QuantileQ[q]
  677.  
  678. iNegativeBinomialQuantile[n_, p_, q_] := Module[{high, low, mid},
  679.     If[q == 1, Return[Infinity]];
  680.     If[N[p^n] < q,
  681.         high = 1;
  682.         While[N[CDF[NegativeBinomialDistribution[n, p], high]] < q, high *= 2];
  683.         low = high/2;
  684.         While[high - low > 1,
  685.             mid = (high+low)/2;
  686.             If[N[CDF[NegativeBinomialDistribution[n, p], mid]] < q,
  687.                 low = mid,
  688.                 high = mid,
  689.                 low = high
  690.             ]
  691.         ]; high,
  692.     (* else *)
  693.         0,
  694.     (* indeterminate *)
  695.         Fail
  696.     ]
  697. ]
  698.  
  699. NegativeBinomialDistribution/: Random[
  700.                 NegativeBinomialDistribution[r_, p_]] :=
  701.        Quantile[NegativeBinomialDistribution[r,p],Random[]]
  702.  
  703. (* Poisson Distribution *)
  704.  
  705. PoissonPQ[mu_] := If[Positive[mu], True,
  706.     Message[PoissonDistribution::parameter, mu]; False, True]
  707.  
  708. PoissonDistribution::parameter =
  709. "Parameter `1` is expected to be greater than zero."
  710.  
  711. PoissonDistribution/: Domain[PoissonDistribution[___]]:={0,1,2, "..."}
  712.  
  713. PoissonDistribution/: PDF[PoissonDistribution[mu_], x_] :=
  714.     If[!Negative[x], If[IntegerQ[x], Exp[-mu] mu^x / x!, 0], 0] /;
  715.         PoissonPQ[mu]
  716.  
  717. PoissonDistribution/: CDF[PoissonDistribution[mu_], x_] :=
  718.     Which[
  719.       Negative[x], 0,
  720.       x == Infinity, 1,
  721.       True,
  722.     GammaRegularized[Floor[x] + 1, mu, Infinity]
  723.     ]    /; PoissonPQ[mu]
  724.  
  725. PoissonDistribution/: Mean[PoissonDistribution[mu_]] :=
  726.     mu /; PoissonPQ[mu]
  727.  
  728. PoissonDistribution/: Variance[PoissonDistribution[mu_]] :=
  729.     mu /; PoissonPQ[mu]
  730.  
  731. PoissonDistribution/: StandardDeviation[PoissonDistribution[mu_]] :=
  732.     Sqrt[mu] /; PoissonPQ[mu]
  733.  
  734. PoissonDistribution/: Skewness[PoissonDistribution[mu_]] :=
  735.     1/Sqrt[mu] /; PoissonPQ[mu]
  736.  
  737. PoissonDistribution/: Kurtosis[PoissonDistribution[mu_]] :=
  738.     3 + 1/mu /; PoissonPQ[mu]
  739.  
  740. PoissonDistribution/: KurtosisExcess[PoissonDistribution[mu_]] :=
  741.     1/mu /; PoissonPQ[mu]
  742.  
  743. PoissonDistribution/:
  744.         CharacteristicFunction[PoissonDistribution[mu_], t_] :=
  745.     Exp[mu (Exp[I t] - 1)] /; PoissonPQ[mu]
  746.  
  747. PoissonDistribution/: Quantile[PoissonDistribution[mu_], q_] :=
  748.     With[{result = iPoissonQuantile[mu, q]},
  749.         result /; result =!= Fail
  750.     ] /; PoissonPQ[mu] && QuantileQ[q]
  751.  
  752. iPoissonQuantile[mu_, q_] := Module[{high, low, mid},
  753.     If[q == 1, Return[Infinity]];
  754.     If[N[Exp[-mu]] < q,
  755.         high = 1;
  756.         While[N[CDF[PoissonDistribution[mu], high]] < q, high *= 2];
  757.         low = high/2;
  758.         While[high - low > 1,
  759.             mid = (high+low)/2;
  760.             If[N[CDF[PoissonDistribution[mu], mid]] < q,
  761.                 low = mid,
  762.                 high = mid,
  763.                 low = high
  764.             ]
  765.         ]; high,
  766.     (* else *)
  767.         0,
  768.     (* indeterminate *)
  769.         Fail
  770.     ]
  771. ]
  772.  
  773. PoissonDistribution/: Random[PoissonDistribution[mu_]] :=
  774.        Quantile[PoissonDistribution[mu],Random[]]
  775.  
  776. End[]
  777.  
  778. SetAttributes[ BernoulliDistribution , ReadProtected];
  779. SetAttributes[ BinomialDistribution , ReadProtected];
  780. SetAttributes[ DiscreteUniformDistribution , ReadProtected];
  781. SetAttributes[  GeometricDistribution, ReadProtected];
  782. SetAttributes[ HypergeometricDistribution , ReadProtected];
  783. SetAttributes[  LogSeriesDistribution, ReadProtected];
  784. SetAttributes[ NegativeBinomialDistribution , ReadProtected];
  785. SetAttributes[  PoissonDistribution, ReadProtected];
  786.  
  787. Protect[BernoulliDistribution, BinomialDistribution,
  788.     DiscreteUniformDistribution, GeometricDistribution,
  789.     HypergeometricDistribution, LogSeriesDistribution,
  790.     NegativeBinomialDistribution, PoissonDistribution];
  791.  
  792. EndPackage[]
  793.  
  794.