What’s New in Probability and Statistics

Oleksandr Pavlyk 
Wolfram Research Inc.

Join the Conversation #WolframTechConf

New in Wolfram Language

Extent

Covering 4 releases + work in progress

TimelinePlot[{{DateObject[{2014,12}]->"Mahematica 10.0.2",DateObject[{2015,3}]->"Mathematica 10.1",DateObject[{2015,7}]->"Mathematica 10.2",DateObject[{2015,10}]->"Mathematica 10.3",Interval[{DateObject[{2015,11}],Infinity}]->"Mathematica Dev"},{DateObject[{2014,10}]->Style["WTC 2014",Gray],DateObject[{2015,10}]->Style["WTC 2015",Gray]}},PlotLayout->"Overlayed",PlotLegends->{"Mathematica releases","Wolfram Technology Conferences"},PerformanceGoal->"Speed"]//Quiet

Scope

Time Series

New in 10.2    [ see talk by Gosia Konwerska on Thursday, Oct. 22 at 2 pm]

Direct arithmetic on time series

Efficiency improvements, support for quantity arrays,

Support for irregularly sampled time series

MovingMap — new and powerful

Random Matrices

New in 10.3   [ see talk by Hsien-Ching Kao, Oct. 20, 11:30am ]

Support for common matrix ensemble, in Gaussian, circular, Laguerre

Fast simulation

TracyWidomDistribution

Random Point

RandomPoint is on the intersection of spatial statistics and geometry

Spatial statistics is WIP [see talk by Eduardo Serna, and Hsien-Ching Kao, Oct. 20, 4pm ]

Point patterns

Variogram

Kriging

Forthcoming Support for Quantities

Time Series

Arithmetic

appl=TimeSeries[FinancialData["AAPL",{2000}]]

ibm=TimeSeries[FinancialData["IBM",{2000}]]

Visualize both:

DateListLogPlot[{appl,ibm}]

Ratio:

DateListLogPlot[appl/ibm]

QQ-plot of log-returns:

QuantilePlot[Differences[Log[appl]],Differences[Log[ibm]]]

Time Series

Filtering functions meant to work on regular time series operate in terms of time-increment, or in terms of number of neiboring points.

This does not translate well to the settings of irregular time series. Window sizes, cut-off frequencies, scales, etc. need to be given in absolute units.

Filtering with meaningful window parameters

Support for non-uniformly spaced time series in MovinMap and various filters.

ts=TimeSeries[{{1,0.8},{1.3,0.3},{2,0.1},{(Subscript[t, 4]),2},{3.7,1.4},{4.2,2}}]

Manipulate[Block[{ts=TimeSeries[{{1,0.8},{1.3,0.3},{2,0.1},{t,2},{3.7,1.4},{4.2,2}}]},
ListLinePlot[{ts,MeanFilter[ts,.5]},Mesh->All,Epilog->{Directive[Dashed,Gray,Thin],Line[{{t,-1},{t,3}}]},Prolog->{EdgeForm[GrayLevel[0.85]],FaceForm[Directive[GrayLevel[0.95]]],Rectangle[{t-.5,-1},{t+.5,3}]},PlotRange->{Automatic,{-0.2,2.2}},PlotLegends->{"Original ts with one movable point","MeanFilter[ts,.5]"}]
],{{t,3,"(Subscript[t, 4])"},2.2,3.6,Appearance->"Labeled"}]

Requires support for absolute filter specs, not relative to minimum time-increments.

ts1000=With[{n=1000},Block[{t=Sort[RandomReal[1,n]]},
TimeSeries[RandomVariate[WienerProcess[][t]],{t}]]]

Consider filtering out high frequency. The cutoff frequency is specified in meaning absolute units:

Labeled[TabView[Table[TraditionalForm[𝒻==f]->ListLinePlot[{ts1000,LowpassFilter[ts1000,Quantity[f,"per second"]]},PlotStyle->{Directive[Thick,LightGray],Darker[Red]}],{f,{2,5,10,25,50}}]],HoldForm[LowpassFilter[ts1000,Quantity[𝒻,"per second"]]],Top]

Time Series

Updated MovingMap

N.B.: incompatible [q/79931], but enabling change, allowing for proper support of irregular time series

Window size by default is interpreted as absolute window size.

data=TimeSeries[FinancialData["SBUX",{"2010","2015"}]]

Three months moving median (almost the same number of data-points):

med=MovingMap[Median,data,{Quantity[3,"Months"],Center}]

Show[DateListLogPlot[data,PlotStyle->GrayLevel[.7]],DateListLogPlot[med,Joined->True,PlotStyle->Thick]]

