Quantities in Probability and Statistics

Gosia Konwerska
 
Wolfram Research
Champaign, IL - October 19, 2016

Join the Conversation #WolframTechConf

The Plan

1. Quantity Data & Descriptive Measures

2. QuantityDistribution

3. Solvers

4. Derived Distributions

5. Estimation

6. Fun Part

Quantity Data & Descriptive Measures (1/3)

Stars

aux=(EntityClass["100 nearest stars", EntityClass["Star","StarNearest100"], "EntityClass[\"Star\", \"StarNearest100\"]", "stars"])[{"Name","DistanceFromSun"}];
data=DeleteMissing[aux,1,2];

data[[1]]

Extract distances.

distances=data[[All,2]];

Plot histogram of distances.

Histogram[distances,AxesLabel->Automatic]

Compute basic statistics.

stats={Min,Mean,Median,Max};

TableForm[Through[stats[distances]],TableHeadings->{stats}]

The distance to the 3rd nearest star.

RankedMin[distances,3]

Quantity Data & Descriptive Measures (2/3)

Data with mixed units

data={Quantity[2.3,"Meters"],Quantity[12,"Feet"],Quantity[2.4,"Yards"],Quantity[342.3,"Centimeters"]};

Check if the units are compatible.

CompatibleUnitQ[data]

Compute means.

Mean[data]

GeometricMean[data]

Convert data to another unit.

UnitConvert[data,"Furlong"]

Convert data to standard SI unit.

UnitConvert[data]

Quantity Data & Descriptive Measures (3/3)

QuantityArray

Store vectors with units using QuantityArray.

sample=RandomVariate[BinormalDistribution[.2],10]

qa=QuantityArray[sample,{"Meters","Pounds"}]

Compute basic descriptive measures.

stats={Mean,Median,StandardDeviation};

TableForm[Through[stats[qa]],TableHeadings->{stats}]

QuantityDistribution (1/2)

Univariate Model

Create the QuantityDistibution of arrival times in seconds from the distribution of the magnitudes.

𝒟=QuantityDistribution[ExponentialDistribution[1/3],"Seconds"]

Use QuantityMagnitude and QuantityUnit to extract the distribution of the magnitude of the associated random variable and its unit.

QuantityMagnitude[𝒟]

QuantityUnit[𝒟]

Use UnitConvert to obtain quantity distribution in new units.

UnitConvert[𝒟,"Minutes"]

Simulate quantities. Single simulation is a Quantity.

RandomVariate[𝒟]

Vector of simulated variates is a QuantityArray.

sample=RandomVariate[𝒟,10^3]

Histogram[sample,Automatic,"PDF",AxesLabel->Automatic]

Compute CDF at a point.

CDF[𝒟,Quantity[.3,"Minutes"]]

Mean and variance.

Mean[𝒟]

Variance[𝒟]

QuantityDistribution (2/2)

Random Vectors with Different Units

Define a joint height-weight distribution with assumed correlation of 0.65. Use quantities as parameters for the distribution.

hw𝒟=BinormalDistribution[{(Quantity[1.8, "m", meters, "Meters"]),(Quantity[85, "kg", kilograms, "Kilograms"])},{(Quantity[0.15, "m", meters, "Meters"]),(Quantity[9, "kg", kilograms, "Kilograms"])},0.65]

Mean and variance.

Mean[hw𝒟]

Variance[hw𝒟]

Moments.

Moment[hw𝒟,{2,3}]

Cumulant[hw𝒟,{1,1}]

Simulate single observation.

RandomVariate[hw𝒟]

Simulate multiple observations.

data=RandomVariate[hw𝒟,10^4]

Probability density histogram.

Histogram3D[data,40,PDF]

Distribution of height.

MarginalDistribution[hw𝒟,1]

Solvers

Probability & Expectation

QuantityDistribution and Quantity objects are accepted in framework solvers like Probability and Expectation.

Height-Weight Distribution

The joint height-weight distribution with assumed correlation of 0.65.

hw𝒟=BinormalDistribution[{(Quantity[1.8, "m", meters, "Meters"]),(Quantity[85, "kg", kilograms, "Kilograms"])},{(Quantity[0.15, "m", meters, "Meters"]),(Quantity[9, "kg", kilograms, "Kilograms"])},0.65]

w𝒟=MarginalDistribution[hw𝒟,2]

Find the expected weight in pounds.

Expectation[UnitConvert[w,"Pounds"],w\[Distributed]w𝒟]

