Probability & Statistics

A Preview

CellDisplay@Framed[Grid[Partition[#,3]]&@ParallelTable[Plot3D[Evaluate@PDF[BinormalDistribution[ρ],{x,y}],{x,-3,3},{y,-3,3},PlotRange->All,MeshFunctions->{#3&},MeshShading->{None,Green,None,Yellow},PlotPoints->35,PlotLabel->ρ,ImageSize->Small,Ticks->None,PlotRange->{Full,Full,{0,0.5}}],{ρ,{-0.95,-0.75,-0.5,-0.25,0,0.25,0.5,0.75,0.95}}],RoundingRadius->10,Background->Lighter[Gray,0.9]]

Probability & Statistics

High Level Functions

Probability, NProbability

Expectation, NExpectation

Distribution Properties

PDF, CDF, SurvivalFunction (SF), HazardFunction (HF), Quantile (QF)

Moment / MGF, CentralMoment / CMGF, FactorialMoment / FMGF, Cumulant / CGF, CF

RandomNumber, DistributionEstimate

Mean, Variance, StandardDeviation, Skewness, Kurtosis, Covariance, Correlation, ...

Median, Quartiles, QuartileDeviation, QuartileSkewness, ...

Likelihood, LogLikelihood, ...

Parametric Distributions

Distributions constructed from parameters.

Univariate Discrete

Univariate Continuous

Multivariate Discrete

Multivariate Continuous

Meta Distributions

Distributions constructed from other distributions.

TransformDistribution

OrderDistribution

ComponentMixtureDistribution

CensoredDistribution

TruncatedDistribution

MarginalDistribution

IndependentDistribution

CopulaDistribution

User Distributions

Distributions constructed from formulas.

ProbabilityDistribution

Data Distributions

Distributions constructed “directly” from data.

EmpiricalDistribution

HistogramDistribution

SmoothKernelDistribution

Distribution Visualization

Understand and compare distributions.

Plot, Plot3D, DiscretePlot, DiscretePlot3D, ...

Histogram, Histogram3D, DensityHistogram, PairedHistogram, ...

QuantilePlot, ProbabilityPlot, ...

BoxWhiskerChart, DistributionChart, ...

Probability

Probability[event,…] in terms of PDF f:

FormattedDefinitionGrid[{{HoldForm@((∑_(ξ=-∞)^∞)Boole[event]f(ξ) ),HoldForm@((∑_((Subscript[ξ, 1])=-∞)^∞)(∑_((Subscript[ξ, n])=-∞)^∞)Boole[event] f((Subscript[ξ, 1]),…,(Subscript[ξ, n])))},{HoldForm@((∫_(-∞)^∞)Boole[event]f(ξ) \[DifferentialD]ξ),HoldForm@((∫_(-∞)^∞)(∫_(-∞)^∞)Boole[event] f((Subscript[ξ, 1]),…,(Subscript[ξ, n]))\[DifferentialD](Subscript[ξ, n])\[DifferentialD](Subscript[ξ, 1]))}}]

Univariate Discrete

Discrete univariate probabilities:

Probability[X≥1,X∈BinomialDistribution[n,p]]

Nonlinear conditions:

Probability[(((X-3))^2)≥X,X∈ZipfDistribution[ρ]]

Logical combinations:

Probability[k^2>10||k<2,k∈PoissonDistribution[λ]]

An airplane engine fails with probability p. Compute the probability that no more then 2 out of 4 fail.

prob1=Probability[X≤2,X∈BinomialDistribution[4,p]]

Compute the probability that no more then 1 out of 2 fail.

prob2=Probability[X≤1,X∈BinomialDistribution[2,p]]

When are four engines better then two engines?

Reduce[FunctionExpand@prob1>prob2&&0<p<1,p]

Plot[{Tooltip[prob1,"4 engines"],Tooltip[prob2,"2 engines"]},{p,0,1}]

Univariate Continuous

Find the probability between k standard deviations:

Table[Probability[μ-k σ<X<μ+k σ,X∈ NormalDistribution[μ,σ]],{k,5}]

N[%]

Typically can compute probabilities for polynomial inequalities:

Probability[(X^3)+2X+1<3,X∈ExponentialDistribution[λ]]

And to some extent for transcendental inequalities:

prob=Probability[Exp[-X]>2(X^2)+1,X∈NormalDistribution[μ,σ]]

N[prob]

Plot3D[prob,{μ,-3,3},{σ,1/2,2}]

A company manufactures nails with length normally distributed and a mean of 0.5 inches. If  50% of the produced nails have lengths between 0.495 and 0.505, find the standard deviation:

Probability[0.495<X<0.505,X∈NormalDistribution[0.5,σ]]

Solve[%==0.5,σ]

A battery has a lifespan that is exponentially distributed with parameter 3000. Find the probability that a random battery has a lifespan of less than 2500 hours:

Probability[X<2500,X∈ ExponentialDistribution[1/3000]]

N[%]

Multivariate

Continuous multivariate probabilities:

Probability[X<2&&Y<1,{X,Y}∈BinormalDistribution[1/2]]

NProbability[X<2&&Y<1,{X,Y}∈BinormalDistribution[1/2]]

NProbability[X≤1/2&&Y≤1,{X,Y}∈DirichletDistribution[{1,2,3}]]

Conditional Probability

Compute conditional probabilities:

Probability[ConditionalEvent[x^2>12,x>2],x∈PoissonDistribution[λ]]

Probability[ConditionalEvent[x^2>12,x>2],x∈ExponentialDistribution[λ]]

Multivariate:

Probability[ConditionalEvent[x^2>12,y<3],{x,y}∈BinormalDistribution[1/2]]

N[%]

Probability[ConditionalEvent[x^2>12,y∈Reals],{x,y}∈BinormalDistribution[ρ]]

Probability[x^2>12,x∈NormalDistribution[]]

Expectation

Expectation[g,…] in terms of PDF f:

FormattedDefinitionGrid[{{HoldForm@((∑_(ξ=-∞)^∞)g[ξ]f(ξ) ),HoldForm@((∑_((Subscript[ξ, 1])=-∞)^∞)(∑_((Subscript[ξ, n])=-∞)^∞)g[(Subscript[ξ, 1]),…,(Subscript[ξ, n])] f((Subscript[ξ, 1]),…,(Subscript[ξ, n])))},{HoldForm@((∫_(-∞)^∞)g[ξ]f(ξ) \[DifferentialD]ξ),HoldForm@((∫_(-∞)^∞)(∫_(-∞)^∞)g[(Subscript[ξ, 1]),…,(Subscript[ξ, n])] f((Subscript[ξ, 1]),…,(Subscript[ξ, n]))\[DifferentialD](Subscript[ξ, n])\[DifferentialD](Subscript[ξ, 1]))}}]

Basic Examples

Compute moments:

Expectation[k^2,k∈PoissonDistribution[μ]]

Even nonstandard ones:

Expectation[(2^(-k)),k∈PoissonDistribution[μ]]

A basketball player shoots free throws until he hits 4 of them. His probability of scoring in any one of them is 0.7. Find the number of shots the player is expected to shoot:

Expectation[X+4,X∈NegativeBinomialDistribution[4,0.7]]

Assume the probability of fouling for each minute interval is 0.1 independently. A basketball player fouls out after 6 fouls. Find expected playing time until foul out:

Expectation[X+6,X∈NegativeBinomialDistribution[6,0.1]]

Here the normal moments don’t exist:

Expectation[Exp[-Abs[X]],X∈CauchyDistribution[1,2]]

N[%]

NExpectation[Exp[-Abs[X]],X∈CauchyDistribution[1,2]]

Distribution Properties – Distribution Functions

Several different functions and their interpretation as probabilities:

CellDisplay@Style[TabView[{"PDF"->DefinitionGrid[{{HoldForm@Probability[X==k,…],HoldForm@Probability[(Subscript[X, 1])==(Subscript[k, 1])∧⋯∧(Subscript[X, n])==(Subscript[k, n]),…]},{HoldForm@Limit[((Probability[x-ϵ<X≤x,…])/ϵ),ϵ->0],HoldForm@Limit[((Probability[(Subscript[x, 1])-(Subscript[ϵ, 1])<(Subscript[X, 1])≤(Subscript[x, 1])∧⋯∧(Subscript[x, n])-(Subscript[ϵ, n])<(Subscript[X, n])≤(Subscript[x, n]),…])/((Subscript[ϵ, 1])⋯ (Subscript[ϵ, n]))),(Subscript[ϵ, 1])->0,…,(Subscript[ϵ, n])->0]}}],
"CDF"->DefinitionGrid[{{HoldForm@Probability[X≤k,…],HoldForm@Probability[(Subscript[X, 1])≤(Subscript[k, 1])∧⋯∧(Subscript[X, n])≤(Subscript[k, n]),…]},{HoldForm@Probability[X≤x,…],HoldForm@Probability[(Subscript[X, 1])≤(Subscript[x, 1])∧⋯∧(Subscript[X, n])≤(Subscript[x, n]),…]}}],
"SF"->DefinitionGrid[{{HoldForm@Probability[X>k,…],HoldForm@Probability[(Subscript[X, 1])>(Subscript[k, 1])∧⋯∧(Subscript[X, n])>(Subscript[k, n]),…]},{HoldForm@Probability[X>x,…],HoldForm@Probability[(Subscript[X, 1])>(Subscript[x, 1])∧⋯∧(Subscript[X, n])>(Subscript[x, n]),…]}}],
"HF"->DefinitionGrid[{{HoldForm@Probability[X==k|X≥k,…],HoldForm@Probability[(Subscript[X, 1])==(Subscript[k, 1])∧⋯∧(Subscript[X, n])==(Subscript[k, n])|(Subscript[X, 1])≥(Subscript[k, 1])∧⋯∧(Subscript[X, n])≥(Subscript[k, n]),…]},{HoldForm@Limit[((Probability[x<X≤x+ϵ|X≥x,…])/ϵ),ϵ->0],HoldForm@Limit[(1/((Subscript[ϵ, 1])⋯ (Subscript[ϵ, n])))Probability[(Subscript[x, 1])<(Subscript[X, 1])≤(Subscript[x, 1])+(Subscript[ϵ, 1])∧⋯∧(Subscript[x, n])<(Subscript[X, n])≤(Subscript[x, n])+(Subscript[ϵ, n])|(Subscript[X, 1])≥(Subscript[x, 1])∧⋯∧(Subscript[X, n])≥(Subscript[x, n]),…],(Subscript[ϵ, 1])->0,…,(Subscript[ϵ, n])->0]}}]}],"Text"]

DistributionPlot[NormalDistribution[],{-3,3}]//CellDisplay

Basic Examples

Continuous univariate:

DistributionPlot[GompertzMakehamDistribution[1,2],{0,2}]

Discrete univariate:

DistributionPlot[BinomialDistribution[10,2/3],{0,11}]

Continuous bivariate:

DistributionPlot3D[BinormalDistribution[1/2],{-3,3},{-3,3}]

Discrete bivariate:

DistributionPlot3D[MultinomialDistribution[10,{1/3,2/3}],{0,10},{0,10}]

These can be accessed individually for each distribution:

PDF[ExponentialDistribution[λ],x]

CDF[ExponentialDistribution[λ],x]

SurvivalFunction[ExponentialDistribution[λ],x]

HazardFunction[ExponentialDistribution[λ],x]

Examples

For discrete distributions the PDF is a sequence:

PDF[BinomialDistribution[n,p],k]

Multivariate distributions take a vector of arguments:

PDF[BinormalDistribution[ρ],{x,y}]

Manipulate[Plot3D[PDF[BinormalDistribution[ρ],{x,y}],{x,-2,2},{y,-2,2},MeshFunctions->{#3&},PlotPoints->35],{ρ,-0.9,0.9},Deployed->True]

Multivariate distributions may have matrix and vector parameters:

Σ={{2,1},{1,2}};

PDF[MultivariateTDistribution[Σ,ν],{x,y}]

Can easily prove relations between distributions:

PDF[StudentTDistribution[ν],x]

Limit[PDF[StudentTDistribution[ν],x],ν->∞]

PDF[NormalDistribution[0,1],x]

If no closed form formula exist you can still evaluate numerically:

PDF[StableDistribution[0,1,1,0,1],x]

PDF[StableDistribution[0,1,1,0,1],0.1]

Plot[PDF[StableDistribution[0,1,1,0,1],x],{x,-5,5}]

Distribution Properties – Moments

Moments and generating functions ... starting from moments:

CellDisplay@TabView[{"Moment"->MomentDefinitionGrid[{{HoldForm[(Subscript[μ, n])==Expectation[(X^n),…]],HoldForm[MGF[t]==GeneratingFunction[((Subscript[μ, n])/(n!)),n,t]]},{HoldForm[(Subscript[μ, (Subscript[n, 1]),…,(Subscript[n, k])])==Expectation[((Subscript[X, 1])^(Subscript[n, 1]))⋯ ((Subscript[X, k])^(Subscript[n, k])),…]],HoldForm[MGF[(Subscript[t, 1]),…,(Subscript[t, k])]==GeneratingFunction[((Subscript[μ, (Subscript[n, 1]),…,(Subscript[n, k])])/((Subscript[n, 1])!⋯ (Subscript[n, k])!)),{(Subscript[n, 1]),…,(Subscript[n, k])},{(Subscript[t, 1]),…,(Subscript[t, k])}]]}}],"CentralMoment"->MomentDefinitionGrid[{{HoldForm[(Subscript[(OverTilde[μ]), n])==Expectation[(((X-μ))^n),X∈𝒟]],HoldForm[CMGF[t]==GeneratingFunction[((Subscript[(OverTilde[μ]), n])/(n!)),n,t]]},{HoldForm[(Subscript[(OverTilde[μ]), (Subscript[n, 1]),…,(Subscript[n, k])])==Expectation[((((Subscript[X, 1])-(Subscript[μ, (Subscript[X, 1])])))^(Subscript[n, 1]))⋯ ((((Subscript[X, k])-(Subscript[μ, (Subscript[X, k])])))^(Subscript[n, k])),{(Subscript[X, 1]),…,(Subscript[X, k])}∈𝒟]],HoldForm[CMGF[(Subscript[t, 1]),…,(Subscript[t, k])]==GeneratingFunction[((Subscript[(OverTilde[μ]), (Subscript[n, 1]),…,(Subscript[n, k])])/((Subscript[n, 1])!⋯ (Subscript[n, k])!)),{(Subscript[n, 1]),…,(Subscript[n, k])},{(Subscript[t, 1]),…,(Subscript[t, k])}]]}}],"FactorialMoment"->MomentDefinitionGrid[{{HoldForm[(Subscript[(OverBar[μ]), n])==Expectation[FactorialPower[X,n],X∈𝒟]],HoldForm[FMGF[t]==GeneratingFunction[((Subscript[(OverBar[μ]), n])/(n!)),n,t-1]]},{HoldForm[(Subscript[(OverBar[μ]), (Subscript[n, 1]),…,(Subscript[n, k])])==Expectation[FactorialPower[(Subscript[X, 1]),(Subscript[n, 1])]⋯ FactorialPower[(Subscript[X, k]),(Subscript[n, k])],{(Subscript[X, 1]),…,(Subscript[X, k])}∈𝒟]],HoldForm[FMGF[(Subscript[t, 1]),…,(Subscript[t, k])]==GeneratingFunction[((Subscript[(OverBar[μ]), (Subscript[n, 1]),…,(Subscript[n, k])])/((Subscript[n, 1])!⋯ (Subscript[n, k])!)),{(Subscript[n, 1]),…,(Subscript[n, k])},{(Subscript[t, 1])-1,…,(Subscript[t, k])-1}]]}}]}]

Generating functions and moments ... starting from generating functions:

CellDisplay@Style[TabView[{"MGF"->MomentDefinitionGrid[{{HoldForm[(Subscript[μ, n])==Derivative[n][MGF][0]],HoldForm[MGF[t]==Expectation[Exp[t x],x∈𝒟]]},{HoldForm[(Subscript[μ, (Subscript[n, 1]),…,(Subscript[n, k])])==Derivative[(Subscript[n, 1]),…,(Subscript[n, k])][MGF][0,…,0]],HoldForm[MGF[(Subscript[t, 1]),…,(Subscript[t, k])]==Expectation[Exp[(Subscript[t, 1])(Subscript[x, 1])+⋯+(Subscript[t, k])(Subscript[x, k])],x∈𝒟]]}}],"CMGF"->MomentDefinitionGrid[{{HoldForm[(Subscript[(OverTilde[μ]), n])==Derivative[n][CMGF][0]],HoldForm[CMGF[t]==Expectation[Exp[t (x-μ)],x∈𝒟]]},{HoldForm[(Subscript[(OverTilde[μ]), (Subscript[n, 1]),…,(Subscript[n, k])])==Derivative[(Subscript[n, 1]),…,(Subscript[n, k])][CMGF][0,…,0]],HoldForm[CMGF[(Subscript[t, 1]),…,(Subscript[t, k])]==Expectation[Exp[(Subscript[t, 1])((Subscript[x, 1])-(Subscript[μ, (Subscript[x, 1])]))+⋯+(Subscript[t, k])((Subscript[x, k])-(Subscript[μ, (Subscript[x, k])]))],x∈𝒟]]}}],"FMGF"->MomentDefinitionGrid[{{HoldForm[(Subscript[(OverBar[μ]), n])==Derivative[n][FMGF][1]],HoldForm[FMGF[t]==Expectation[(t^x),x∈𝒟]]},{HoldForm[(Subscript[(OverBar[μ]), (Subscript[n, 1]),…,(Subscript[n, k])])==Derivative[(Subscript[n, 1]),…,(Subscript[n, k])][FMGF][1,…,1]],HoldForm[FMGF[(Subscript[t, 1]),…,(Subscript[t, k])]==Expectation[((Subscript[t, 1])^(Subscript[x, 1]))⋯ ((Subscript[t, k])^(Subscript[x, k])),x∈𝒟]]}}],
"CGF"->MomentDefinitionGrid[{{HoldForm[(Subscript[κ, n])==Derivative[n][CGF][0]],HoldForm[CGF[t]==Log[MGF[t]]]},{HoldForm[(Subscript[κ, (Subscript[n, 1]),…,(Subscript[n, k])])==Derivative[(Subscript[n, 1]),…,(Subscript[n, k])][CGF][0,…,0]],HoldForm[CGF[(Subscript[t, 1]),…,(Subscript[t, k])]==Log[MGF[(Subscript[t, 1]),…,(Subscript[t, k])]]]}}],
"CF"->MomentDefinitionGrid[{{…,HoldForm[CF[t]==Expectation[Exp[I t x],x∈𝒟]]},{…,HoldForm[CF[(Subscript[t, 1]),…,(Subscript[t, k])]==Expectation[Exp[I ((Subscript[t, 1])(Subscript[x, 1])+⋯+(Subscript[t, k])(Subscript[x, k]))],x∈𝒟]]}}]}],"Text"]

Basic Examples

Univariate moments:

Table[Moment[NormalDistribution[μ,σ],k],{k,4}]

Table[CentralMoment[NormalDistribution[μ,σ],k],{k,4}]

Table[FactorialMoment[NormalDistribution[μ,σ],k],{k,4}]

Table[Cumulant[NormalDistribution[μ,σ],k],{k,4}]

All at once:

MomentTabView[NormalDistribution[μ,σ],5]

Multivariate moments:

Table[Moment[BinormalDistribution[ρ],{m,n}],{m,3},{n,3}]//MatrixForm

Univariate generating functions:

MomentGeneratingFunction[ExponentialDistribution[λ],t]

CentralMomentGeneratingFunction[ExponentialDistribution[λ],t]

FactorialMomentGeneratingFunction[ExponentialDistribution[λ],t]

CumulantGeneratingFunction[ExponentialDistribution[λ],t]

CharacteristicFunction[ExponentialDistribution[λ],t]

MGFTabView[ExponentialDistribution[λ]]

Examples

Build a method of moments estimator:

{(Subscript[m, 1])==Moment[NormalDistribution[μ,σ],1],(Subscript[m, 2])==Simplify@Moment[NormalDistribution[μ,σ],2]}

data=RandomNumber[NormalDistribution[1,2],(10^4)];

NSolve[{Moment[data,1]==Moment[NormalDistribution[μ,σ],1],Moment[data,2]==Moment[NormalDistribution[μ,σ],2],σ>0},{μ,σ}]

Guess the closed form for moments:

Table[Moment[ExponentialDistribution[λ],k],{k,5}]

FindSequenceFunction[%,k]

Moment[ExponentialDistribution[λ],k]

FullSimplify[%%-%]

Find the moments from an MGF:

MomentGeneratingFunction[ExponentialDistribution[λ],t]

Table[Derivative[n][Function[t,(λ/(-t+λ))]][0],{n,5}]

SeriesCoefficient[n!(λ/(-t+λ)),{t,0,n}]

Moment Converters – Formal Moments

Convert between any types of moments:

CellDisplay@Manipulate[TraditionalForm[f[n]]->TraditionalForm@MomentConvert[f[n],g],{{f,Moment,"source"},{Moment,CentralMoment,FactorialMoment,Cumulant}},{{g,CentralMoment,"target"},{Moment,CentralMoment,FactorialMoment,Cumulant}},{{n,2},1,20,1}]

Basic Examples

Represent a cumulant in terms of other moments:

MomentConvert[Cumulant[4],"Moment"]

TraditionalForm[%]

MomentConvert[Cumulant[4],"CentralMoment"]//TraditionalForm

MomentConvert[Cumulant[4],"FactorialMoment"]//TraditionalForm

Convert any polynomial combination of moments:

MomentConvert[%,Cumulant]

Manipulate[TraditionalForm[f[n]]->TraditionalForm@MomentConvert[f[n],g],{{f,Moment,"source"},{Moment,CentralMoment,FactorialMoment,Cumulant}},{{g,CentralMoment,"target"},{Moment,CentralMoment,FactorialMoment,Cumulant}},{n,1,20,1}]

Convert multivariate moments:

MomentConvert[Cumulant[{2,4}],"CentralMoment"]//TraditionalForm

Manipulate[TraditionalForm[f[{m,n}]]->TraditionalForm@MomentConvert[f[{m,n}],g],{{f,Moment,"source"},{Moment,CentralMoment,FactorialMoment,Cumulant}},{{g,CentralMoment,"target"},{Moment,CentralMoment,FactorialMoment,Cumulant}},{m,1,20,1},{n,1,20,1}]

Evaluate a polynomial combination of moments for a distribution:

MomentEvaluate[Moment[1]CentralMoment[4],NormalDistribution[μ,σ]]

Moment Converters – Formal Sample Moments

Convert to any type of sample moment:

CellDisplay@Manipulate[TraditionalForm[f[m]]->TraditionalForm@MomentConvert[f[m],{"SampleEstimator",n}],{{f,CentralMoment,"source"},{Moment,CentralMoment,FactorialMoment,Cumulant}},{{m,2},1,20,1}]

Basic Examples

For sample moments we use power sums (Subscript[S, r])==(∑_(i=1)^n)((Subscript[x, i])^r)

Manipulate[PowerSymmetricPolynomial[r,Array[(Subscript[x, #])&,{10}]],{r,1,10,1}]

Convert to a sample estimator:

MomentConvert[Moment[4],"SampleEstimator"]//TraditionalForm

MomentConvert[Moment[4],{"SampleEstimator",n}]//TraditionalForm

Manipulate[TraditionalForm[f[m]]->TraditionalForm@MomentConvert[f[m],{"SampleEstimator",n}],{{f,Moment,"source"},{Moment,CentralMoment,FactorialMoment,Cumulant}},{m,1,20,1}]

Convert any polynomial combination of moments:

MomentConvert[Cumulant[4]((CentralMoment[2])^2),{"SampleEstimator",n}]//TraditionalForm

Convert multivariate moment combinations:

MomentConvert[CentralMoment[{1,2}],{"SampleEstimator",n}]//TraditionalForm

Examples

Convert formal sample sums to moments (viewing all data points as independent identically distributed):

MomentConvert[Moment[2],"SampleEstimator"]

MomentConvert[%,"Moment"]

For central moments and cumulants using sample sums lead to biased estimators: [Ever wondered about that factor (1/(n-1)) ... for sample variance?]

MomentConvert[CentralMoment[2],{"SampleEstimator",n}]//TraditionalForm

MomentConvert[%,{"CentralMoment",n}]//TraditionalForm

Can systematically generate unbiased estimators:

MomentConvert[CentralMoment[2],{"UnbiasedSampleEstimator",n}]

TraditionalForm[%]

MomentConvert[%,{"CentralMoment",n}]

These estimators are represented in terms of special symmetric polynomials: (Subscript[A, (Subscript[r, 1]),…,(Subscript[r, k])])==∑((Subscript[x, (Subscript[i, 1])])^(Subscript[r, 1]))⋯ ((Subscript[x, (Subscript[i, k])])^(Subscript[r, k]))

AugmentedSymmetricPolynomial[{1,2},{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])}]

Can convert between these and power sums (typically easier to implement power sums for evaluating on large datasets):

MomentConvert[AugmentedSymmetricPolynomial[{1,2}],"PowerSymmetricPolynomial"]//TraditionalForm

MomentConvert[%,"AugmentedSymmetricPolynomial"]//TraditionalForm

MomentConvert[CentralMoment[2],{"UnbiasedSampleEstimator",n},"PowerSymmetricPolynomial"]//FunctionExpand//TraditionalForm

Find an unbiased estimator of product of two cumulants:

est=MomentConvert[Cumulant[2]Cumulant[3],"UnbiasedSampleEstimator"];

mean=MomentConvert[est,"Cumulant"]

Compute the variance for the estimator:

MomentConvert[(((est-mean))^2),{"Cumulant",n}]//TraditionalForm

Get the asymptotics:

Series[%,{n,∞,1}]//TraditionalForm

Distribution Properties – Simulation

(*
data=RandomNumber[NormalDistribution[0,1],10^3];
datap=ListLinePlot[data,Ticks->None,ImageSize->Tiny];
distp=Plot[{0,PDF[NormalDistribution[0,1],x]},{x,-3,3},Filling->Axis,Ticks->None,ImageSize->Tiny];
*)
FramedStatement[HoldForm[RandomNumber[,…]⟹ ]]

Basic Examples

Given a distribution we can simulate it:

data=RandomReal[GumbelDistribution[1,3],10^4];

Show[
Histogram[data,25,"ProbabilityDensity"],
Plot[PDF[GumbelDistribution[1,3],x],{x,-15,8},PlotStyle->Thick]]

The average number of accidents in a city is 100 per day. Simulate the daily accidents:

RandomInteger[PoissonDistribution[100],30]

ListPlot[%,Filling->Axis]

Parametric Distributions – Estimation

(*
data=RandomNumber[NormalDistribution[0,1],10^3];
datap=ListLinePlot[data,Ticks->None,ImageSize->Tiny];
distp=Plot[{0,PDF[NormalDistribution[0,1],x]},{x,-3,3},Filling->Axis,Ticks->None,ImageSize->Tiny];
distfamp=Plot[Evaluate[Table[PDF[NormalDistribution@@θ,x],{θ,{{-1,1/2},{0,1},{1,2/3},{1/2,2},{0,1/3},{2/3,1}}}]],{x,-3,3},Filling->Axis,Ticks->None,PlotRange->All,ImageSize->Tiny];
*)


FramedStatement[HoldForm[DistributionEstimate[,,…]⟹]]

Basic Examples

And given outcomes we can estimate the distribution:

data=RandomInteger[PoissonDistribution[100],10^3];

DistributionEstimate[data,PoissonDistribution[μ]]

There are several different means of comparing data to models, quantile-quantile plots compares the quantiles:

QuantilePlot[data,PoissonDistribution[100.281`]]

Examples

Get the distributions word length for some languages:

languages={"Dutch","English","French","German","Hebrew","Polish","Russian","Spanish","Swedish"};

worddata=Table[StringLength/@DictionaryLookup[{l,All}],{l,languages}];

Estimate binomial models:

binom=Table[DistributionEstimate[i,BinomialDistribution[n,p],MaxIterations->1000],{i,worddata}]

Compare results:

comp=Table[Show[Histogram[worddata[[i]],{Range[25]-1/2},"ProbabilityDensity",PlotLabel->languages[[i]]],DiscretePlot[PDF[binom[[i]],x],{x,0,25},PlotRange->All,PlotStyle->PointSize[.02]]],{i,Length[languages]}];

GraphicsGrid[Partition[comp,3],ImageSize->Large]

Parametric Distributions

Distributions constructed from parameters:

({{"Univariate"}, {"Continuous"}})({{"Univariate"}, {"Discrete"}})({{"Multivariate"}, {"Continuous"}})({{"Multivariate"}, {"Discrete"}})//FrameDisplay//CellDisplay

Coverage

Compared to Mathematica 7:

Compared to some close competitors:

Tooltip[,({{({{({{"Univariate"}, {"Continuous"}})}, {({{({{"", "Mathematica7", "Mathematica8", "Maple", "Matlab", "R"}, {"Weibull", " ", " ", " ", " ", " "}, {"Wald", " ", " ", " ", " ", " "}, {"Uniform", " ", " ", " ", " ", " "}, {"Student-T", " ", " ", " ", " ", " "}, {"Standard Normal", " ", " ", " ", " ", " "}, {"Standard Gamma", " ", " ", " ", " ", " "}, {"Standard Exponential", " ", " ", " ", " ", " "}, {"Sech Square", " ", " ", " ", " ", " "}, {"Normal", " ", " ", " ", " ", " "}, {"Log Normal", " ", " ", " ", " ", " "}, {"Logistic", " ", " ", " ", " ", " "}, {"Inverse Gaussian", " ", " ", " ", " ", " "}, {"Generalized Pareto", " ", " ", " ", " ", " "}, {"Generalized Gamma", " ", " ", " ", " ", " "}, {"Gamma", " ", " ", " ", " ", " "}, {"F-Ratio", " ", " ", " ", " ", " "}, {"Exponential", " ", " ", " ", " ", " "}, {"Error", " ", " ", " ", " ", " "}, {"Erlang", " ", " ", " ", " ", " "}, {"Chi-square", " ", " ", " ", " ", " "}, {"Beta", " ", " ", " ", " ", " "}, {"Rayleigh", " ", " ", " ", " ", " "}, {"Noncentral Student-T", " ", " ", " ", " ", " "}, {"Noncentral F-Ratio", " ", " ", " ", " ", " "}, {"Noncentral Chi-square", " ", " ", " ", " ", " "}, {"Doubly Noncentral F", " ", " ", " ", " ", " "}, {"Triangular", " ", " ", " ", " ", " "}, {"Shifted Pareto", " ", " ", " ", " ", " "}, {"Power", " ", " ", " ", " ", " "}, {"Pareto IV", " ", " ", " ", " ", " "}, {"Pareto III", " ", " ", " ", " ", " "}, {"Pareto II", " ", " ", " ", " ", " "}, {"Pareto I", " ", " ", " ", " ", " "}, {"Pareto", " ", " ", " ", " ", " "}, {"Laplace", " ", " ", " ", " ", " "}, {"Inverse Pareto", " ", " ", " ", " ", " "}, {"Doubly Pareto Uniform", " ", " ", " ", " ", " "}, {"Chotikapanich Generalized Pareto", " ", " ", " ", " ", " "}, {"Cauchy", " ", " ", " ", " ", " "}, {"Bradford", " ", " ", " ", " ", " "}, {"Maxwell", " ", " ", " ", " ", " "}, {"Gumbel Maximum", " ", " ", " ", " ", " "}, {"Gumbel", " ", " ", " ", " ", " "}, {"Extreme Value", " ", " ", " ", " ", " "}}), ({{"", "Mathematica8", "Mathematica7", "R", "Maple", "Matlab"}, {"Inverse Gamma", " ", " ", " ", " ", " "}, {"Scaled Inverse Chi-square", " ", " ", " ", " ", " "}, {"Inverse Chi-square", " ", " ", " ", " ", " "}, {"Half-Normal", " ", " ", " ", " ", " "}, {"Chi", " ", " ", " ", " ", " "}, {"vonMises", " ", " ", " ", " ", " "}, {"Noncentral Beta", " ", " ", " ", " ", " "}, {"Moyal", " ", " ", " ", " ", " "}, {"Johnson SU", " ", " ", " ", " ", " "}, {"Johnson SB", " ", " ", " ", " ", " "}, {"Johnson", " ", " ", " ", " ", " "}, {"Generalized Extreme Value", " ", " ", " ", " ", " "}, {"Rice", " ", " ", " ", " ", " "}, {"Pearson VII", " ", " ", " ", " ", " "}, {"Pearson VI", " ", " ", " ", " ", " "}, {"Pearson V", " ", " ", " ", " ", " "}, {"Pearson IV", " ", " ", " ", " ", " "}, {"Pearson III", " ", " ", " ", " ", " "}, {"Pearson I", " ", " ", " ", " ", " "}, {"Pearson", " ", " ", " ", " ", " "}, {"Nakagami", " ", " ", " ", " ", " "}, {"Log Logistic", " ", " ", " ", " ", " "}, {"Fatigue Life", " ", " ", " ", " ", " "}, {"Birnbaum Saunders", " ", " ", " ", " ", " "}, {"Arcsin", " ", " ", " ", " ", " "}, {"Tukey", " ", " ", " ", " ", " "}, {"Generalized Gaussian", " ", " ", " ", " ", " "}, {"Exponential Power", " ", " ", " ", " ", " "}, {"Wigner Semicircle", " ", " ", " ", " ", " "}, {"Upper Truncated Pareto", " ", " ", " ", " ", " "}, {"Tukey Lambda", " ", " ", " ", " ", " "}, {"Truncated Normal", " ", " ", " ", " ", " "}, {"Truncated Exponential", " ", " ", " ", " ", " "}, {"Student-Z", " ", " ", " ", " ", " "}, {"Stable", " ", " ", " ", " ", " "}, {"Skew Normal", " ", " ", " ", " ", " "}, {"Singh Maddala", " ", " ", " ", " ", " "}, {"Sech", " ", " ", " ", " ", " "}, {"RS Generalized Tukey Lambda", " ", " ", " ", " ", " "}, {"PERT", " ", " ", " ", " ", " "}, {"Paralogistic", " ", " ", " ", " ", " "}, {"Normal Mixture", " ", " ", " ", " ", " "}, {"Map Airy", " ", " ", " ", " ", " "}, {"Makeham", " ", " ", " ", " ", " "}}), ({{"", "Mathematica8", "Mathematica7", "Matlab", "R", "Maple"}, {"Log Gamma", " ", " ", " ", " ", " "}, {"Lindley", " ", " ", " ", " ", " "}, {"Levy Skew Alpha Stable", " ", " ", " ", " ", " "}, {"Kumaraswamy", " ", " ", " ", " ", " "}, {"Irwin Hall", " ", " ", " ", " ", " "}, {"Inverse Burr", " ", " ", " ", " ", " "}, {"Hyperbolic Secant", " ", " ", " ", " ", " "}, {"Histogram", " ", " ", " ", " ", " "}, {"Gompertz Makeham", " ", " ", " ", " ", " "}, {"Gompertz", " ", " ", " ", " ", " "}, {"Gamma Mixture", " ", " ", " ", " ", " "}, {"FMKL Generalized Tukey Lambda", " ", " ", " ", " ", " "}, {"Fisher-Z", " ", " ", " ", " ", " "}, {"Exponential Mixture", " ", " ", " ", " ", " "}, {"ExponentialGamma", " ", " ", " ", " ", " "}, {"Dagum", " ", " ", " ", " ", " "}, {"Censored Normal", " ", " ", " ", " ", " "}, {"Censored Gumbel", " ", " ", " ", " ", " "}, {"Burr XII", " ", " ", " ", " ", " "}, {"Burr III", " ", " ", " ", " ", " "}, {"Borel Tanner", " ", " ", " ", " ", " "}, {"Beta Prime", " ", " ", " ", " ", " "}, {"Beta Mixture", " ", " ", " ", " ", " "}, {"Benini", " ", " ", " ", " ", " "}, {"Bates", " ", " ", " ", " ", " "}, {"Quadratic Form", " ", " ", " ", " ", " "}, {"Piecewise", " ", " ", " ", " ", " "}, {"Wrapped Up Cauchy", " ", " ", " ", " ", " "}, {"Wrapped Normal", " ", " ", " ", " ", " "}, {"Wilcoxon", " ", " ", " ", " ", " "}, {"Tweedie", " ", " ", " ", " ", " "}, {"Skew Laplace", " ", " ", " ", " ", " "}, {"Power Normal", " ", " ", " ", " ", " "}, {"Normal Scores", " ", " ", " ", " ", " "}, {"Maximum F-Ratio", " ", " ", " ", " ", " "}, {"Kruskall Wallis", " ", " ", " ", " ", " "}, {"Kendall Tau", " ", " ", " ", " ", " "}, {"Hjort", " ", " ", " ", " ", " "}, {"Generalized Weibull", " ", " ", " ", " ", " "}, {"Generalized Logistic", " ", " ", " ", " ", " "}, {"Gamma Count", " ", " ", " ", " ", " "}, {"Cramer vonMises", " ", " ", " ", " ", " "}, {"Burr I", " ", " ", " ", " ", " "}})}})}}), ({{({{"Univariate"}, {"Discrete"}})}, {({{"", "Mathematica8", "Mathematica7", "Maple", "R", "Matlab"}, {"Poisson Gamma", " ", " ", " ", " ", " "}, {"Poisson", " ", " ", " ", " ", " "}, {"Negative Binomial", " ", " ", " ", " ", " "}, {"Hypergeometric", " ", " ", " ", " ", " "}, {"Geometric", " ", " ", " ", " ", " "}, {"Coin Tossing", " ", " ", " ", " ", " "}, {"Binomial", " ", " ", " ", " ", " "}, {"Bernoulli", " ", " ", " ", " ", " "}, {"Discrete Uniform", " ", " ", " ", " ", " "}, {"Beta Binomial", " ", " ", " ", " ", " "}, {"Zipf", " ", " ", " ", " ", " "}, {"Zeta", " ", " ", " ", " ", " "}, {"Log Series", " ", " ", " ", " ", " "}, {"Beta Negative Binomial", " ", " ", " ", " ", " "}, {"Zero Inflated Poisson", " ", " ", " ", " ", " "}, {"Zero Altered Poisson", " ", " ", " ", " ", " "}, {"Yule Simon", " ", " ", " ", " ", " "}, {"Yule", " ", " ", " ", " ", " "}, {"Waring Yule", " ", " ", " ", " ", " "}, {"Waring", " ", " ", " ", " ", " "}, {"Truncated Pareto", " ", " ", " ", " ", " "}, {"Skellam", " ", " ", " ", " ", " "}, {"Polya Aeppli", " ", " ", " ", " ", " "}, {"Poisson Mixture", " ", " ", " ", " ", " "}, {"Order", " ", " ", " ", " ", " "}, {"Borel", " ", " ", " ", " ", " "}, {"Tweedie", " ", " ", " ", " ", " "}, {"Power Variance Function Poisson", " ", " ", " ", " ", " "}, {"Multiplicative Poisson", " ", " ", " ", " ", " "}, {"Multiplicative Binomial", " ", " ", " ", " ", " "}, {"Generalized Hypergeometric", " ", " ", " ", " ", " "}, {"Double Poisson", " ", " ", " ", " ", " "}, {"Double Binomial", " ", " ", " ", " ", " "}})}}), ({{({{({{"Multivariate"}, {"Continuous"}})}, {({{"", "Mathematica8", "Mathematica7", "Matlab", "R", "Maple"}, {"Multivariate T", " ", " ", " ", " ", " "}, {"Multinormal", " ", " ", " ", " ", " "}, {"Bivariate Normal", " ", " ", " ", " ", " "}, {"Hotelling T-Square", " ", " ", " ", " ", " "}, {"Bivariate Poisson", " ", " ", " ", " ", " "}, {"T-Copula", " ", " ", " ", " ", " "}, {"Gumbel Copula", " ", " ", " ", " ", " "}, {"Gaussian Mixture", " ", " ", " ", " ", " "}, {"Gaussian Copula", " ", " ", " ", " ", " "}, {"Frank Copula", " ", " ", " ", " ", " "}, {"Clayton Copula", " ", " ", " ", " ", " "}, {"Uniform Sum", " ", " ", " ", " ", " "}, {"Truncated Bivariate Normal", " ", " ", " ", " ", " "}, {"Standard Bivariate Normal", " ", " ", " ", " ", " "}, {"Exponential Mixtures", " ", " ", " ", " ", " "}, {"Copula Student-T", " ", " ", " ", " ", " "}, {"Bivariate Farlie Gumbel Morgenstern", " ", " ", " ", " ", " "}, {"Binormal", " ", " ", " ", " ", " "}, {"Wishart", " ", " ", " ", " ", " "}, {"Inverse Wishart", " ", " ", " ", " ", " "}, {"Simplex", " ", " ", " ", " ", " "}, {"Dirichlet", " ", " ", " ", " ", " "}})}})}, {({{({{"Multivariate"}, {"Discrete"}})}, {({{"", "Mathematica8", "Mathematica7", "R", "Matlab", "Maple"}, {"Multinomial", " ", " ", " ", " ", " "}, {"Negative Multinomial II", " ", " ", " ", " ", " "}, {"Negative Multinomial I", " ", " ", " ", " ", " "}, {"Negative Multinomial", " ", " ", " ", " ", " "}, {"Multivariate Poisson", " ", " ", " ", " ", " "}, {"Multivariate Hypergeometric", " ", " ", " ", " ", " "}, {"Compound Poisson", " ", " ", " ", " ", " "}})}})}})}})]//FrameDisplay//CellDisplay

We intend to make documentation for each one extensive:

{BinomialDistribution, PoissonDistribution, NormalDistribution, ExponentialDistribution}

Meta Distributions

FramedMediumStatement@Grid[{{TransformDistribution,"functions of random variables"},
{OrderDistribution,"(n^th) element of random vectors"},
{CensoredDistribution,"interval censored distribution"},
{TruncatedDistribution,"interval truncated distribution"},
{ComponentMixtureDistribution,"mixture of distributions"},
{IndependentDistribution,"independ multivariate components"},
{CopulaDistribution,"multivariate distribution from kernel and marginals"},
{MarginalDistribution,"lower dimensional distribution"}},
Alignment->{Left,Baseline}]

TruncatedDistribution

Distributions truncated to an interval or cube:

(*
distp=Plot[{0,PDF[NormalDistribution[0,1],x]},{x,-3,3},Filling->Axis,Ticks->None,PlotRange->{0,0.5},ImageSize->Tiny];
truncp=Plot[{0,PDF[TruncatedDistribution[{-1,2},NormalDistribution[0,1]],x]},{x,-3,3},Filling->Axis,Ticks->None,PlotRange->{0,0.5},ImageSize->Tiny];
*)



FramedStatement[HoldForm[TruncatedDistribution[{(Subscript[x, min]),(Subscript[x, max])},]⟹ ]]

Basic Examples

Truncate the values to an interval — or condition on the interval:

𝒟=TruncatedDistribution[{-1,1},NormalDistribution[0,1]]

PDF[𝒟,x]

Plot[{PDF[𝒟,x],PDF[NormalDistribution[0,1],x]},{x,-2,2},Filling->Axis]

CensoredDistribution

Distributions censored to an interval or cube:

(*
distp=Plot[{0,CDF[NormalDistribution[0,1],x]},{x,-3,3},Filling->Axis,Ticks->None,ImageSize->Tiny];
censp=Plot[Evaluate@{0,CDF[CensoredDistribution[{-1,∞},NormalDistribution[0,1]],x]},{x,-3,3},Exclusions->None,Filling->Axis,Ticks->None,ImageSize->Tiny];
*)

FramedStatement[HoldForm[CensoredDistribution[{(Subscript[x, min]),(Subscript[x, max])},]⟹]]

Basic Examples

A left-censored distribution:

𝒟=CensoredDistribution[{a,∞},NormalDistribution[0,1]]

CDF[𝒟,x]

Manipulate[Plot[%,{x,-3,3},Exclusions->None,Filling->Axis],{{a,-1},-3,3}]

A right-censored distribution:

𝒟=CensoredDistribution[{-∞,b},NormalDistribution[0,1]]

CDF[𝒟,x]

Manipulate[Plot[%,{x,-3,3},Exclusions->None,Filling->Axis],{{b,1},-3,3}]

IndependentDistribution

Distribution built from product or lower dimensional distributions:

(*
ind1p=Plot[{0,PDF[NormalDistribution[0,1],x]},{x,-2,2},Filling->Axis,Ticks->None,ImageSize->Tiny];
ind2p=Plot[{0,PDF[LaplaceDistribution[0,1],x]},{x,-2,2},Filling->Axis,Ticks->None,ImageSize->Tiny];
indp=Plot3D[{0,PDF[IndependentDistribution[NormalDistribution[0,1],LaplaceDistribution[0,1]],{x,y}]},{x,-2,2},{y,-2,2},PlotStyle->Lighter[ColorData[1,2],0.5],Filling->Axis,FillingStyle->Lighter[ColorData[1,2],0.8],MeshFunctions->{#3&},Lighting->"Neutral",Ticks->None,PlotRange->All,ImageSize->Tiny];
*)



FramedStatement[HoldForm[IndependentDistribution[,]⟹ ]]

Basic Examples

IndependentDistribution – distribution for independent random variables

𝒟=IndependentDistribution[NormalDistribution[0,1],LaplaceDistribution[0,1]]

PDF[𝒟,{x,y}]

Plot3D[%,{x,-2,2},{y,-2,2}]

MarginalDistribution

Lower dimensional distribution by extractig marginal distribution:

(*
distp=Show[Plot3D[CDF[BinormalDistribution[1/2],{x,y}],{x,-2,2},{y,-2,2},PlotStyle->Lighter[ColorData[1,2],0.5],Filling->Axis,FillingStyle->Lighter[ColorData[1,2],0.8],MeshFunctions->{#3&},Lighting->"Neutral",Ticks->None,PlotRange->All,ImageSize->Tiny],ParametricPlot3D[{{u,2,CDF[BinormalDistribution[1/2],{u,2}]},{2,u,CDF[BinormalDistribution[1/2],{2,u}]}},{u,-2,2},PlotStyle->{Yellow,Green}]/.Line[pts_]->Tube[pts,0.05],PlotRange->{0,1}];
margp=Plot[{0,CDF[BinormalDistribution[1/2],{u,2}]},{u,-2,2},Filling->Axis,Ticks->None,ImageSize->Tiny];
*)


FramedStatement[HoldForm[MarginalDistribution[,…]⟹ ]]

Basic Examples

Suppose we have a multivariate distribution built from univariate:

𝒟=IndependentDistribution[NormalDistribution[μ,σ],ExponentialDistribution[λ],ChiSquareDistribution[ν]];

These are 1D marginal distributions:

Table[MarginalDistribution[𝒟,i],{i,3}]//Column

These are 2D marginal distributions:

Table[MarginalDistribution[𝒟,Complement[{1,2,3},{i}]],{i,3}]//Column

In some case marginals are known to be other special distributions:

{MarginalDistribution[BinormalDistribution[ρ],1],MarginalDistribution[BinormalDistribution[ρ],2]}

{MarginalDistribution[MultinomialDistribution[n, {p,q}], 1],MarginalDistribution[MultinomialDistribution[n, {p,q}], 2]}

In other cases there is no special result, but it still a valid distribution:

𝒟=MarginalDistribution[MultinomialDistribution[4, {1/2,1/3,1/6}], {1,2}]

PDF[𝒟, {x,y}]

DiscretePlot3D[%,{x,0,4},{y,0,4},PlotRange->All]

Meta Distributions – ComponentMixtureDistribution

Distribution built from a convex combination of component distributions:

(*
mix1p=Plot[{0,PDF[NormalDistribution[1,1/2],x]},{x,-1,3},Filling->Axis,Ticks->None,ImageSize->Tiny,PlotRange->{0,2.5}];
mix2p=Plot[{0,PDF[NormalDistribution[2,1/6],x]},{x,-1,3},Filling->Axis,Ticks->None,ImageSize->Tiny,PlotRange->{0,2.5}];
mixp=Plot[{0,PDF[ComponentMixtureDistribution[{3/5,2/5}, {NormalDistribution[1,1/2],NormalDistribution[2,1/6]}],x]},{x,-1,3},Filling->Axis,Ticks->None,ImageSize->Tiny,PlotRange->{0,2.5}];
*)



FramedStatement[HoldForm[ComponentMixtureDistribution[{(Subscript[λ, 1]),(Subscript[λ, 2])},{,}]⟹ ]]

Basic Examples

𝒟=ComponentMixtureDistribution[{3/5,2/5}, {NormalDistribution[1,1/2],NormalDistribution[2,1/6]}]

PDF[𝒟,x]

Plot[{PDF[𝒟,x],PDF[NormalDistribution[1,1/2],x],PDF[NormalDistribution[2,1/6],x]},{x,-1,3},Filling->Axis,PlotRange->All]

Examples

Can produce rich variety of distribution shapes with relatively few components:

Gaussian = NormalDistribution[];
SkewedUnimodal = ComponentMixtureDistribution[{1/5,1/5,3/5},{NormalDistribution[],NormalDistribution[1/2,2/3],NormalDistribution[13/12,5/9]}];
StrongSkew = ComponentMixtureDistribution[ConstantArray[1/8,8],Table[NormalDistribution[3((2/3)^i - 1),(2/3)^i],{i,0,7}]];
KurtoticUnimodal = ComponentMixtureDistribution[{2/3,1/3},{NormalDistribution[],NormalDistribution[0,1/10]}];
Outlier = ComponentMixtureDistribution[{1/10,9/10},{NormalDistribution[],NormalDistribution[0,1/10]}];
Bimodal = ComponentMixtureDistribution[{1/2,1/2},{NormalDistribution[-1,2/3],NormalDistribution[1,2/3]}];
SeparatedBimodal = ComponentMixtureDistribution[{1/2,1/2},{NormalDistribution[-3/2,1/2],NormalDistribution[3/2,1/2]}];
SkewedBimodal = ComponentMixtureDistribution[{3/4,1/4},{NormalDistribution[],NormalDistribution[3/2,1/3]}];
Trimodal = ComponentMixtureDistribution[{9/20,9/20,1/10},{NormalDistribution[-6/5,3/5],NormalDistribution[6/5,3/5],NormalDistribution[0,1/4]}];
Claw = ComponentMixtureDistribution[{1/2,1/10,1/10,1/10,1/10,1/10},Join[{NormalDistribution[]},Table[NormalDistribution[i/2 - 1, 1/10],{i,0,4}]]];
DoubleClaw = ComponentMixtureDistribution[Join[{49/100,49/100},ConstantArray[1/350,7]],Join[{NormalDistribution[-1,2/3],NormalDistribution[1,2/3]},Table[NormalDistribution[(i-3)/2,1/100],{i,0,6}]]];
AsymmetricClaw = ComponentMixtureDistribution[Join[{1/2},Table[(2^(1-i))/31,{i,-2,2}]],Join[{NormalDistribution[]},Table[NormalDistribution[i + 1/2, (2^-i)/10],{i,-2,2}]]];
AsymmetricDoubleClaw = ComponentMixtureDistribution[{46/100,46/100,1/300,1/300,1/300,7/300,7/300,7/300},Join[Join[Table[NormalDistribution[2i - 1,2/3],{i,0,1}],Table[NormalDistribution[-i/2,1/100],{i,1,3}]],Table[NormalDistribution[i/2,7/100],{i,1,3}]]];
SmoothComb = ComponentMixtureDistribution[Table[(2^(5-i))/63,{i,0,5}],Table[NormalDistribution[(65-(96(1/2)^i))/21,(32/63)/2^i],{i,0,5}]];
DiscreteComb = ComponentMixtureDistribution[{2/7,2/7,2/7,1/21,1/21,1/21},Join[Table[NormalDistribution[(12i - 15)/7,2/7],{i,0,2}],Table[NormalDistribution[2i/7,1/21],{i,8,10}]]];

dists={Gaussian,SkewedUnimodal,StrongSkew,KurtoticUnimodal,Outlier,Bimodal,SeparatedBimodal,SkewedBimodal,Trimodal,Claw,DoubleClaw,AsymmetricClaw,AsymmetricDoubleClaw,SmoothComb,DiscreteComb};

plots=Table[Plot[PDF[𝒟,x],{x,-3,3},Filling->Axis,FillingStyle->Directive[Opacity[0.4],Lighter[ColorData[1,2],0.5]],PlotStyle->ColorData[1,2],Exclusions->None,Evaluated->True,PlotPoints->100,PlotRange->All,AxesOrigin->{0,0}],{𝒟,dists}];

GraphicsGrid[Partition[plots,3],ImageSize->Large]

Meta Distributions – OrderDistribution

Distribution built by ordering the values from several distribution and picking the (k^th) one:

(*
dist=Plot[{0,PDF[NormalDistribution[],x]},{x,-3,3},Filling->Axis,Ticks->None,ImageSize->Tiny];
odistp=Plot[{0,PDF[OrderDistribution[{NormalDistribution[],2},2],x]},{x,-3,3},Filling->Axis,Ticks->None,ImageSize->Tiny];
*)


FramedStatement[HoldForm[OrderDistribution[{,2},2]⟹ ]]

Basic Examples

Lets define a generalization of Min and Max called Nth:

Nth[list_List,k_]:=Sort[list][[k]]

Permutations[{1,2,3}]

Map[Nth[#,2]&,%]

Lets look at the data from picking out the 2nd of 5 elements:

data=Map[Nth[#,2]&,RandomNumber[LaplaceDistribution[0,1],{10^3,5}]];

hist=Histogram[ data,Automatic,"ProbabilityDensity"]

We can compute the distribution exactly:

PDF[OrderDistribution[{LaplaceDistribution[0,1],5},2],x]

pdf=Plot[%,{x,-4,1}]

Show[hist,pdf]

These are the other ones:

Plot[Evaluate@Table[PDF[OrderDistribution[{LaplaceDistribution[0,1],5},i],x],{i,5}],{x,-2,2},Filling->Axis]

Meta Distributions – TransformDistribution

Distribution constructed by transforming random variables:

(*
distp=Plot[{0,CDF[UniformDistribution[{0,1}],x]},{x,0,2},Filling->Axis,Ticks->None,ImageSize->Tiny];
𝒟1=TransformDistribution[(X^2),X∈UniformDistribution[{0,1}]];
tdistp=Plot[Evaluate@{0,CDF[𝒟1,x]},{x,0,2},Filling->Axis,Ticks->None,ImageSize->Tiny];
*)


FramedStatement[HoldForm[TransformDistribution[g[x],x∈]⟹ ]]

Basic Examples

TransformDistribution – Distribution of f[X] where X∈𝒟

𝒟=TransformDistribution[(X^2),X∈UniformDistribution[{0,1}]]

Get all the properties:

{PDF[𝒟,x],CDF[𝒟,x]}

{Plot[Evaluate@PDF[𝒟,x],{x,0,2},Filling->Axis],Plot[Evaluate@CDF[𝒟,x],{x,0,2},Filling->Axis,PlotRange->{0,1}]}

Random numbers:

Histogram@RandomNumber[𝒟,(10^4)]

Many special relations for distributions:

TransformDistribution[X+Y,X∈PoissonDistribution[(Subscript[μ, 1])]∧Y∈PoissonDistribution[(Subscript[μ, 2])]]

TransformDistribution[(X^2)+(Y^2),X∈NormalDistribution[]∧Y∈NormalDistribution[]]

TransformDistribution[Exp[X],X∈NormalDistribution[μ,σ]]

Meta Distributions – CopulaDistribution

Distribution built from bounded kernel distribution and marginal distributions:

(*
subkp=Plot3D[Evaluate@CDF[CopulaDistribution["Maximal",{UniformDistribution[{0,1}],UniformDistribution[{0,1}]}],{x,y}],{x,0,1},{y,0,1},PlotStyle->Lighter[ColorData[1,2],0.5],Filling->Axis,FillingStyle->Lighter[ColorData[1,2],0.8],MeshFunctions->{#3&},Lighting->"Neutral",Ticks->None,ImageSize->Tiny];
marg1p=Plot[{0,CDF[ExponentialDistribution[2],x]},{x,-1,3},Filling->Axis,Ticks->None,ImageSize->Tiny];
marg2p=Plot[{0,CDF[NormalDistribution[],x]},{x,-3,3},Filling->Axis,Ticks->None,ImageSize->Tiny];
copp=Plot3D[Evaluate@CDF[CopulaDistribution["Maximal",{ExponentialDistribution[2],NormalDistribution[]}],{x,y}],{x,-1,3},{y,-3,3},PlotStyle->Lighter[ColorData[1,2],0.5],Filling->Axis,FillingStyle->Lighter[ColorData[1,2],0.8],MeshFunctions->{#3&},Lighting->"Neutral",Ticks->None,PlotPoints->35,ImageSize->Tiny];
*)



FramedStatement[HoldForm[CopulaDistribution[,{,}]⟹ ]]

Basic Examples

CopulaDistribution – a multivariate CDF C[(Subscript[F, 1]),…,(Subscript[F, n])] with subkernel C and marginal distributions (Subscript[F, i])[(Subscript[x, i])]

Fmax=Refine[CDF[CopulaDistribution["Maximal",{UniformDistribution[{0,1}],UniformDistribution[{0,1}]}],{x,y}],0<x<1&&0<y<1]

Fmin=Refine[CDF[CopulaDistribution["Minimal",{UniformDistribution[{0,1}],UniformDistribution[{0,1}]}],{x,y}],0<x<1&&0<y<1]

minmax=Plot3D[{Fmax,Fmin},{x,0,1},{y,0,1},PlotLabel->"Frechet Lower/Upper Bound",ImageSize->Medium,PlotStyle->Opacity[0.5],Mesh->10,MeshFunctions->{#3&},MeshShading->{Yellow,None,Green,None}]

Lots of different copula families to model dependency structures:

ParallelTable[Plot3D[Evaluate@CDF[CopulaDistribution[c,{UniformDistribution[{0,1}],UniformDistribution[{0,1}]}],{x,y}],{x,0,1},{y,0,1},PlotLabel->c],{c,{"Minimal","Maximal","Product",{"Frank",1/2},{"Clayton",5},{"GumbelHougaard",2},{"FGM",1/3},{"AMH",1/2},{"Binormal",1/2}}}]

Use all possible univariate marginal distributions:

ParallelTable[Plot3D[Evaluate@CDF[CopulaDistribution[{"Frank",1/2},{UniformDistribution[{0,1}],d}],{x,y}],{x,0,1},{y,-3,3},PlotLabel->d],{d,{UniformDistribution[{0,1}],ExponentialDistribution[2],NormalDistribution[0,1],LaplaceDistribution[0,1],GumbelDistribution[1,2]}}]

User Distributions

Distributions constructed from formulas:

(*
𝒟=ProbabilityDistribution[(Piecewise[{{(E^(x λ)) p λ, x<0}, {(E^(-x λ)) (1-p) λ, x≥0}, {0, True}}]),{x,-∞,∞},Assumptions->λ>0∧0<p<1];
Block[{p=1/3},Plot[Evaluate@Table[PDF[𝒟,x],{λ,{0.5,1}}],{x,-3,3},Filling->Axis,Ticks->None,ImageSize->Tiny]];
*)


FramedStatement[HoldForm[ProbabilityDistribution[(Piecewise[{{(E^(x λ)) p λ, x<0}, {(E^(-x λ)) (1-p) λ, x≥0}}]),…]⟹ ]]

Basic Examples

Define a new distribution by giving the PDF:

𝒟=ProbabilityDistribution[(((Sqrt[2])/π))(1/(1+(x^4))),{x,-∞,∞}]

Behaves just like any other distribution:

Moment[𝒟,2]

CDF[𝒟,x]

Probability[X^2<3,X∈𝒟]

Can include parameters, in this case for a two-sided exponential distribution:

𝒟=ProbabilityDistribution[(Piecewise[{{(E^(x λ)) p λ, x<0}, {(E^(-x λ)) (1-p) λ, x≥0}}]),{x,-∞,∞},Assumptions->λ>0∧0<p<1];

CDF[𝒟,x]

Probability[-1/2<X<1/2,X∈𝒟]

{Block[{p=1/2},Plot[Evaluate@Table[PDF[𝒟,x],{λ,{0.5,1,2}}],{x,-3,3},Filling->Axis,PlotRange->All]],Block[{p=1/3},Plot[Evaluate@Table[PDF[𝒟,x],{λ,{0.5,1,2}}],{x,-3,3},Filling->Axis,PlotRange->All]]}

A reflected gamma distribution:

ReflectedGammaDistribution[α_]=ProbabilityDistribution[Piecewise[{{PDF[GammaDistribution[α,1],x]/2,x≥0},{PDF[GammaDistribution[α,1],-x]/2,x<0}}],{x,-∞,∞},Assumptions->α>0];

pdf=PDF[ReflectedGammaDistribution[α],x]

cdf=CDF[ReflectedGammaDistribution[α],x]

CellDisplay@Manipulate[Plot[{(Piecewise[{{(1/2) ((Piecewise[{{(((E^(-x)) (x^(-1+α)))/(Gamma[α])), x>0}, {0, True}}])), x≥0}, {(1/2) ((Piecewise[{{(((E^x) (((-x))^(-1+α)))/(Gamma[α])), -x>0}, {0, True}}])), x<0}, {0, True}}]),(Piecewise[{{(1/2), x==0}, {((Gamma[α,-x])/(2 Gamma[α])), x<0}, {((2 Gamma[α]-Gamma[α,x])/(2 Gamma[α])), True}}])},{x,-7,7},PlotRange->{0,1},ImageSize->Large],{α,0.02,3}]

Gumbel bivariate exponential:

𝒟=ProbabilityDistribution[Exp[-2(x+y)](Exp[x+y]+α(Exp[x]-2)(Exp[y]-2)),{x,0,∞},{y,0,∞},Assumptions->-1<α<1];

CDF[𝒟,{x,y}]

Probability[(X^2)Y<1,{X,Y}∈𝒟]

N[%]

Manipulate[Plot3D[(E^(-2 (x+y))) ((E^(x+y))+(-2+(E^x)) (-2+(E^y)) α),{x,0,2},{y,0,2},MeshFunctions->{#3&},MeshShading->{Green,None,Yellow,None},PlotPoints->35,ImageSize->Large,PlotRange->{0,1}],{α,-0.98,0.98}]

Data Distributions – HistogramDistribution

Distributions constructed from data by binning:

(* 
𝒟ref=ComponentMixtureDistribution[{3/4,1/4},{NormalDistribution[],NormalDistribution[3/2,1/3]}];
datap=ListLinePlot[RandomNumber[𝒟ref,10^3],Ticks->None,ImageSize->Tiny];
𝒟dist=Quiet@HistogramDistribution[RandomNumber[𝒟ref,10^4],{-3,3,0.1}];
distp=Plot[{0,PDF[𝒟dist,x]},{x,-3,3},Filling->Axis,PlotPoints->200,Ticks->None,ImageSize->Tiny]; 
*)



FramedStatement[HoldForm[HistogramDistribution[,…]⟹ ]]

Basic Examples

Normally we think of histogram as a graphical technique for approximately visualizing distributions:

Histogram[RandomNumber[NormalDistribution[],10^4]]

Histogram3D[RandomNumber[NormalDistribution[],{10^4,2}]]

But we can just as well think of it as an actual distribution:

𝒟=HistogramDistribution[RandomNumber[NormalDistribution[],10^4],5]

All the properties taht we have for a parametric distribution works, including distribution related functions:

{PDF[𝒟,x],CDF[𝒟,x],SurvivalFunction[𝒟,x]}

Plot[%,{x,-5,5},PlotPoints->100,Filling->Axis]

Quantile function:

Quantile[𝒟,p]

Plot[%,{p,0,1}]

Moment related functions:

Table[Moment[𝒟,n],{n,5}]

Table[CentralMoment[𝒟,n],{n,5}]

Table[FactorialMoment[𝒟,n],{n,5}]

MomentGeneratingFunction[𝒟,t]

Plot[%,{t,-2,2}]

Probabilities and Expectations:

NProbability[X^2<X,X∈𝒟]

NExpectation[Exp[-Abs[X]],X∈𝒟]

Random numbers:

RandomNumber[𝒟,5]

Histogram@RandomNumber[𝒟,10^4]

With more bins in the histogram distribution we get a closer approximation to the underlying distribution:

𝒟=HistogramDistribution[RandomNumber[NormalDistribution[],10^4]];

Plot[PDF[𝒟,x],{x,-5,5},PlotPoints->150,Filling->Axis]

Histogram[RandomNumber[𝒟,10^4]]

Data Distributions – HistogramDistribution

(*
𝒟ref=ComponentMixtureDistribution[{1/3,2/3},{BinormalDistribution[{-1,-1},{1,1},1/2],BinormalDistribution[{1,1},{1,1},-1/2]}];
sdata=Graphics[{ColorData[1,1],Point@RandomNumber[𝒟ref,10^3]},Ticks->None,ImageSize->Tiny]
data=RandomNumber[𝒟ref,10^4];
𝒟dist=Quiet@HistogramDistribution[RandomNumber[𝒟ref,10^4],{{-3,3,0.5},{-3,3,0.5}}];
distp=Histogram3D[{{{0,0},{0,1}},data},{{-3,3,0.5},{-3,3,0.5}},Ticks->None]
*)


FramedStatement[HoldForm[HistogramDistribution[,…]⟹ ]]

Basic Examples

Work in any number of dimensions:

𝒟=HistogramDistribution[RandomNumber[NormalDistribution[],{10^4,2}]];

All properties work:

{Plot3D[PDF[𝒟,{x,y}],{x,-3,3},{y,-3,3},PlotPoints->50,MaxRecursion->0],Plot3D[CDF[𝒟,{x,y}],{x,-3,3},{y,-3,3}],Plot3D[SurvivalFunction[𝒟,{x,y}],{x,-3,3},{y,-3,3}]}

Moment related properties:

Table[Moment[𝒟,{i,j}],{i,3},{j,3}]

Table[CentralMoment[𝒟,{i,j}],{i,3},{j,3}]

Table[FactorialMoment[𝒟,{i,j}],{i,3},{j,3}]

Probabilities & Expectation:

NProbability[X^2<1&&Y>-1/2,{X,Y}∈𝒟]

Random numbers:

RandomNumber[𝒟,5]

Histogram3D[RandomNumber[𝒟,10^4]]

Marginal distributions:

𝒟x=MarginalDistribution[𝒟,1]

Plot[PDF[𝒟x,x],{x,-3,3},PlotPoints->100,Filling->Axis]

Data Distributions – SmoothKernelDistribution

Distributions constructed from data by smoothing:

(* 
𝒟ref=ComponentMixtureDistribution[{3/4,1/4},{NormalDistribution[],NormalDistribution[3/2,1/3]}];
datap=ListLinePlot[RandomNumber[𝒟ref,500],Ticks->None,ImageSize->Tiny];
𝒟data=SmoothKernelDistribution[RandomNumber[𝒟ref,10^5]];
distp=Plot[{0,PDF[𝒟data,x]},{x,-3,3},Filling->Axis,Ticks->None,ImageSize->Tiny];
*)


FramedStatement[HoldForm[SmoothKernelDistribution[,…]⟹ ]]

Basic Examples

Suppose we have data coming from the following underlying distribution:

𝒟ref=ComponentMixtureDistribution[{3/4,1/4},{NormalDistribution[],NormalDistribution[3/2,1/3]}];

Plot[PDF[𝒟ref,x],{x,-3,3}]

data=RandomNumber[𝒟ref,(10^6)];

We can generate a smooth distribution estimate through (∑_(i=1)^n)K[((x-(Subscript[x, i]))/h)] where K is the smoothing kernel and h is the bandwidth:

𝒟data=SmoothKernelDistribution[data];

Plot[{PDF[𝒟ref,x],PDF[𝒟data,x]},{x,-3,3}]

Resulting distribution functions:

Table[Plot[f[𝒟data,x],{x,-3,3},PlotLabel->f],{f,{PDF,CDF,SurvivalFunction,HazardFunction}}]

Plot[Quantile[𝒟data,p],{p,0,1}]

Moment related properties:

Table[Moment[𝒟data,n],{n,4}]

Table[CentralMoment[𝒟data,n],{n,4}]

Table[FactorialMoment[𝒟data,n],{n,4}]

Table[Cumulant[𝒟data,n],{n,4}]

Probabilities & Expectations:

NProbability[(X^2)<1,X∈𝒟data]

NExpectation[Abs[X-Mean[𝒟data]],X∈𝒟data]

Random numbers:

RandomNumber[𝒟data,5]

Histogram[RandomNumber[𝒟data,10^5]]

Data Distributions – SmoothKernelDistribution

(* 𝒟ref=ComponentMixtureDistribution[{1/3,2/3},{BinormalDistribution[{-1,-1},{1,1},1/2],BinormalDistribution[{1,1},{1,1},-1/2]}];
datap=Graphics[{PointSize[Tiny],Point@RandomNumber[𝒟ref,10^3]},Ticks->None,ImageSize->Tiny];
𝒟dist=SmoothKernelDistribution[RandomNumber[𝒟ref,10^5]];
distp=Plot3D[Evaluate@PDF[𝒟dist,{x,y}],{x,-3,3},{y,-3,3},PlotStyle->Lighter[ColorData[1,2],0.5],Filling->Axis,FillingStyle->Lighter[ColorData[1,2],0.8],Lighting->"Neutral",Mesh->None,PlotPoints->50,Evaluated->True,ImageSize->Tiny,Ticks->None]; *)

FramedStatement[HoldForm[SmoothKernelDistribution[,…]⟹ ]]

Basic Examples

Work in any number of dimensions:

𝒟ref=ComponentMixtureDistribution[{1/3,2/3},{BinormalDistribution[{-1,-1},{1,1},1/2],BinormalDistribution[{1,1},{1,1},-1/2]}];

Plot3D[PDF[𝒟ref,{x,y}],{x,-3,3},{y,-3,3},PlotPoints->35]

data=RandomNumber[𝒟ref,(10^5)];

Get a smooth density distribution:

𝒟data=SmoothKernelDistribution[data];

Plot3D[PDF[𝒟data,{x,y}],{x,-3,3},{y,-3,3},PlotPoints->35]

Moment related:

Table[Moment[𝒟data,{i,j}],{i,3},{j,3}]

Table[CentralMoment[𝒟data,{i,j}],{i,3},{j,3}]

Table[Cumulant[𝒟data,{i,j}],{i,3},{j,3}]

Probabilities & Expectations:

NProbability[X<0&&Y>2,{X,Y}∈𝒟data]

NExpectation[Exp[X Y],{X,Y}∈𝒟data]

Random numbers:

RandomNumber[𝒟data,5]

Examples

Old Faithful geyser data: {waiting time, duration}

OldFaithfulData = {{3.6`,79},{1.8`,54},{3.333`,74},{2.283`,62},{4.533`,85},{2.883`,55},{4.7`,88},{3.6`,85},{1.95`,51},{4.35`,85},{1.833`,54},{3.917`,84},{4.2`,78},{1.75`,47},{4.7`,83},{2.167`,52},{1.75`,62},{4.8`,84},{1.6`,52},{4.25`,79},{1.8`,51},{1.75`,47},{3.45`,78},{3.067`,69},{4.533`,74},{3.6`,83},{1.967`,55},{4.083`,76},{3.85`,78},{4.433`,79},{4.3`,73},{4.467`,77},{3.367`,66},{4.033`,80},{3.833`,74},{2.017`,52},{1.867`,48},{4.833`,80},{1.833`,59},{4.783`,90},{4.35`,80},{1.883`,58},{4.567`,84},{1.75`,58},{4.533`,73},{3.317`,83},{3.833`,64},{2.1`,53},{4.633`,82},{2,59},{4.8`,75},{4.716`,90},{1.833`,54},{4.833`,80},{1.733`,54},{4.883`,83},{3.717`,71},{1.667`,64},{4.567`,77},{4.317`,81},{2.233`,59},{4.5`,84},{1.75`,48},{4.8`,82},{1.817`,60},{4.4`,92},{4.167`,78},{4.7`,78},{2.067`,65},{4.7`,73},{4.033`,82},{1.967`,56},{4.5`,79},{4,71},{1.983`,62},{5.067`,76},{2.017`,60},{4.567`,78},{3.883`,76},{3.6`,83},{4.133`,75},{4.333`,82},{4.1`,70},{2.633`,65},{4.067`,73},{4.933`,88},{3.95`,76},{4.517`,80},{2.167`,48},{4,86},{2.2`,60},{4.333`,90},{1.867`,50},{4.817`,78},{1.833`,63},{4.3`,72},{4.667`,84},{3.75`,75},{1.867`,51},{4.9`,82},{2.483`,62},{4.367`,88},{2.1`,49},{4.5`,83},{4.05`,81},{1.867`,47},{4.7`,84},{1.783`,52},{4.85`,86},{3.683`,81},{4.733`,75},{2.3`,59},{4.9`,89},{4.417`,79},{1.7`,59},{4.633`,81},{2.317`,50},{4.6`,85},{1.817`,59},{4.417`,87},{2.617`,53},{4.067`,69},{4.25`,77},{1.967`,56},{4.6`,88},{3.767`,81},{1.917`,45},{4.5`,82},{2.267`,55},{4.65`,90},{1.867`,45},{4.167`,83},{2.8`,56},{4.333`,89},{1.833`,46},{4.383`,82},{1.883`,51},{4.933`,86},{2.033`,53},{3.733`,79},{4.233`,81},{2.233`,60},{4.533`,82},{4.817`,77},{4.333`,76},{1.983`,59},{4.633`,80},{2.017`,49},{5.1`,96},{1.8`,53},{5.033`,77},{4,77},{2.4`,65},{4.6`,81},{3.567`,71},{4,70},{4.5`,81},{4.083`,93},{1.8`,53},{3.967`,89},{2.2`,45},{4.15`,86},{2,58},{3.833`,78},{3.5`,66},{4.583`,76},{2.367`,63},{5,88},{1.933`,52},{4.617`,93},{1.917`,49},{2.083`,57},{4.583`,77},{3.333`,68},{4.167`,81},{4.333`,81},{4.5`,73},{2.417`,50},{4,85},{4.167`,74},{1.883`,55},{4.583`,77},{4.25`,83},{3.767`,83},{2.033`,51},{4.433`,78},{4.083`,84},{1.833`,46},{4.417`,83},{2.183`,55},{4.8`,81},{1.833`,57},{4.8`,76},{4.1`,84},{3.966`,77},{4.233`,81},{3.5`,87},{4.366`,77},{2.25`,51},{4.667`,78},{2.1`,60},{4.35`,82},{4.133`,91},{1.867`,53},{4.6`,78},{1.783`,46},{4.367`,77},{3.85`,84},{1.933`,49},{4.5`,83},{2.383`,71},{4.7`,80},{1.867`,49},{3.833`,75},{3.417`,64},{4.233`,76},{2.4`,53},{4.8`,94},{2,55},{4.15`,76},{1.867`,50},{4.267`,82},{1.75`,54},{4.483`,75},{4,78},{4.117`,79},{4.083`,78},{4.267`,78},{3.917`,70},{4.55`,79},{4.083`,70},{2.417`,54},{4.183`,86},{2.217`,50},{4.45`,90},{1.883`,54},{1.85`,54},{4.283`,77},{3.95`,79},{2.333`,64},{4.15`,75},{2.35`,47},{4.933`,86},{2.9`,63},{4.583`,85},{3.833`,82},{2.083`,57},{4.367`,82},{2.133`,67},{4.35`,74},{2.2`,54},{4.45`,83},{3.567`,73},{4.5`,73},{4.15`,88},{3.817`,80},{3.917`,71},{4.45`,83},{2,56},{4.283`,79},{4.767`,78},{4.533`,84},{1.85`,58},{4.25`,83},{1.983`,43},{2.25`,60},{4.75`,75},{4.117`,81},{2.15`,46},{4.417`,90},{1.817`,46},{4.467`,74}};

𝒟data=SmoothKernelDistribution[OldFaithfulData];

Plot3D[PDF[𝒟data,{x,y}],{x,1,5.5},{y,35,105},PlotRange->All,ColorFunction->"ThermometerColors"]

One More Thing ...

Meet the team ... and take the opportunity to interact:

team={pavlyk=, gosiak=, darreng=, andyr=, dkapadia=,henryk=};

FramedStatement@Grid@{{Labeled[pavlyk,Text@"Sasha Pavlyk",Top],Labeled[gosiak,Text@"Gosia Konwerska",Top],Labeled[darreng,Text@"Darren Glosemeyer",Top],Labeled[andyr,Text@"Andy Ross",Top],Labeled[dkapadia,Text@"Devendra Kapadia",Top],Labeled[henryk,Text@"Henry Kwong",Top]}}

Initialization

TraditionalForm Typesetting

Setup TraditionalForm rules for moments:

Unprotect[Moment,CentralMoment,FactorialMoment,Cumulant,PowerSymmetricPolynomial,AugmentedSymmetricPolynomial];

Format[Moment[i_List],TraditionalForm]:=Subscript["μ",Sequence@@i]

Format[Moment[i_],TraditionalForm]:=Subscript["μ",i]

Format[CentralMoment[i_List],TraditionalForm]:=Subscript["(OverTilde[μ])",Sequence@@i]

Format[CentralMoment[i_],TraditionalForm]:=Subscript["(OverTilde[μ])",i]

Format[FactorialMoment[i_List],TraditionalForm]:=Subscript["(OverBar[μ])",Sequence@@i]

Format[FactorialMoment[i_],TraditionalForm]:=Subscript["(OverBar[μ])",i]

Format[Cumulant[i_List],TraditionalForm]:=Subscript["κ",Sequence@@i]

Format[Cumulant[i_],TraditionalForm]:=Subscript["κ",i]

Format[PowerSymmetricPolynomial[i_],TraditionalForm]:=Subscript["S",i]

Format[AugmentedSymmetricPolynomial[{i__List}],TraditionalForm]:=Subscript["A",i]

Format[AugmentedSymmetricPolynomial[i_List],TraditionalForm]:=Subscript["A",i]

Protect[Moment,CentralMoment,FactorialMoment,Cumulant,PowerSymmetricPolynomial,AugmentedSymmetricPolynomial];

Examples:

{Moment[n],Moment[{1,2}]}//TraditionalForm

{CentralMoment[n],CentralMoment[{1,2}]}//TraditionalForm

{FactorialMoment[n],FactorialMoment[{1,2}]}//TraditionalForm

{Cumulant[n],Cumulant[{1,2}]}//TraditionalForm

{PowerSymmetricPolynomial[n],PowerSymmetricPolynomial[{1,2}]}//TraditionalForm

{AugmentedSymmetricPolynomial[{1,2}],AugmentedSymmetricPolynomial[{{1,2},{3,4}}]}//TraditionalForm

Display Boxes

CellDisplay: print cell without Out[nn]

CellDisplay[expr_]:=CellPrint[ExpressionCell[expr,"Output",Background->None,ShowCellLabel->False,CellFrame->False]]

FrameDisplay: frame object with standard settings for this document

FrameDisplay[expr_]:=Framed[expr,RoundingRadius->10,Background->Lighter[Gray,0.9],FrameMargins->10,FrameStyle->GrayLevel[0.5]]

Use this for definitions:

DefinitionGrid[{{duv_,dmv_},{cuv_,cmv_}}]:=
FrameDisplay@Grid[
{{SpanFromLeft,Text[Style["Univariate",FontFamily->"Helvetica"]],Text[Style["Multivariate",FontFamily->"Helvetica"]]},
{Text[Style["Discrete",FontFamily->"Helvetica"]],TraditionalForm@duv,TraditionalForm@dmv},
{Text[Style["Continuous",FontFamily->"Helvetica"]],TraditionalForm@cuv,TraditionalForm@cmv}},
Spacings->{2,4},
Dividers->{{False,True,True,False},{False,True,True,False}}
]

FormattedDefinitionGrid[{{duv_,dmv_},{cuv_,cmv_}}]:=
FrameDisplay@Style[Grid[
{{SpanFromLeft,Text[Style["Univariate",FontFamily->"Helvetica"]],Text[Style["Multivariate",FontFamily->"Helvetica"]]},
{Text[Style["Discrete",FontFamily->"Helvetica"]],duv,dmv},
{Text[Style["Continuous",FontFamily->"Helvetica"]],cuv,cmv}},
Spacings->{2,4},
Dividers->{{False,True,True,False},{False,True,True,False}}
],"Text"]

For moment definitions:

MomentDefinitionGrid[{{uvm_,uvgf_},{mvm_,mvgf_}}]:=
FrameDisplay@
Style[Grid[
{{SpanFromLeft,Text[Style["Moment",FontFamily->"Helvetica"]],Text[Style["Generating Function",FontFamily->"Helvetica"]]},
{Text[Style["Univariate",FontFamily->"Helvetica"]],TraditionalForm@uvm,TraditionalForm@uvgf},
{Text[Style["Multivariate",FontFamily->"Helvetica"]],TraditionalForm@mvm,TraditionalForm@mvgf}},
Spacings->{2,4},
Dividers->{{False,True,True,False},{False,True,True,False}}
],"Text"]

For moments:

MomentGrid[list_,type_]:=
Grid[Transpose@{list},Alignment->{Left,Automatic},Background->{None,{{Lighter[Blend[Switch[type,Moment,{Blue,Green},CentralMoment,{Red,Orange},FactorialMoment,{Magenta,Cyan},Cumulant,{Purple,Red}]],.8],White}}},Dividers->{None,{Darker[Gray,.6],{False},Darker[Gray,.6]}},ItemSize->30]

MomentTabView[dist_,n_Integer?Positive]:=
CellDisplay@TabView[{"Moment"->MomentGrid[Table[Moment[dist,k],{k,n}],Moment],"CentralMoment"->MomentGrid[Table[CentralMoment[dist,k],{k,n}],CentralMoment],"FactorialMoment"->MomentGrid[Table[FactorialMoment[dist,k],{k,n}],FactorialMoment],"Cumulant"->MomentGrid[Table[Cumulant[dist,k],{k,n}],Cumulant]}]

For moment generating functions:

MGFTabView[dist_]:=
CellDisplay@TabView[{"MGF"->MomentGeneratingFunction[dist,t],"CMGF"->CentralMomentGeneratingFunction[dist,t],"FMGF"-> FactorialMomentGeneratingFunction[dist,t],"CGF"->CumulantGeneratingFunction[dist,t],"CF"->CharacteristicFunction[dist,t]}]

FramedStatement:

FramedStatement[expr_]:=
CellDisplay@FrameDisplay@Text[Style[expr,Large,FontFamily->"Helvetica"]]

FramedMediumStatement[expr_]:=
CellDisplay@FrameDisplay@Text[Style[expr,FontSize->16,FontFamily->"Helvetica"]]

Various visualization functions:

DistributionPlot, DistributionPlot3D, ListPlot1D, ...

Notes/Todos: 
• Should probably store default plot range for each distribution. The quantile based heuristic works fine in 90% of cases only. 
• Could consider adding support for multiple distributions. 
• Could consider using a PlotLayout option which is either “CommonXAxis” or “Grid” ... 
• Could consider controls for what elements to plot.

DataPlot1D: Useful ways of displaying point clouds.

ListPlot1D[data_List]:=DataPlot1D[data,"JitteredPoints"]

ListPlot1D[data_List,"Points",opts___?OptionQ]:=
Graphics[{PointSize[Tiny],Point@Developer`ToPackedArray@Transpose[{data,ConstantArray[0.,Length[data]]}]},FilterOptions[Graphics,opts],AspectRatio->1/10]

ListPlot1D[data_List,"JitteredPoints",opts___?OptionQ]:=
Graphics[{PointSize[Tiny],Point@Transpose[{data,RandomReal[{-1,1},Length[data]]}]},FilterOptions[Graphics,opts],AspectRatio->1/10]

ListPlot1D[data_List,"Lines",opts___?OptionQ]:=
Graphics[{Thickness[Tiny],Line@Developer`ToPackedArray@Table[{{data[[i]],-1.},{data[[i]],1}},{i,Length[data]}]},FilterOptions[Graphics,opts],AspectRatio->1/10]

DistribubtionPlot related functions:

Clear[DistributionPlot,DistributionPlot3D]

Statistics`Library`ContinuousUnivariateDistributionQ[ExponentialDistribution[2]];

A single univariate distribution:

DistributionPlot[dist_]/;Statistics`Library`UnivariateDistributionQ[dist]:=DistributionPlot[dist,Quantile[dist,{0.05,0.95}]]

Utilities:

BivariateContinuousDistributionQ[dist_]:=Statistics`Library`ContinuousMultivariateDistributionQ[dist]&&Statistics`Library`DistributionDimensionality[dist]==2

BivariateDiscreteDistributionQ[dist_]:=Statistics`Library`DiscreteMultivariateDistributionQ[dist]&&Statistics`Library`DistributionDimensionality[dist]==2

FilterOptions:

FilterOptions[s_Symbol,opts___?OptionQ]:=Sequence@@FilterRules[{opts},Options[s]]

DistributeDefinitions[FilterOptions]

Defaults:

DistributionPlot[dist_,opts___?OptionQ]/;Statistics`Library`UnivariateDistributionQ[dist]:=DistributionPlot[dist,Quantile[dist,{0.05,0.95}],opts]

Continuous univariate distributions:

DistributionPlot[dist_,{xmin_,xmax_},opts___?OptionQ]/;Statistics`Library`ContinuousUnivariateDistributionQ[dist]:=
Module[{x,pdf,cdf,sf,hf,qf},
{pdf,cdf,sf,hf,qf}=Parallelize@ 
{Plot[PDF[dist,x],{x,xmin,xmax},PlotRange->{0,Automatic},PlotLabel->"PDF",PlotStyle->ColorData[1,1],Filling->Axis,FillingStyle->Lighter[ColorData[1,1],0.7],Evaluate@FilterOptions[Plot,opts]],Plot[CDF[dist,x],{x,xmin,xmax},PlotRange->{0,1},PlotLabel->"CDF",PlotStyle->ColorData[1,2],Filling->Axis,FillingStyle->Lighter[ColorData[1,2],0.7],Evaluate@FilterOptions[Plot,opts]],Plot[SurvivalFunction[dist,x],{x,xmin,xmax},PlotRange->{0,1},PlotLabel->"SF",PlotStyle->ColorData[1,3],Filling->Axis,FillingStyle->Lighter[ColorData[1,3],0.7],Evaluate@FilterOptions[Plot,opts]],
Plot[HazardFunction[dist,x],{x,xmin,xmax},PlotRange->{0,Automatic},PlotLabel->"HF",PlotStyle->ColorData[1,4],Filling->Axis,FillingStyle->Lighter[ColorData[1,4],0.7],Evaluate@FilterOptions[Plot,opts]],
Plot[Quantile[dist,p],{p,0,1},PlotLabel->"QF",Evaluate@FilterOptions[Plot,opts]]};

Framed[GraphicsGrid[{{pdf,cdf},{sf,hf} (*,{qf,SpanFromBoth}*)},ImageSize->Large],RoundingRadius->10,Background->Lighter[Gray,0.9]]
]

Discrete univariate distributions:

DistributionPlot[dist_,{xmin_,xmax_},opts___?OptionQ]/;Statistics`Library`DiscreteUnivariateDistributionQ[dist]:=
Module[{k,x,kmin=Floor[xmin],kmax=Ceiling[xmax],pdf,cdf,sf,hf,qf},
DistributeDefinitions[kmin,kmax];
{pdf,cdf,sf,hf,qf}=Parallelize@
{DiscretePlot[PDF[dist,k],{k,kmin,kmax},PlotRange->{{xmin,xmax},Automatic},PlotLabel->"PDF",PlotStyle->ColorData[1,1],Filling->Axis,FillingStyle->Lighter[ColorData[1,1],0.7],Evaluate@FilterOptions[DiscretePlot,opts]],
Plot[CDF[dist,x],{x,xmin,xmax},PlotRange->{0,1},PlotLabel->"CDF",PlotStyle->ColorData[1,2],Filling->Axis,FillingStyle->Lighter[ColorData[1,2],0.7],PlotPoints->100,Evaluate@FilterOptions[Plot,opts]],
Plot[SurvivalFunction[dist,x],{x,xmin,xmax},PlotRange->{0,1},PlotLabel->"SF",PlotStyle->ColorData[1,3],Filling->Axis,FillingStyle->Lighter[ColorData[1,3],0.7],PlotPoints->100,Evaluate@FilterOptions[Plot,opts]],
DiscretePlot[HazardFunction[dist,k],{k,kmin,kmax},PlotRange->{{xmin,xmax},Automatic},PlotLabel->"HF",PlotStyle->ColorData[1,4],Filling->Axis,FillingStyle->Lighter[ColorData[1,4],0.7],Evaluate@FilterOptions[DiscretePlot,opts]],
Plot[Quantile[dist,p],{p,0,1},PlotLabel->"QF",Evaluate@FilterOptions[Plot,opts]]
};

Framed[GraphicsGrid[{{pdf,cdf},{sf,hf} (*,{qf,SpanFromBoth}*)},ImageSize->Large],RoundingRadius->10,Background->Lighter[Gray,0.9]]
]

Continuous bivariate distributions:

DistributionPlot3D[dist_,{xmin_,xmax_},{ymin_,ymax_},opts___?OptionQ]/;BivariateContinuousDistributionQ[dist]:=
Module[{x,y,pdf,cdf,sf,hf},
{pdf,cdf,sf,hf}=Parallelize@
{Plot3D[PDF[dist,{x,y}],{x,xmin,xmax},{y,ymin,ymax},PlotRange->{0,Automatic},PlotLabel->"PDF",PlotStyle->Lighter[ColorData[1,1],0.5],Filling->Axis,FillingStyle->Lighter[ColorData[1,1],0.8],MeshFunctions->{#3&},Lighting->"Neutral",Evaluate@FilterOptions[Plot3D,opts]],
Plot3D[CDF[dist,{x,y}],{x,xmin,xmax},{y,ymin,ymax},PlotRange->{0,1},PlotLabel->"CDF",PlotStyle->Lighter[ColorData[1,2],0.5],Filling->Axis,FillingStyle->Lighter[ColorData[1,2],0.8],MeshFunctions->{#3&},Lighting->"Neutral",Evaluate@FilterOptions[Plot3D,opts]],
Plot3D[SurvivalFunction[dist,{x,y}],{x,xmin,xmax},{y,ymin,ymax},PlotRange->{0,1},PlotLabel->"SF",PlotStyle->Lighter[ColorData[1,3],0.5],Filling->Axis,FillingStyle->Lighter[ColorData[1,3],0.8],MeshFunctions->{#3&},Lighting->"Neutral",Evaluate@FilterOptions[Plot3D,opts]],
Plot3D[HazardFunction[dist,{x,y}],{x,xmin,xmax},{y,ymin,ymax},PlotRange->{0,Automatic},PlotLabel->"HF",PlotStyle->Lighter[ColorData[1,4],0.5],Filling->Axis,FillingStyle->Lighter[ColorData[1,4],0.8],MeshFunctions->{#3&},Lighting->"Neutral",Evaluate@FilterOptions[Plot3D,opts]]};

Framed[GraphicsGrid[{{pdf,cdf},{sf,hf}},ImageSize->Large],RoundingRadius->10,Background->Lighter[Gray,0.9]]
]

Discrete bivariate distributions:

DistributionPlot3D[dist_,{xmin_,xmax_},{ymin_,ymax_},opts___?OptionQ]/;BivariateDiscreteDistributionQ[dist]:=
Module[{x,y,n,m,nmin=Floor[xmin],nmax=Ceiling[xmax],mmin=Floor[ymin],mmax=Ceiling[ymax],pdf,cdf,sf,hf},
DistributeDefinitions[nmin,nmax,mmin,mmax];
{pdf,cdf,sf,hf}=Parallelize@
{DiscretePlot3D[PDF[dist,{n,m}],{n,nmin,nmax},{m,mmin,mmax},PlotRange->{{xmin,xmax},{ymin,ymax},Automatic},PlotLabel->"PDF",PlotStyle->Lighter[ColorData[1,1],0.5],Filling->Axis,FillingStyle->Lighter[ColorData[1,1],0.8],Lighting->"Neutral",Evaluate@FilterOptions[DiscretePlot3D,opts]],
Plot3D[CDF[dist,{x,y}],{x,xmin,xmax},{y,ymin,ymax},PlotRange->{0,1},PlotLabel->"CDF",PlotStyle->Lighter[ColorData[1,2],0.5],Filling->Axis,FillingStyle->Lighter[ColorData[1,2],0.8],MeshFunctions->{#3&},Lighting->"Neutral",PlotPoints->50,MaxRecursion->0,Evaluate@FilterOptions[Plot3D,opts]],
Plot3D[SurvivalFunction[dist,{x,y}],{x,xmin,xmax},{y,ymin,ymax},PlotRange->{0,1},PlotLabel->"SF",PlotStyle->Lighter[ColorData[1,3],0.5],Filling->Axis,FillingStyle->Lighter[ColorData[1,3],0.8],MeshFunctions->{#3&},Lighting->"Neutral",PlotPoints->50,MaxRecursion->0,Evaluate@FilterOptions[Plot3D,opts]],
DiscretePlot3D[HazardFunction[dist,{n,m}],{n,nmin,nmax},{m,mmin,mmax},PlotRange->{{xmin,xmax},{ymin,ymax},Automatic},PlotLabel->"HF",PlotStyle->Lighter[ColorData[1,4],0.5],Filling->Axis,FillingStyle->Lighter[ColorData[1,4],0.8],Lighting->"Neutral",Evaluate@FilterOptions[DiscretePlot3D,opts]]};

Framed[GraphicsGrid[{{pdf,cdf},{sf,hf}},ImageSize->Large],RoundingRadius->10,Background->Lighter[Gray,0.9]]
]

Distribution Converter

PDF[ExtremeValueDistribution[μ,σ],x]== PDF[MaxAsymptoticDistribution[μ,σ,0],x]//Simplify

PDF[FrechetDistribution[α,σ],x]== PDF[MaxAsymptoticDistribution[1/σ,1/(α σ),1/α],x]//Simplify[#,σ>0&&α>0]&

PDF[WeibullDistribution[α,σ],x]==PDF[MinAsymptoticDistribution[σ,σ/α,-1/α],x]//Simplify[#,σ>0&&α>0]&

PDF[GumbelDistribution[μ,σ],x]==PDF[MinAsymptoticDistribution[μ,σ,0],x]//Simplify

DistributionConvert:

DistributionConvert[ExtremeValueDistribution[μ_,σ_]]:=MaxAsymptoticDistribution[μ,σ,0]

DistributionConvert[FrechetDistribution[α_,σ_]]:=MaxAsymptoticDistribution[1/σ,1/(α σ),1/α]

DistributionConvert[WeibullDistribution[α_,σ_]]:=MinAsymptoticDistribution[σ,σ/α,-1/α]

DistributionConvert[GumbelDistribution[μ_,σ_]]:=MinAsymptoticDistribution[μ,σ,0]

Extra:

//FramedStatement

Parametric Distributions – Stable

[NormalDistribution, StableDistribution] [Add plot of convergence to stable ... ]

[Applications: Asymptotic distributions]

Parametric Distributions – Extreme

[GumbelDistribution, FrechetDistribution, WeibullDistribution, MinAsymptoticDistribution, MaxAsymptoticDistribution]

[Application: Asymptotic limits, Reliability, ...] [Add plot of convergence to extreme value distribution ... ]

Extreme value distributions are asymptotic distributions for Max[(Subscript[X, 1]),…,(Subscript[X, n])] or Min[(Subscript[X, 1]),…,(Subscript[X, n])];

• Type I (Gumbel): Unbounded tail, decrease exponentially
• Type II (Frechet): Unbounded tail, decrease slower then exponentially (i.e. not all moments exist)
• Type III (Weibull): Bounded tail

We have introduced:

With min and max there are six cases all put together in

With[{β=2},Plot[Evaluate@Table[PDF[GumbelDistribution[α,β],x],{α,{-1,2,4}}],{x,-8,8},Filling->Axis]]

With[{β=1.5},Plot[Evaluate@Table[PDF[FrechetDistribution[α,β],x],{α,{0.5,1,2}}],{x,0,3},Filling->Axis]]

With[{β=2},Plot[Evaluate@Table[PDF[WeibullDistribution[α,β],x],{α,{0.5,2,4}}],{x,0,5},Filling->Axis]]

All the current ones are special cases:

DistributionConvert[ExtremeValueDistribution[α,β]]

DistributionConvert[GumbelDistribution[α,β]]

DistributionConvert[FrechetDistribution[α,β]]

DistributionConvert[WeibullDistribution[α,β]]

PDF[ExtremeValueDistribution[μ,σ],x]== PDF[MaxAsymptoticDistribution[μ,σ,0],x]//Simplify

PDF[FrechetDistribution[α,σ],x]== PDF[MaxAsymptoticDistribution[1/σ,1/(α σ),1/α],x]//Simplify[#,σ>0&&α>0]&

PDF[WeibullDistribution[α,σ],x]==PDF[MinAsymptoticDistribution[σ,σ/α,-1/α],x]//Simplify[#,σ>0&&α>0]&

PDF[GumbelDistribution[μ,σ],x]==PDF[MinAsymptoticDistribution[μ,σ,0],x]//Simplify

Table[Plot[{CDF[MaxAsymptoticDistribution[-1,1,-1],x],CDF[OrderDistribution[{UniformDistribution[{0,1}],n},n],(1+x/n)]},{x,-2,2}],{n,{2,4,8,16,32}}]

Table[Plot[{CDF[MaxAsymptoticDistribution[-1,1,-1],x],CDF[OrderDistribution[{UniformDistribution[{0,1}],n},n],(1+x/n)]},{x,-2,2}]


Show[Plot[CDF[ExtremeValueDistribution[0,1],x],{x,-2,2},PlotStyle->{Red,Thick,Dashed},PlotRange->All],Plot[CDF[OrderDistribution[{NormalDistribution[0,1],n},n],InverseCDF[NormalDistribution[],1-1/(n)]+x (InverseCDF[NormalDistribution[],1-1/(n E)]-InverseCDF[NormalDistribution[],1-1/(n)])],{x,-2,2},Evaluated->True,PlotRange->All]],{n,{2,4,8,16,32}}]

Table[Plot[CDF[OrderDistribution[{HalfNormalDistribution[1],n},1],x],{x,0,3},PlotRange->All],{n,10}]

Median[OrderDistribution[{HalfNormalDistribution[1],n},1]]

Median[OrderDistribution[{UniformDistribution[{0,1}],n},n]]

Median[MaxAsymptoticDistribution[μ,σ,ξ]]