Output time series inherits subset of time-stamps of parent, and these are irregular:

SubsetQ[data["Dates"],med["Dates"]]

{RegularlySampledQ[data],RegularlySampledQ[med]}

Three months moving median with regularly spaced bi-weekly windows, thus also thinning the data :

medr=MovingMap[Median,data,{{3,"Month"},Center,{2,"Week"}}]

RegularlySampledQ[medr]

Show[DateListLogPlot[data,PlotStyle->Directive[Thick,GrayLevel[.7]],PlotLegends->{Row[{"Original", HoldForm[data]},": "]}],
DateListLogPlot[{med,medr},Joined->True,PlotStyle->{Directive[Opacity[.5],Thickness[.016]],Automatic},PlotTheme->"VibrantColors",PlotLegends->{Row[{"Unthinned",HoldForm[med]},": "],Row[{"Thinned", HoldForm[medr]},": "]}],ImageSize->Large]

MovingMap

Access to time-stamps

MovingMap provides function with access to time-stamps of data-points within the window.

It also allows the function to return the time-stamp to associate with the value for a window.

MovingFunction[ Function[{times, vals}, …], ts ]

Consider well-known representation for the OrnsteinUhlenbeckProcess[μ,σ,θ,(Subscript[x, 0])]  as a time-changed  Wiener process.

x[t]==μ+((Subscript[x, 0])-μ)Exp[-θ t] + σ Exp[-θ t] w[ ((Exp[2 θ t]-1)/(2θ)) ]

t/.Solve[((Exp[2 θ t]-1)/(2θ))==τ∧θ>0∧τ>0,t,Reals]//Normal

Function takes values and times of Wiener process path data-points within a window:

ouPoint[vVec_,τVec_,{μ_,σ_,θ_,x0_}]:=Module[{tVec=(1/(2θ))Log[1+2 θ τVec],α=(((1+2 θ τVec))^(-1/2))},tVec->μ+(x0-μ)α+σ α vVec]

E.g., for Wiener process with value v==0.1 at τ==0.2 the function maps the OU process with value x==1.474 at time t==0.784

ouPoint[{0.2},{0.1},{1,2,3,1.2}]

wp=RandomFunction[WienerProcess[],{0,15,0.01},5];

Use relative sized window that contains one point