Use option TargetUnits to find the expected weight in pounds.

Expectation[w,w\[Distributed]w𝒟,TargetUnits->"Pounds"]

Find the expected absolute deviation from the mean height.

Expectation[Abs[h-First[Mean[hw𝒟]]],{h,w}\[Distributed]hw𝒟]

Compute the conditional probability of a person’s height being over 1.75 meters, given that the person’s weight is above 72 kilograms.

NProbability[h>Quantity[1.75,"m"]\[Conditioned]w>Quantity[72,"kg"],{h,w}\[Distributed]hw𝒟]

Apple Production

Poland’s rate of apple production:

rate = EntityValue[Entity["Country","Poland"],EntityProperty["Country","AgricultureProduction",{"AgricultureProduct" ->"Apple"}]]//N

Assumed distribution of actual production:

dist = RayleighDistribution[rate]

Plot[PDF[dist,Quantity[x,QuantityUnit[dist]]],{x,0,10^7}]

Compute probability that daily apple production exceeds weight of 15 International Space Stations:

const=(Entity["International Space Station", Entity["Satellite","25544"], "Entity[\"Satellite\", \"25544\"]", "satellite"])[EntityProperty["Satellite", "DryMass"]];

Probability[ Quantity[1, "Day"]*r >15* const, Distributed[r, dist]]

Derived Distributions

Derived distributions allow QuantityDistribution input.

CensoredDistribution
MarginalDistribution
MixtureDistribution
OrderDistribution
ParameterMixtureDistribution
ProductDistribution
QuantityDistribution
TransformedDistribution
TruncatedDistribution

TransformedDistribution

The joint height-weight distribution with assumed correlation of 0.65.

hw𝒟=BinormalDistribution[{(Quantity[1.8, "m", meters, "Meters"]),(Quantity[85, "kg", kilograms, "Kilograms"])},{(Quantity[0.15, "m", meters, "Meters"]),(Quantity[9, "kg", kilograms, "Kilograms"])},0.65]

The distribution of the body mass index (BMI) for the assumed population can be found using TransformedDistribution.

bmi𝒟=TransformedDistribution[w/h^2,{h,w}\[Distributed]hw𝒟]

Compute numerically the mean value of the BMI distribution.

NExpectation[x,x\[Distributed]bmi𝒟]

Estimate the probability density function of the BMI distribution using sample histogram.

sample=RandomVariate[bmi𝒟,10^5];
Histogram[sample,Automatic,"PDF",AxesLabel->Automatic]

Estimation (1/2)

Parametric Models

The time series of the temperatures in Chicago in 2015.

(*TimeSeries[WeatherData["Chicago","Temperature",{2015}],MissingDataMethod->"Interpolation"]*)

tempsChicago=TimeSeries[];

hist=Histogram[tempsChicago,30,"PDF",AxesLabel->Automatic]

Extract the values.

temps=Values[tempsChicago]

Find parametric estimate for the temperature distribution. The resulting model takes on the units of the data.

ed=EstimatedDistribution[temps,SkewNormalDistribution[a,b,c]]

Specify the units of the estimated model.

edF=EstimatedDistribution[temps,QuantityDistribution[SkewNormalDistribution[a,b,c],"DegreesFahrenheit"]]

Check that both models are equivalent.

UnitConvert[ed,"DegreesFahrenheit"]

Estimation (2/2)

Nonparametric Estimators

Temperature

Nonparametric estimate for the temperature distribution. The data of Chicago temperatures from 2015.

tempsChicago

temps=Values[tempsChicago]

Compute nonparametric estimated distribution using SmoothKernelDistribution.

skd=SmoothKernelDistribution[temps]

Compare the parametric and nonparametric estimates with the data histogram.

Show[hist,Plot[{PDF[ed,Quantity[x,"DegreesCelsius"]],PDF[skd,Quantity[x,"DegreesCelsius"]]},{x,-20,42},PlotLegends->{ed,skd}]]

Wind

The time series of wind speeds from Chicago, IL, over the period of two years.

(*TimeSeries[WeatherData["Chicago","WindSpeed",{DateObject[{2014,1,1}],DateObject[{2015,12,31}]}],MissingDataMethod->Automatic]*)

windChicago=TimeSeries[];
winds=windChicago["Values"]

Use a nonparametric model of the turbine's power output as a function of wind speed to estimate the average power output for a GE 1.5 MW wind turbine installed at the location. Any numeric bounds for SmoothKernelDistribion must be provided with appropriate units.

