
(*==============================================================================

     distributions.tl

     This file is part of the TREEFINDER software.


                                    COPYRIGHT (C) 1997-2008 BY GANGOLF JOBB

==============================================================================*)
(*

 Synopsis:

  PDF[distribution,x]
  CDF[distribution,x]
  Quantile[distribution,p]
  Random[distribution]

*)
(*============================================================================*)

Create[SampleD], Protect[SampleD], (* SampleD[_data] *)
Create[NormalD], Protect[NormalD], (* NormalD[mu,sigma] *)
Create[StudentTD], Protect[StudentTD], (* StudentTD[n] *)
Create[ChiSquareD], Protect[ChiSquareD], (* ChiSquareD[n] *)
Create[FRatioD], Protect[FRatioD], (* FRatioD[n1,n2] *)
Create[BetaD], Protect[BetaD], (* BetaD[a,b] *)
Create[CauchyD], Protect[CauchyD], (* CauchyD[a,b] *)
Create[ChiD], Protect[ChiD], (* ChiD[n] *)
Create[ExponentialD], Protect[ExponentialD], (* ExponentialD[lambda] *)
Create[ExtremeValueD], Protect[ExtremeValueD], (* ExtremeValueD[alpha,beta] *)
Create[GammaD], Protect[GammaD], (* GammaD[alpha,beta] *)
Create[HalfNormalD], Protect[HalfNormalD], (* HalfNormalD[theta] *)
Create[LaplaceD], Protect[LaplaceD], (* LaplaceD[mu,beta] *)
Create[LogisticD], Protect[LogisticD], (* LogisticD[mu,beta] *)
Create[LogNormalD], Protect[LogNormalD], (* LogNormalD[mu,sigma] *)
Create[ParetoD], Protect[ParetoD], (* ParetoD[k,alpha] *)
Create[RayleighD], Protect[RayleighD], (* RayleighD[sigma] *)
Create[UniformD], Protect[UniformD], (* UniformD[min,max] *)
Create[WeibullD], Protect[WeibullD], (* WeibullD[alpha,beta] *)
Create[BernoulliD], Protect[BernoulliD], (* BernoulliD[p] *)
Create[BinomialD], Protect[BinomialD], (* BinomialD[n,p] *)
Create[DiscreteUniformD], Protect[DiscreteUniformD], (* DiscreteUniformD[n] *)
Create[GeometricD], Protect[GeometricD], (* GeometricD[p] *)
Create[HypergeometricD], Protect[HypergeometricD], (* HypergeometricD[n,nsucc,ntot] *)
Create[LogSeriesD], Protect[LogSeriesD], (* LogSeriesD[theta] *)
Create[NegativeBinomialD], Protect[NegativeBinomialD], (* NegativeBinomialD[n,p] *)
Create[PoissonD], Protect[PoissonD], (* PoissonD[mu] *)

(*============================================================================*)