oup=MovingMap[ouPoint[#1,#2,{1,0.3,5.3,2}]&,wp,Quantity[1,"Event"]]

mf[t_]=Mean[OrnsteinUhlenbeckProcess[1,0.3,5.3,2][t]]

Show[Plot[mf[t],{t,0,oup["MaximumTime"]},PlotStyle->{Black,Thin}],ListPlot[oup],PlotRange->All]

Random matrices

Matrix Distributions

Guide page guide/MatrixDistributions

MatrixNormalDistribution

WishartMatrixDistribution

GaussianUnitaryMatrixDistribution

CircularRealMatrixDistribution

Matrix distributions support efficient sampling:

RandomVariate[GaussianOrthogonalMatrixDistribution[5000]];//AbsoluteTiming

Compare with direct sampling

(#+Transpose[#])/2&[RandomVariate[NormalDistribution[],{5000,5000}]];//AbsoluteTiming

Support for properties

Matrix distributions are interfaced to P&S framework through MatrixPropertyDistribution

logcondD=MatrixPropertyDistribution[ Log10[Max[Eigenvalues[mat]]/Min[Eigenvalues[mat]]],mat\[Distributed]WishartMatrixDistribution[23,{{2,1/2},{1/2,3}}]]

SmoothHistogram[RandomVariate[logcondD,10^5],Filling->Axis,PlotTheme->"Detailed"]

NExpectation[10^x,x\[Distributed]logcondD]

NProbability[x>0.6,x\[Distributed]logcondD]

The matrix property supported by MatrixPropertyDistribution can be scalar- or vector-valued

di=MatrixPropertyDistribution[ Diagonal[(u^\[Transpose]).m.u],{m\[Distributed]WishartMatrixDistribution[23,{{2,1/2},{1/2,3}}],u\[Distributed]CircularRealMatrixDistribution[2]}];

RandomVariate[di]

Histogram3D[RandomVariate[di,10^5],Automatic,PDF]

Tracy-Widom distribution

Externally requested universal distribution, related to the maximal eigenvalue of random matrices, and not easy to compute .

Plot[{PDF[TracyWidomDistribution[1],x],PDF[TracyWidomDistribution[2],x],PDF[TracyWidomDistribution[4],x]},{x,-5,3},Filling->Axis,PlotLegends->{β==1,β==2,β==4}]

Occurring in many, sometimes unexpected, branches of mathematics.

losLength[li_]:=Length[LongestOrderedSequence[li]]

n=50000;
list=Range[n];
lengths=Table[losLength[RandomSample[list]],{400}];

Lengths of the longest increasing subsequence in random permutations are well approximated by Tracy-Widom law:

QuantilePlot[(1/(n^(1/6)))(lengths-2.0(Sqrt[n])),TracyWidomDistribution[2]]

DistributionFitTest[(1/(n^(1/6)))(lengths-2.0(Sqrt[n])),TracyWidomDistribution[2]]

Noteworthy misc. additions/changes

Shifted Gompertz Distribution

ShiftedGompertzDistribution — life-time distribution used in marketing research

CDF[ShiftedGompertzDistribution[ξ,λ],t]

New multivariate normality hypothesis test

BaringhausHenzeTest — hypothesis test for multivariate normality, a.k.a. BHEP test

data=RandomVariate[ProductDistribution[BinormalDistribution[2/3],BinormalDistribution[-1/2]],121];

BaringhausHenzeTest[data,"TestDataTable"]

Automatically used in distribution fit test:

DistributionFitTest[data,Automatic,{"TestDataTable",All}]

Inhomogeneous Poisson Process

InhomogeneousPoissonProcess — Poisson process with variable intensity

ListStepPlot[λdata={{0,8}, {1,9},{2,7},{4,6},{5,12},{6,14},{7,11},{8,6},{9,4},{10,3},{11,8},{12,15},{13,12},{14,10},{15,8},{16,6},{17,12},{18,17},{19,15},{20,12},{21,6},{22,5},{23,7}},PlotLabel->"IPP intensity",ImageSize->Small]

λ[t_]=Interpolation[λdata,InterpolationOrder->0][t];
𝒫=InhomogeneousPoissonProcess[λ[t],t];

ListStepPlot[RandomFunction[𝒫,{0,23,1},4]]

Enhancements for ProbabilityDistribution

Automatic normalization

prd=ProbabilityDistribution[x z+y(1-x),{x,0,1},{y,0,1},{z,0,1},Method->"Normalize"]

Dedicated methods for polynomial, both 1D and multivariae

pts=RandomVariate[prd,2567];

Table[If[i==j,Histogram[pts[[All,i]],Automatic,PDF],QuantilePlot[pts[[All,i]],pts[[All,j]]]],{i,3},{j,3}]//Grid

Sampling methods

ProbabilityDistribution[x z+y(1-x),{x,0,1},{y,0,1},{z,0,1},Method->"Normalize"]

RandomVariate[prd,10,Method->"AcceptanceRejection"]

RLink supports R-3.2

Needs["RLink`"]

InstallR["RHomeLocation"->FileNameJoin[{"C:","Program Files","R","R-3.2.2"}]]

REvaluate["R.Version()"]

UninstallR[]

Spatial Analysis: Random Point

Uniform sampling from geometric regions

RandomPoint[Sphere[]]

%∈Sphere[]

Special regions

Graphics[Point[RandomPoint[StadiumShape[{{0,0},{2,3}},2],10^3]],ImageSize->Small]

RandomPoint[Simplex[IdentityMatrix[{11,10}]]]

Formula Regions

ℛ=ImplicitRegion[(((-1+x))^2)+(y^2)==1&&(x^2)+(y^2)+(z^2)==4,{x,y,z}];
pts=RandomPoint[ℛ,10^3];
Graphics3D[{PointSize[Tiny],Point[pts]},ImageSize->Medium]

ℛ=ParametricRegion[{(3+Cos[v])Cos[u],(3+Cos[v])Sin[u],Sin[v]},{{u,0,2Pi},{v,0,2Pi}}];
pts=RandomPoint[ℛ,4000,{{-4,4},{-4,4},{-1,1}}];
Graphics3D[{PointSize[Tiny],Point[pts]},Boxed->False,ImageSize->Medium]

Mesh Regions

2D

t=Text[Style["π",Bold]];

mr=DiscretizeGraphics[t,_Text,MaxCellMeasure->0.1];

pts=RandomPoint[mr,3000];

Graphics[{PointSize[Tiny],Point[pts]}]

3D

mm=Import[FileNameJoin[{NotebookDirectory[],"sasha.ply"}],{"PLY","MeshRegion"}];

pts=RandomPoint[m,3000];

{Graphics3D[{PointSize[Tiny],Point[pts]},ImageSize->Medium],MeshRegion[mm,PlotTheme->"SmoothShading",ImageSize->Medium]}

Spatial Analysis

Kriging

Covered in the talk by Eduardo Serna and Hsien-Ching Kao, today at 4pm

Spatial median

Spatial median:.

SpatialMedian[{(Subscript[p, 1]),…,(Subscript[p, m])}]==ArgMin[ (∑_(k=1)^m)d[p,(Subscript[p, k])],p]

cities=CityData[{All,"Montana"}];
locations=CityData[#,"Coordinates"]&/@cities;

unweighted=GeoPosition[SpatialMedian[locations]]

populations=QuantityMagnitude[CityData[#,"Population"]&/@cities];

weighted=GeoPosition[SpatialMedian[locations,populations]]

Legended[GeoGraphics[{Gray,PointSize[0.006],Point/@(GeoPosition/@locations),Red,PointSize[0.02],Point[unweighted],Black,Point[weighted],EdgeForm[Black],FaceForm[Red],Polygon[LinguisticAssistant]},ImageSize->Large,Frame->True],PointLegend[{Gray,Red,Black},{"Cities","Unweighted spatial median","Weighted spatial median"}]]

Quantities

Sources

Increasing number of built-in data collections give results using quantities

MortalityData[\[LeftAssociation]"Age"->(DateObject["Sun 4 Jul 1976", DateObject[{1976,7,4}]]),"Gender"->"Male"\[RightAssociation],"LifeExpectancy"]

HumanGrowthData[Association["Age"->Quantity[8,"Years"],"Gender"->(Entity["male", Entity["Gender","Male"], "Entity[\"Gender\", \"Male\"]", "gender"])],"Weight",Quantity[90,"Percent"]]

WeatherData[{"Champaign","Illinois","UnitedStates"},"Temperature"]

QuantityArray

qa = QuantityArray[RandomVariate[GammaDistribution[120,1],100],"Meters"]

Mean[%]

Kurtosis[qa]

QuantilePlot[qa,FrameLabel->Automatic]

Distributions

Determine parameters as quantities:

{fitBoy}=Solve[Quantile[NormalDistribution[μ,σ],{0.25,0.75}]==(HumanGrowthData[Association["Age"->LinguisticAssistant,"Gender"->(Entity["male", Entity["Gender","Male"], "Entity[\"Gender\", \"Male\"]", "gender"])],"Weight",#]&/@{LinguisticAssistant,LinguisticAssistant}),{μ,σ}]

{fitGirl}=Solve[Quantile[NormalDistribution[μ,σ],{0.25,0.75}]==(HumanGrowthData[Association["Age"->LinguisticAssistant,"Gender"->(Entity["female", Entity["Gender","Female"], "Entity[\"Gender\", \"Female\"]", "gender"])],"Weight",#]&/@{LinguisticAssistant,LinguisticAssistant}),{μ,σ}]

Use these in distributions:

wdBoy=NormalDistribution[μ,σ]/.fitBoy

wdGirl=NormalDistribution[μ,σ]/.fitGirl

Form a derived distribution, assuming independence of distributions:

TransformedDistribution[x-y,{x\[Distributed]wdBoy,y\[Distributed]wdGirl}]

Sample:

data=RandomVariate[wdBoy,1000]

Estimate:

EstimatedDistribution[UnitConvert[data,"Kilograms"],NormalDistribution[m,s]]

UnitConvert[wdBoy,"Kilograms"]

Compute with

Expectation[w\[Conditioned]w>Quantity[26,"Kilograms"],w\[Distributed]wdBoy]

NProbability[b>g,{b\[Distributed]wdBoy,g\[Distributed]wdGirl}]

UnitConvert[%,"Percent"]

Expectation[x,Quantity[x,"Ounces"]\[Distributed]wdBoy]

Expectation[QuantityMagnitude[x,"Ounces"],x\[Distributed]wdBoy]

Find quantiles

Quartiles[wdBoy]

(HumanGrowthData[Association["Age"->LinguisticAssistant,"Gender"->(Entity["male", Entity["Gender","Male"], "Entity[\"Gender\", \"Male\"]", "gender"])],"Weight",#]&/@{LinguisticAssistant,LinguisticAssistant,LinguisticAssistant})

Labeled[Legended[Framed@ImageCollage[{,,,,},ImagePadding->1,Background->GrayLevel[.9]],Grid[Transpose[{Thumbnail[#,30]&/@{,,,,},Style[#,"Times"]&/@{"Gosia Konwerska","Devendra Kapadia","Hsien-Ching Kao","Eduardo Serna","Oleksandr Pavlyk"}}],Alignment->Left]],Style["Probability & Statistics Team","Section"],{Top,Top}]