ws𝒟=SmoothKernelDistribution[winds,Automatic,{"Bounded",(Quantity[0, "km""/""h", kilometers per hour, (("Kilometers")/("Hours"))]),"Gaussian"}]

turbine=Interpolation[QuantityArray[{{0.`,0.`},{0.5`,0.`},{1.`,0.`},{1.5`,0.`},{2.`,0.`},{2.5`,0.`},{3.`,0.`},{3.5`,0.`},{4.`,36.`},{4.5`,66.`},{5.`,104.`},{5.5`,150.`},{6.`,205.`},{6.5`,269.`},{7.`,344.`},{7.5`,428.`},{8.`,528.`},{8.5`,644.`},{9.`,774.`},{9.5`,926.5`},{10.`,1079.`},{10.5`,1211.`},{11.`,1342.`},{11.5`,1401.`},{12.`,1460.`},{12.5`,1477.`},{13.`,1494.`},{13.5`,1500.`},{30.`,1500.`}},{"Meters"/"Seconds","Kilowatts"}]//Normal,InterpolationOrder->1];

NExpectation[turbine[v],v\[Distributed]ws𝒟]

An Interesting Example To Rule Them All

Windchill Index

The temperatures in Chicago during recent winters.

(*TemporalData@Table[TimeSeries[WeatherData["Chicago","Temperature",{{year,12,1},{year+1,3,31}}],MissingDataMethod->"Interpolation"],{year,2000,2015}]*)

winterTemps=TemporalData[];

DateListPlot[winterTemps,FrameLabel->Automatic]

Extract the values and join into single QuantityArray.

temps=Join@@winterTemps["ValueList"]

Histogram[temps,AxesLabel->Automatic]

The wind speeds in Chicago during recent winters.

(*TemporalData@Table[TimeSeries[WeatherData["Chicago","WindSpeed",{{year,12,1},{year+1,3,31}}],MissingDataMethod->Automatic],{year,2000,2015}]*)

winterWinds=TemporalData[];

DateListPlot[winterWinds,FrameLabel->Automatic]

Extract the values and join into single QuantityArray.

winds=Join@@winterWinds["ValueList"]

Histogram[winds,AxesLabel->Automatic]

Nonparametric estimate of winter temperatures distribution.

temp𝒟=SmoothKernelDistribution[temps]

Parametric estimate of winter wind speeds distribution.

wind𝒟=EstimatedDistribution[winds, WeibullDistribution[a,b,-1]]

Histograms and PDFs of fitted distributions.

{Show[Histogram[temps,Automatic,PDF,AxesLabel->Automatic],Plot[PDF[temp𝒟,Quantity[x,"DegreesCelsius"]],{x,-40,40}],PlotLabel->"Temperatures",ImageSize->Medium],Show[Histogram[winds,Automatic,PDF,AxesLabel->Automatic],Plot[PDF[wind𝒟,Quantity[x,"kmph"]],{x,0,70}],PlotLabel->"Winds",ImageSize->Medium]}

Distribution of the windchill calculated from the temperature and wind speed distributions. The formula assumes quantities come in imperial units so the unit conversions are required.

windchill=TransformedDistribution[35.74+0.6215 QuantityMagnitude[Ta,"DegreesFahrenheit"]+(0.4275 QuantityMagnitude[Ta,"DegreesFahrenheit"] -35.75) ((QuantityMagnitude[V,"Miles"/"Hours"])^(+0.16)),{V\[Distributed]wind𝒟,Ta\[Distributed]temp𝒟}]

DistributionParameterQ[windchill]

Plot smooth kernel density estimate of the wind chill index in Chicago from the population sample.

SmoothHistogram[RandomVariate[windchill,10^5],Filling->Axis,AxesLabel->{"°F"},Epilog->{Directive[Dashed,Thick,Red],InfiniteLine[{{32,0},{32,1}}]}]

Thank You!

Questions?

gosiak@wolfram.com

Related Talks

Date Visualization

Nirmal Malapaka, Wednesday Oct 19, 5:30pm, Spartan/Golden Gopher

Time Series Processing

Gosia Konwerska, Thursday Oct 20, 1:00pm, Spartan/Golden Gopher

What’s New in Probability and Statistics

Hsien-Ching Kao, Thursday Oct 20, 8:30am, Spartan/Golden Gopher