PDF[distribution,x]:=<<

 If[!IsNumber[x],Return[ Put[{distribution,x},0,'PDF] ]];

 If[

  distribution==SampleD[_data],
  Scope[ {data:=Get[distribution,1],prob,step},
   If[ !IsNumericVector[data]||Length[data]<2 ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[Rtol[1e-13,IsInvariable[data]],Throw["PDF: SampleD: No variation in data."]];
   prob=0.0; step=1.0/(Length[data]-1); data=Table[{i,prob};prob+=step,(=i,Sort[data]=)];
   InterpolateDerivative[FunctionTable[data],x]
  ],

  distribution==NormalD[_mu,_sigma],
  Scope[ {mu:=Get[distribution,1],sigma:=Get[distribution,2]},
   If[ !IsNumber[mu]||!IsNumber[sigma] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[sigma<=0,Throw["PDF: NormalD: Second parameter is expected to be positive."]];
   Exp[-((x-mu)/sigma)^2/2]/(sigma*Sqrt[2*Pi])
  ],

  distribution==StudentTD[_n],
  Scope[ {n:=Get[distribution,1]},
   If[ !IsNumber[n] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[n<=0,Throw["PDF: StudentTD: Parameter is expected to be positive."]];
   1/(Sqrt[n]*Beta[n/2,1/2])*Sqrt[n/(n+x^2)]^(n+1)
  ],

  distribution==ChiSquareD[_n],
  Scope[ {n:=Get[distribution,1]},
   If[ !IsNumber[n] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[n<=0,Throw["PDF: ChiSquareD: Parameter is expected to be positive."]];
   If[x<0,Return[0]];
   x^(n/2-1)/(Exp[x/2]*Sqrt[2]^n*Gamma[n/2])
  ],

  distribution==FRatioD[_n1,_n2],
  Scope[ {n1:=Get[distribution,1],n2:=Get[distribution,2]},
   If[ !IsNumber[n1]||!IsNumber[n2] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[n1<=0||n2<=0,Throw["PDF: FRatioD: Parameters are expected to be positive."]];
   If[x<0,Return[0]];
   n1^(n1/2)*n2^(n2/2)*x^(n1/2-1)/((n2+n1*x)^((n1+n2)/2)*Beta[n1/2,n2/2])
  ],

  distribution==BetaD[_a,_b],
  Scope[ {a:=Get[distribution,1],b:=Get[distribution,2]},
   If[ !IsNumber[a]||!IsNumber[b] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[a<=0||b<=0,Throw["PDF: BetaD: Parameters are expected to be positive."]];
   If[x<0||x>1,Return[0]];
   x^(a-1)*(1-x)^(b-1)/Beta[a,b]
  ],

  distribution==CauchyD[_a,_b],
  Scope[ {a:=Get[distribution,1],b:=Get[distribution,2]},
   If[ !IsNumber[a]||!IsNumber[b] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[b<=0,Throw["PDF: CauchyD: Second parameter is expected to be positive."]];
   1/(b*Pi*(1+((x-a)/b)^2))
  ],

  distribution==ChiD[_n],
  Scope[ {n:=Get[distribution,1]},
   If[ !IsNumber[n] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[n<=0,Throw["PDF: ChiD: Parameter is expected to be positive."]];
   If[x<0,Return[0]];
   x^(n-1)*Exp[-x^2/2]/(2^(n/2-1)*Gamma[n/2])
  ],

  distribution==ExponentialD[_lambda],
  Scope[ {lambda:=Get[distribution,1]},
   If[ !IsNumber[lambda] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[lambda<=0,Throw["PDF: ExponentialD: Parameter is expected to be positive."]];
   If[x<0,Return[0]];
   lambda*Exp[-lambda*x]
  ],

  distribution==ExtremeValueD[_alpha,_beta],
  Scope[ {alpha:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[alpha]||!IsNumber[beta] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[beta<=0,Throw["PDF: ExtremeValueD: Second parameter is expected to be positive."]];
   Exp[-Exp[(alpha-x)/beta]+(alpha-x)/beta]/beta
  ],

  distribution==GammaD[_alpha,_beta],
  Scope[ {alpha:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[alpha]||!IsNumber[beta] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[alpha<=0||beta<=0,Throw["PDF: GammaD: Parameters are expected to be positive."]];
   If[x<0,Return[0]];
   x^(alpha-1)*Exp[-x/beta]/(beta^alpha*Gamma[alpha])
  ],

  distribution==HalfNormalD[_theta],
  Scope[ {theta:=Get[distribution,1]},
   If[ !IsNumber[theta] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[theta<=0,Throw["PDF: HalfNormalD: Parameter is expected to be positive."]];
   If[x<0,Return[0]];
   2*theta*Exp[- theta^2*x^2/Pi]/Pi
  ],

  distribution==LaplaceD[_mu,_beta],
  Scope[ {mu:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[mu]||!IsNumber[beta] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[beta<=0,Throw["PDF: LaplaceD: Second parameter is expected to be positive."]];
   (Exp[-Sign[x-mu]*(x-mu)/beta]/(2*beta))
  ],

  distribution==LogisticD[_mu,_beta],
  Scope[ {mu:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[mu]||!IsNumber[beta] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[beta<=0,Throw["PDF: LogisticD: Second parameter is expected to be positive."]];
   Exp[-(x-mu)/beta]/(beta*(1+Exp[-(x-mu)/beta])^2)
  ],

  distribution==LogNormalD[_mu,_sigma],
  Scope[ {mu:=Get[distribution,1],sigma:=Get[distribution,2]},
   If[ !IsNumber[mu]||!IsNumber[sigma] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[sigma<=0,Throw["PDF: LogNormalD: Second parameter is expected to be positive."]];
   If[x<=0,Return[0]];
   Exp[-(Ln[x]-mu)^2/(2*sigma^2)]/(Sqrt[2*Pi]*sigma*x)
  ],

  distribution==ParetoD[_k,_alpha],
  Scope[ {k:=Get[distribution,1],alpha:=Get[distribution,2]},
   If[ !IsNumber[k]||!IsNumber[alpha] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[k<=0||alpha<=0,Throw["PDF: ParetoD: Parameters are expected to be positive."]];
   If[x<0,Return[0]];
   alpha*k^alpha*x^(-alpha-1)
  ],

  distribution==RayleighD[_sigma],
  Scope[ {sigma:=Get[distribution,1]},
   If[ !IsNumber[sigma] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[sigma<=0,Throw["PDF: RayleighD: Parameter is expected to be positive."]];
   If[x<0,Return[0]];
   x*Exp[-x^2/(2*sigma^2)]/sigma^2
  ],

  distribution==UniformD[_min,_max],
  Scope[ {min:=Get[distribution,1],max:=Get[distribution,2]},
   If[ !IsNumber[min]||!IsNumber[max] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[max-min<=0,Throw["PDF: UniformD: Second parameter is expected to be greater than the first."]];
   (Sign[x-min]-Sign[x-max])/(2*(max-min))
  ],

  distribution==WeibullD[_alpha,_beta],
  Scope[ {alpha:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[alpha]||!IsNumber[beta] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[alpha<=0||beta<=0,Throw["PDF: WeibullD: Parameters are expected to be positive."]];
   If[x<0,Return[0]];
   alpha*x^(alpha-1)*Exp[-(x/beta)^alpha]/beta^alpha
  ],

  distribution==BernoulliD[_p],
  Scope[ {p:=Get[distribution,1]},
   If[ !IsNumber[p]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[p<0||p>1,Throw["PDF: BernoulliD: Parameter is expected to be a probability."]];
   If[x<0||x>1,Return[0]];
   If[x==0,1-p,p]
  ],

  distribution==BinomialD[_n,_p],
  Scope[ {n:=Get[distribution,1],p:=Get[distribution,2]},
   If[ !IsInteger[n]||!IsNumber[p]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[n<1,Throw["PDF: BinomialD: First parameter is expected to be a positive integer."]];
   If[p<0||p>1,Throw["PDF: BinomialD: Second parameter is expected to be a probability."]];
   If[x<0||x>n,Return[0]];
   Binomial[n,x]*p^x*(1-p)^(n-x)
  ],

  distribution==DiscreteUniformD[_n],
  Scope[ {n:=Get[distribution,1]},
   If[ !IsInteger[n]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[n<1,Throw["PDF: DiscreteUniformD: Parameter is expected to be a positive integer."]];
   If[x<1||x>n,Return[0]];
   1/n
  ],

  distribution==GeometricD[_p],
  Scope[ {p:=Get[distribution,1]},
   If[ !IsNumber[p]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[p<0||p>1,Throw["PDF: GeometricD: Parameter is expected to be a probability."]];
   If[x<0,Return[0]];
   p*(1-p)^x
  ],

  distribution==HypergeometricD[_n,_nsucc,_ntot],
  Scope[ {n:=Get[distribution,1],nsucc:=Get[distribution,2],ntot:=Get[distribution,3]},
   If[ !IsInteger[n]||!IsInteger[nsucc]||!IsInteger[ntot]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[n<0,Throw["PDF: HypergeometricD: First parameter is expected to be non-negative."]];
   If[nsucc<0,Throw["PDF: HypergeometricD: Second parameter is expected to be non-negative."]];
   If[ntot<1,Throw["PDF: HypergeometricD: Third parameter is expected to be positive."]];
   If[n>ntot,Throw["PDF: HypergeometricD: First parameter must be less or equal than the third."]];
   If[nsucc>ntot,Throw["PDF: HypergeometricD: Second parameter must be less or equal than the third."]];
   If[x<Max[{0,n-ntot+nsucc}]||x>Min[{n,nsucc}],Return[0]];
   (Binomial[nsucc,x]*Binomial[ntot-nsucc,n-x])/Binomial[ntot,n]
  ],

  distribution==LogSeriesD[_theta],
  Scope[ {theta:=Get[distribution,1]},
   If[ !IsNumber[theta]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[theta<=0||theta>=1,Throw["PDF: LogSeriesD: Parameter is expected to be less than 1 and non-negative."]];
   If[x<=0,Return[0]];
   theta^x/(-x*Ln[1-theta])
  ],

  distribution==NegativeBinomialD[_n,_p],
  Scope[ {n:=Get[distribution,1],p:=Get[distribution,2]},
   If[ !IsInteger[n]||!IsNumber[p]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[n<=0,Throw["PDF: NegativeBinomialD: First parameter is expected to be positive."]];
   If[p<=0||p>=1,Throw["PDF: NegativeBinomialD: Parameter is expected to be a non-negative probability."]];
   If[x<0,Return[0]];
   Binomial[x+n-1,n-1]*p^n*(1-p)^x
  ],

  distribution==PoissonD[_mu],
  Scope[ {mu:=Get[distribution,1]},
   If[ !IsNumber[mu]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'PDF] ]];
   If[mu<=0,Throw["PDF: PoissonD: Parameter is expected to be positive."]];
   If[x<0,Return[0]];
   Exp[-mu]*mu^x/Factorial[x]
  ],

  Put[{distribution,x},0,'PDF]

 ]

>>,
SetMask[PDF,"11"], Protect[PDF],

(*============================================================================*)

CDF[distribution,x]:=<<

 If[!IsNumber[x],Return[ Put[{distribution,x},0,'CDF] ]];

 If[

  distribution==SampleD[_data],
  Scope[ {data:=Get[distribution,1],prob,step},
   If[ !IsNumericVector[data]||Length[data]<2 ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[Rtol[1e-13,IsInvariable[data]],Return[x<Mean[data]?0.0?1.0]];
   prob=0.0; step=1.0/(Length[data]-1); data=Table[{i,prob};prob+=step,(=i,Sort[data]=)];
   Interpolate[FunctionTable[data],x]
  ],

  distribution==NormalD[_mu,_sigma],
  Scope[ {mu:=Get[distribution,1],sigma:=Get[distribution,2]},
   If[ !IsNumber[mu]||!IsNumber[sigma] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[sigma<=0,Throw["CDF: NormalD: Second parameter is expected to be positive."]];
   (Erf[(x-mu)/(Sqrt[2]*sigma)]+1)/2
  ],

  distribution==StudentTD[_n],
  Scope[ {n:=Get[distribution,1]},
   If[ !IsNumber[n] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[n<=0,Throw["CDF: StudentTD: Parameter is expected to be positive."]];
   (1+Sign[x]*BetaQ[n/2,1/2,n/(n+x^2)])/2
  ],

  distribution==ChiSquareD[_n],
  Scope[ {n:=Get[distribution,1]},
   If[ !IsNumber[n] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[n<=0,Throw["CDF: ChiSquareD: Parameter is expected to be positive."]];
   If[x<0,Return[0]];
   GammaP[n/2,x/2]
  ],

  distribution==FRatioD[_n1,_n2],
  Scope[ {n1:=Get[distribution,1],n2:=Get[distribution,2]},
   If[ !IsNumber[n1]||!IsNumber[n2] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[n1<=0||n2<=0,Throw["CDF: FRatioD: Parameters are expected to be positive."]];
   If[x<0,Return[0]];
   BetaQ[n2/2,n1/2,n2/(n2+n1*x)]
  ],

  distribution==BetaD[_a,_b],
  Scope[ {a:=Get[distribution,1],b:=Get[distribution,2]},
   If[ !IsNumber[a]||!IsNumber[b] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[a<=0||b<=0,Throw["CDF: BetaD: Parameters are expected to be positive."]];
   If[x<0,Return[0]];
   If[x>1,Return[1]];
   BetaP[a,b,x]
  ],

  distribution==CauchyD[_a,_b],
  Scope[ {a:=Get[distribution,1],b:=Get[distribution,2]},
   If[ !IsNumber[a]||!IsNumber[b] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[b<=0,Throw["CDF: CauchyD: Second parameter is expected to be positive."]];
   Atan[(x-a)/b]/Pi+1/2
  ],

  distribution==ChiD[_n],
  Scope[ {n:=Get[distribution,1]},
   If[ !IsNumber[n] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[n<=0,Throw["CDF: ChiD: Parameter is expected to be positive."]];
   If[x<0,Return[0]];
   GammaP[n/2,x^2/2]
  ],

  distribution==ExponentialD[_lambda],
  Scope[ {lambda:=Get[distribution,1]},
   If[ !IsNumber[lambda] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[lambda<=0,Throw["CDF: ExponentialD: Parameter is expected to be positive."]];
   If[x<0,Return[0]];
   1-Exp[-lambda*x]
  ],

  distribution==ExtremeValueD[_alpha,_beta],
  Scope[ {alpha:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[alpha]||!IsNumber[beta] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[beta<=0,Throw["CDF: ExtremeValueD: Second parameter is expected to be positive."]];
   Exp[-Exp[-(x-alpha)/beta]]
  ],

  distribution==GammaD[_alpha,_beta],
  Scope[ {alpha:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[alpha]||!IsNumber[beta] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[alpha<=0||beta<=0,Throw["CDF: GammaD: Parameters are expected to be positive."]];
   If[x<0,Return[0]];
   GammaP[alpha,x/beta]
  ],

  distribution==HalfNormalD[_theta],
  Scope[ {theta:=Get[distribution,1]},
   If[ !IsNumber[theta] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[theta<=0,Throw["CDF: HalfNormalD: Parameter is expected to be positive."]];
   If[x<0,Return[0]];
   Erf[theta*x/Sqrt[Pi]]
  ],

  distribution==LaplaceD[_mu,_beta],
  Scope[ {mu:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[mu]||!IsNumber[beta] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[beta<=0,Throw["CDF: LaplaceD: Second parameter is expected to be positive."]];
   ((1+Sign[x-mu]*(1-Exp[-Sign[x-mu]*(x-mu)/beta]))/2)
  ],

  distribution==LogisticD[_mu,_beta],
  Scope[ {mu:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[mu]||!IsNumber[beta] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[beta<=0,Throw["CDF: LogisticD: Second parameter is expected to be positive."]];
   1/(1+Exp[-(x-mu)/beta])
  ],

  distribution==LogNormalD[_mu,_sigma],
  Scope[ {mu:=Get[distribution,1],sigma:=Get[distribution,2]},
   If[ !IsNumber[mu]||!IsNumber[sigma] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[sigma<=0,Throw["CDF: LogNormalD: Second parameter is expected to be positive."]];
   If[x<=0,Return[0]];
   (Erf[(-mu + Ln[x])/(Sqrt[2]*sigma)]+1)/2
  ],

  distribution==ParetoD[_k,_alpha],
  Scope[ {k:=Get[distribution,1],alpha:=Get[distribution,2]},
   If[ !IsNumber[k]||!IsNumber[alpha] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[k<=0||alpha<=0,Throw["CDF: ParetoD: Parameters are expected to be positive."]];
   If[x<0,Return[0]];
   1-(k/x)^alpha
  ],

  distribution==RayleighD[_sigma],
  Scope[ {sigma:=Get[distribution,1]},
   If[ !IsNumber[sigma] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[sigma<=0,Throw["CDF: RayleighD: Parameter is expected to be positive."]];
   If[x<0,Return[0]];
   1-Exp[-x^2/(2*sigma^2)]
  ],

  distribution==UniformD[_min,_max],
  Scope[ {min:=Get[distribution,1],max:=Get[distribution,2]},
   If[ !IsNumber[min]||!IsNumber[max] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[max-min<=0,Throw["CDF: UniformD: Second parameter is expected to be greater than the first."]];
   (1+Sign[x-min]*(x-min)/(max-min)-(max-x)*Sign[max-x]/(max-min))/2
  ],

  distribution==WeibullD[_alpha,_beta],
  Scope[ {alpha:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[alpha]||!IsNumber[beta] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[alpha<=0||beta<=0,Throw["CDF: WeibullD: Parameters are expected to be positive."]];
   If[x<0,Return[0]];
   1-Exp[-(x/beta)^alpha]
  ],

  distribution==BernoulliD[_p],
  Scope[ {p:=Get[distribution,1]},
   If[ !IsNumber[p]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[p<0||p>1,Throw["CDF: BernoulliD: Parameter is expected to be a probability."]];
   If[x<0,Return[0]];
   If[x>=1,Return[1]];
   1-p
  ],

  distribution==BinomialD[_n,_p],
  Scope[ {n:=Get[distribution,1],p:=Get[distribution,2]},
   If[ !IsInteger[n]||!IsNumber[p]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[n<1,Throw["CDF: BinomialD: First parameter is expected to be a positive integer."]];
   If[p<0||p>1,Throw["CDF: BinomialD: Second parameter is expected to be a probability."]];
   If[x<0,Return[0]];
   If[x>=n,Return[1]];
   If[x==0,Return[(1-p)^n]];
   BetaP[n-x,x+1,1-p]
  ],

  distribution==DiscreteUniformD[_n],
  Scope[ {n:=Get[distribution,1]},
   If[ !IsInteger[n]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[n<1,Throw["CDF: DiscreteUniformD: Parameter is expected to be a positive integer."]];
   If[x<1,Return[0]];
   If[x>=n,Return[1]];
   x/n
  ],

  distribution==GeometricD[_p],
  Scope[ {p:=Get[distribution,1]},
   If[ !IsNumber[p]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[p<0||p>1,Throw["CDF: GeometricD: Parameter is expected to be a probability."]];
   If[x<0,Return[0]];
   1-(1-p)^(x+1)
  ],

  distribution==HypergeometricD[_n,_nsucc,_ntot],
  Scope[ {n:=Get[distribution,1],nsucc:=Get[distribution,2],ntot:=Get[distribution,3]},
   If[ !IsInteger[n]||!IsInteger[nsucc]||!IsInteger[ntot]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[n<0,Throw["CDF: HypergeometricD: First parameter is expected to be non-negative."]];
   If[nsucc<0,Throw["CDF: HypergeometricD: Second parameter is expected to be non-negative."]];
   If[ntot<1,Throw["CDF: HypergeometricD: Third parameter is expected to be positive."]];
   If[n>ntot,Throw["CDF: HypergeometricD: First parameter must be less or equal than the third."]];
   If[nsucc>ntot,Throw["CDF: HypergeometricD: Second parameter must be less or equal than the third."]];
   If[x<Max[{0,n-ntot+nsucc}],Return[0]];
   If[x>=Min[{n,nsucc}],Return[1]];
   Sum[
    (Binomial[nsucc,k]*Binomial[ntot-nsucc,n-k])/Binomial[ntot,n],
    (=k,Max[{0,n-ntot+nsucc}]:Min[{x,n,nsucc}]=)
   ]
  ],

  distribution==LogSeriesD[_theta],
  Scope[ {theta:=Get[distribution,1]},
   If[ !IsNumber[theta]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[theta<=0||theta>=1,Throw["CDF: LogSeriesD: Parameter is expected to be less than 1 and non-negative."]];
   If[x<=0,Return[0]];
   Sum[theta^k/(-k*Ln[1-theta]),(=k,1:x=)]
  ],

  distribution==NegativeBinomialD[_n,_p],
  Scope[ {n:=Get[distribution,1],p:=Get[distribution,2]},
   If[ !IsInteger[n]||!IsNumber[p]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[n<=0,Throw["CDF: NegativeBinomialD: First parameter is expected to be positive."]];
   If[p<=0||p>=1,Throw["CDF: NegativeBinomialD: Parameter is expected to be a non-negative probability."]];
   If[x<0,Return[0]];
   BetaP[n,x+1,p]
  ],

  distribution==PoissonD[_mu],
  Scope[ {mu:=Get[distribution,1]},
   If[ !IsNumber[mu]||!IsInteger[x] ,Return[ Put[{distribution,x},0,'CDF] ]];
   If[mu<=0,Throw["CDF: PoissonD: Parameter is expected to be positive."]];
   If[x<0,Return[0]];
   GammaQ[x+1,mu]
  ],

  Put[{distribution,x},0,'CDF]

 ]

>>,
SetMask[CDF,"11"], Protect[CDF],

(*============================================================================*)

Quantile[distribution,p]:=<<

 If[!IsNumber[p],Return[ Put[{distribution,p},0,'Quantile] ]];

 If[p<0||p>1,Throw["Quantile: Second argument is expected to be a probability."]];

 If[

  distribution==SampleD[_data],
  Scope[ {data:=Get[distribution,1],prob,step},
   If[ !IsNumericVector[data]||Length[data]<2 ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[Rtol[1e-13,IsInvariable[data]],Return[AsReal[Mean[data]]]];
   prob=0.0; step=1.0/(Length[data]-1); data=Table[{i,prob};prob+=step,(=i,Sort[data]=)];
   Interpolate[InverseTable[data],p]
  ],

  distribution==NormalD[_mu,_sigma],
  Scope[ {mu:=Get[distribution,1],sigma:=Get[distribution,2]},
   If[ !IsNumber[mu]||!IsNumber[sigma] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[sigma<=0,Throw["Quantile: NormalD: Second parameter is expected to be positive."]];
   mu+sigma*Iphi[p]
  ],

  distribution==StudentTD[_n],
  Scope[ {n:=Get[distribution,1]},
   If[ !IsNumber[n] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[n<=0,Throw["Quantile: StudentTD: Parameter is expected to be positive."]];
   Sign[2*p-1]*Sqrt[n/BetaX[n/2,1/2,1-(2*p-1)*Sign[2*p-1]]-n]
  ],

  distribution==ChiSquareD[_n],
  Scope[ {n:=Get[distribution,1]},
   If[ !IsNumber[n] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[n<=0,Throw["Quantile: ChiSquareD: Parameter is expected to be positive."]];
   2*GammaY[n/2,1-p]
  ],

  distribution==FRatioD[_n1,_n2],
  Scope[ {n1:=Get[distribution,1],n2:=Get[distribution,2]},
   If[ !IsNumber[n1]||!IsNumber[n2] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[n1<=0||n2<=0,Throw["Quantile: FRatioD: Parameters are expected to be positive."]];
   (n2/BetaX[n2/2,n1/2,1-p]-n2)/n1
  ],

  distribution==BetaD[_a,_b],
  Scope[ {a:=Get[distribution,1],b:=Get[distribution,2]},
   If[ !IsNumber[a]||!IsNumber[b] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[a<=0||b<=0,Throw["Quantile: BetaD: Parameters are expected to be positive."]];
   BetaX[a,b,p]
  ],

  distribution==CauchyD[_a,_b],
  Scope[ {a:=Get[distribution,1],b:=Get[distribution,2]},
   If[ !IsNumber[a]||!IsNumber[b] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[b<=0,Throw["Quantile: CauchyD: Second parameter is expected to be positive."]];
   a+b*Tan[Pi*(p-1/2)]
  ],

  distribution==ChiD[_n],
  Scope[ {n:=Get[distribution,1]},
   If[ !IsNumber[n] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[n<=0,Throw["Quantile: ChiD: Parameter is expected to be positive."]];
   Sqrt[2]*Sqrt[GammaY[n/2,1-p]]
  ],

  distribution==ExponentialD[_lambda],
  Scope[ {lambda:=Get[distribution,1]},
   If[ !IsNumber[lambda] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[lambda<=0,Throw["Quantile: ExponentialD: Parameter is expected to be positive."]];
   -Ln[1-p]/lambda
  ],

  distribution==ExtremeValueD[_alpha,_beta],
  Scope[ {alpha:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[alpha]||!IsNumber[beta] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[beta<=0,Throw["Quantile: ExtremeValueD: Second parameter is expected to be positive."]];
   alpha-beta*Ln[Ln[1/p]]
  ],

  distribution==GammaD[_alpha,_beta],
  Scope[ {alpha:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[alpha]||!IsNumber[beta] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[alpha<=0||beta<=0,Throw["Quantile: GammaD: Parameters are expected to be positive."]];
   beta*GammaY[alpha,1-p]
  ],

  distribution==HalfNormalD[_theta],
  Scope[ {theta:=Get[distribution,1]},
   If[ !IsNumber[theta] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[theta<=0,Throw["Quantile: HalfNormalD: Parameter is expected to be positive."]];
   Sqrt[Pi]*Ierf[p]/theta
  ],

  distribution==LaplaceD[_mu,_beta],
  Scope[ {mu:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[mu]||!IsNumber[beta] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[beta<=0,Throw["Quantile: LaplaceD: Second parameter is expected to be positive."]];
   mu-beta*Sign[2*p-1]*Ln[1-Sign[2*p-1]*(2*p-1)]
  ],

  distribution==LogisticD[_mu,_beta],
  Scope[ {mu:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[mu]||!IsNumber[beta] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[beta<=0,Throw["Quantile: LogisticD: Second parameter is expected to be positive."]];
   mu-beta*Sign[beta]*Ln[1-p]+beta*Sign[beta]*Ln[p]
  ],

  distribution==LogNormalD[_mu,_sigma],
  Scope[ {mu:=Get[distribution,1],sigma:=Get[distribution,2]},
   If[ !IsNumber[mu]||!IsNumber[sigma] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[sigma<=0,Throw["Quantile: LogNormalD: Second parameter is expected to be positive."]];
   Exp[Ierf[2*p-1]*Sqrt[2]*sigma+mu]
  ],

  distribution==ParetoD[_k,_alpha],
  Scope[ {k:=Get[distribution,1],alpha:=Get[distribution,2]},
   If[ !IsNumber[k]||!IsNumber[alpha] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[k<=0||alpha<=0,Throw["Quantile: ParetoD: Parameters are expected to be positive."]];
   k/(1-p)^(1/alpha)
  ],

  distribution==RayleighD[_sigma],
  Scope[ {sigma:=Get[distribution,1]},
   If[ !IsNumber[sigma] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[sigma<=0,Throw["Quantile: RayleighD: Parameter is expected to be positive."]];
   sigma*Sqrt[Ln[1/(1-p)^2]]
  ],

  distribution==UniformD[_min,_max],
  Scope[ {min:=Get[distribution,1],max:=Get[distribution,2]},
   If[ !IsNumber[min]||!IsNumber[max] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[max-min<=0,Throw["Quantile: UniformD: Second parameter is expected to be greater than the first."]];
   p*max+(1-p)*min
  ],

  distribution==WeibullD[_alpha,_beta],
  Scope[ {alpha:=Get[distribution,1],beta:=Get[distribution,2]},
   If[ !IsNumber[alpha]||!IsNumber[beta] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[alpha<=0||beta<=0,Throw["Quantile: WeibullD: Parameters are expected to be positive."]];
   beta*(-Ln[1-p])^(1/alpha)
  ],

  distribution==BernoulliD[_pp],
  Scope[ {pp:=Get[distribution,1]},
   If[ !IsNumber[pp] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[pp<0||pp>1,Throw["Quantile: BernoulliD: Parameter is expected to be a probability."]];
   If[p<=1-pp,0,1]
  ],

  distribution==BinomialD[_n,_pp],
  Scope[ {n:=Get[distribution,1],pp:=Get[distribution,2]},
   If[ !IsInteger[n]||!IsNumber[pp] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[n<1,Throw["Quantile: BinomialD: First parameter is expected to be a positive integer."]];
   If[pp<0||pp>1,Throw["Quantile: BinomialD: Second parameter is expected to be a probability."]];
   Scope[{low:=0,high:=n,mid},
    If[p==0,Return[low]];
    If[p==1,Return[high]];
    If[p<=PDF[BinomialD[n,pp],low],Return[low]];
    While[high-low>1,
     mid=Floor[(high+low)/2];
     If[CDF[BinomialD[n,pp],mid]<p,low=mid,high=mid]
    ];
    high
   ]
  ],

  distribution==DiscreteUniformD[_n],
  Scope[ {n:=Get[distribution,1]},
   If[ !IsInteger[n] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[n<1,Throw["Quantile: DiscreteUniformD: Parameter is expected to be a positive integer."]];
   Max[{Ceiling[n*p],1}]
  ],

  distribution==GeometricD[_pp],
  Scope[ {pp:=Get[distribution,1]},
   If[ !IsNumber[pp] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[pp<0||pp>1,Throw["Quantile: GeometricD: Parameter is expected to be a probability."]];
   If[pp<p,Ceiling[Ln[1-p]/Ln[1-pp]]-1,0]
  ],

  distribution==HypergeometricD[_n,_nsucc,_ntot],
  Scope[ {n:=Get[distribution,1],nsucc:=Get[distribution,2],ntot:=Get[distribution,3]},
   If[ !IsInteger[n]||!IsInteger[nsucc]||!IsInteger[ntot] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[n<0,Throw["Quantile: HypergeometricD: First parameter is expected to be non-negative."]];
   If[nsucc<0,Throw["Quantile: HypergeometricD: Second parameter is expected to be non-negative."]];
   If[ntot<1,Throw["Quantile: HypergeometricD: Third parameter is expected to be positive."]];
   If[n>ntot,Throw["Quantile: HypergeometricD: First parameter must be less or equal than the third."]];
   If[nsucc>ntot,Throw["Quantile: HypergeometricD: Second parameter must be less or equal than the third."]];
   Scope[{low:=Max[{0,n-ntot+nsucc}],high:=Min[{x,n,nsucc}],mid},
    If[p==0,Return[low]];
    If[p==1,Return[high]];
    If[p<=PDF[HypergeometricD[n,nsucc,ntot],low],Return[low]];
    While[high-low>1,
     mid=Floor[(high+low)/2];
     If[CDF[HypergeometricD[n,nsucc,ntot],mid]<p,low=mid,high=mid]
    ];
    high
   ]
  ],

  distribution==LogSeriesD[_theta],
  Scope[ {theta:=Get[distribution,1]},
   If[ !IsNumber[theta] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[theta<=0||theta>=1,Throw["Quantile: LogSeriesD: Parameter is expected to be less than 1 and non-negative."]];
   Scope[{low:=1,high:=2,mid},
    If[p<=PDF[LogSeriesD[theta],low],Return[low]];
    While[CDF[LogSeriesD[theta],high]<p,high*=2];
    low=high/2;
    While[high-low>1,
     mid=Floor[(high+low)/2];
     If[CDF[LogSeriesD[theta],mid]<p,low=mid,high=mid]
    ];
    high
   ]
  ],

  distribution==NegativeBinomialD[_n,_pp],
  Scope[ {n:=Get[distribution,1],pp:=Get[distribution,2]},
   If[ !IsInteger[n]||!IsNumber[pp] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[n<=0,Throw["Quantile: NegativeBinomialD: First parameter is expected to be positive."]];
   If[pp<=0||pp>=1,Throw["Quantile: NegativeBinomialD: Parameter is expected to be a non-negative probability."]];
   Scope[{low:=0,high:=1,mid},
    If[p<=PDF[NegativeBinomialD[n,pp],low],Return[low]];
    While[CDF[NegativeBinomialD[n,pp],high]<p,high*=2];
    low=high/2;
    While[high-low>1,
     mid=Floor[(high+low)/2];
     If[CDF[NegativeBinomialD[n,pp],mid]<p,low=mid,high=mid]
    ];
    high
   ]
  ],

  distribution==PoissonD[_mu],
  Scope[ {mu:=Get[distribution,1]},
   If[ !IsNumber[mu] ,Return[ Put[{distribution,p},0,'Quantile] ]];
   If[mu<=0,Throw["Quantile: PoissonD: Parameter is expected to be positive."]];
   Scope[{low:=0,high:=1,mid},
    If[p<=PDF[PoissonD[mu],low],Return[low]];
    While[CDF[PoissonD[mu],high]<p,high*=2];
    low=high/2;
    While[high-low>1,
     mid=Floor[(high+low)/2];
     If[CDF[PoissonD[mu],mid]<p,low=mid,high=mid]
    ];
    high
   ]
  ],

  Put[{distribution,p},0,'Quantile]

 ]

>>,
SetMask[Quantile,"11"], Protect[Quantile],

(*============================================================================*)

Random[distribution]:=<<

 If[

  distribution==SampleD[_data], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==NormalD[_mu,_sigma], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==StudentTD[_n], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==ChiSquareD[_n], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==FRatioD[_n1,_n2], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==BetaD[_a,_b], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==CauchyD[_a,_b], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==ChiD[_n], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==ExponentialD[_lambda], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==ExtremeValueD[_alpha,_beta], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==GammaD[_alpha,_beta], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==HalfNormalD[_theta], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==LaplaceD[_mu,_beta], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==LogisticD[_mu,_beta], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==LogNormalD[_mu,_sigma], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==ParetoD[_k,_alpha], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==RayleighD[_sigma], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==UniformD[_min,_max], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==WeibullD[_alpha,_beta], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==BernoulliD[_p], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==BinomialD[_n,_p], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==DiscreteUniformD[_n], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==GeometricD[_p], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==HypergeometricD[_n,_nsucc,_ntot], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==LogSeriesD[_theta], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==NegativeBinomialD[_n,_p], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  distribution==PoissonD[_mu], Quantile[distribution,RandomReal[1e-7,1-1e-7]],

  Put[{distribution},0,'Random]

 ]

>>,
SetMask[Random,"1"], Protect[Random],

(*============================================================================*)

()
