Overview

The goal of this package is to provide an intuitive and mathematically sound interface for taking expectations of random variables and their higher-order functions (i.e., if $X \sim N(0, 1)$, what is $\mathbb{E}[\sin(X)]$?).

The underlying random variables are assumed to be distributions from Distributions.jl. Currently, only univariate distributions are supported.

Installation

To install, run:

julia> using Pkg

julia> pkg"add Expectations Distributions"
┌ Warning: The Pkg REPL mode is intended for interactive use only, and should not be used from scripts. It is recommended to use the functional API instead.
└ @ Pkg.REPLMode /usr/local/julia/share/julia/stdlib/v1.9/Pkg/src/REPLMode/REPLMode.jl:382
   Resolving package versions...
┌ Warning: The active manifest file at `/juliateam/.julia/packages/Expectations/uN97T/docs/Manifest.toml` has an old format that is being maintained.
│ To update to the new format, which is supported by Julia versions ≥ 1.6.2, run `import Pkg; Pkg.upgrade_manifest()` which will upgrade the format without re-resolving.
│ To then record the julia version re-resolve with `Pkg.resolve()` and if there are resolve conflicts consider `Pkg.update()`.
└ @ Pkg.Types /usr/local/julia/share/julia/stdlib/v1.9/Pkg/src/manifest.jl:316
    Updating `~/.julia/packages/Expectations/uN97T/docs/Project.toml`
  [31c24e10] + Distributions v0.25.104
  [2fe49d83] ~ Expectations v1.9.2 `~/.julia/packages/Expectations/uN97T` ⇒ v1.9.2
    Updating `~/.julia/packages/Expectations/uN97T/docs/Manifest.toml`
  [2fe49d83] ~ Expectations v1.9.2 `~/.julia/packages/Expectations/uN97T` ⇒ v1.9.2

julia> using Expectations, Distributions

The Expectation Operator

The key object in this package is an expectation operator, or an object <: Expectation. These include all objects capable of being called on a function; e.g. that support a method function (e::Expectation)(f::Function). You can create these as following:

julia> dist = Normal();

julia> E = expectation(dist)
IterableExpectation{Vector{Float64}, Vector{Float64}}([-9.706235997359522, -8.680837722732212, -7.825051744352812, -7.055396866960293, -6.339997686869597, -5.662381850082874, -5.012600596486518, -4.384020365898052, -3.7718944231592366, -3.1726346394204032  …  3.1726346394204032, 3.7718944231592366, 4.384020365898052, 5.012600596486518, 5.662381850082874, 6.339997686869597, 7.055396866960293, 7.825051744352812, 8.680837722732212, 9.706235997359522], [1.6408069991319633e-21, 1.5855609449663175e-17, 1.6240801299724176e-14, 4.5734258713261955e-12, 5.178459467189773e-10, 2.88217515404765e-8, 8.90908886862124e-7, 1.6579981630673396e-5, 0.00019651294398482591, 0.001544707339866085  …  0.001544707339866085, 0.00019651294398482591, 1.6579981630673396e-5, 8.90908886862124e-7, 2.88217515404765e-8, 5.178459467189773e-10, 4.5734258713261955e-12, 1.6240801299724176e-14, 1.5855609449663175e-17, 1.6408069991319633e-21])

You can also choose algorithms and default parameters (see below for list):

julia> E = expectation(dist, Gaussian; n = 30) # Could have done expectation(dist) or expectation(dist; n = 30)
IterableExpectation{Vector{Float64}, Vector{Float64}}([-9.706235997359522, -8.680837722732212, -7.825051744352812, -7.055396866960293, -6.339997686869597, -5.662381850082874, -5.012600596486518, -4.384020365898052, -3.7718944231592366, -3.1726346394204032  …  3.1726346394204032, 3.7718944231592366, 4.384020365898052, 5.012600596486518, 5.662381850082874, 6.339997686869597, 7.055396866960293, 7.825051744352812, 8.680837722732212, 9.706235997359522], [1.6408069991319633e-21, 1.5855609449663175e-17, 1.6240801299724176e-14, 4.5734258713261955e-12, 5.178459467189773e-10, 2.88217515404765e-8, 8.90908886862124e-7, 1.6579981630673396e-5, 0.00019651294398482591, 0.001544707339866085  …  0.001544707339866085, 0.00019651294398482591, 1.6579981630673396e-5, 8.90908886862124e-7, 2.88217515404765e-8, 5.178459467189773e-10, 4.5734258713261955e-12, 1.6240801299724176e-14, 1.5855609449663175e-17, 1.6408069991319633e-21])

These objects can then be applied to functions:

julia> E(x -> x)
-5.028942338420413e-18

julia> E(x -> x^2)
0.9999999999999984

There is also a convenience function to evaluate expectations directly, without returning the operator:

julia> f = x -> x^2
#5 (generic function with 1 method)

julia> expectation(f, dist)
0.9999999999999984

In general, expectation(f, dist, ...) is equivalent to E(f), where E = expectation(dist, ...).

IterableExpectation

The only concrete subtype of Expectation currently supported is IterableExpectation{NT, WT}. These are expectations for which we have a discrete vector of quadrature nodes and weights, either defined by user fiat, or set algorithmically. These support some additional behavior:

julia> nodeList = nodes(E);

julia> vals = map(x -> x^2, nodeList);

julia> E * vals
0.9999999999999984

julia> (2E) * vals
1.999999999999997

The above behavior, in some sense, puts the "operator" in "expectation operator"; that is, it allows it to move elements of a vector space around, and to be scalar-multiplied.

User-Defined Nodes

There are some situations where we are forced to use a specific set of nodes. In those situations, E = expectation(dist, nodes) will create the relevant object.

Mixture Models

We also have support for univariate mixture models. The MixtureExpectation type is a struct with two fields:

  • expectations, which is a list of IterableExpectation objects.

  • mixtureweights, which is the mixing probabilities over the various model components.

The mixture models are constructed using the default settings for each component distribution. (It still accepts kwargs, which are applied to each.)

julia> d = MixtureModel([Uniform(), Normal(), Gamma()]);

julia> E = expectation(d; n = 30); # n = 30 nodes for each

julia> @show typeof(E)
typeof(E) = MixtureExpectation{Vector{IterableExpectation{Vector{Float64}, Vector{Float64}}}, Vector{Float64}}
MixtureExpectation{Vector{IterableExpectation{Vector{Float64}, Vector{Float64}}}, Vector{Float64}}

julia> E(x -> abs(x))
0.769651856398927

If you want to change this, you should construct each distribution separately, and then chain them together.

julia> E1 = expectation(Uniform())
IterableExpectation{Vector{Float64}, Vector{Float64}}([0.001553257962675192, 0.008165938360126412, 0.019989067515846226, 0.036899976285362845, 0.05871973210397369, 0.08521711880861582, 0.11611128394758691, 0.15107475260334208, 0.18973690850537855, 0.23168792592899007  …  0.7683120740710099, 0.8102630914946214, 0.848925247396658, 0.8838887160524131, 0.9147828811913842, 0.9412802678960264, 0.9631000237146372, 0.9800109324841537, 0.9918340616398735, 0.9984467420373249], [0.003984096248083297, 0.009233234155545493, 0.01439235394166165, 0.019399596284813594, 0.024201336415297054, 0.02874657810880951, 0.03298711494109026, 0.036877987368852574, 0.040377947614710114, 0.04344989360054148  …  0.04344989360054148, 0.040377947614710114, 0.036877987368852574, 0.03298711494109026, 0.02874657810880951, 0.024201336415297054, 0.019399596284813594, 0.01439235394166165, 0.009233234155545493, 0.003984096248083297])

julia> E2 = expectation(Normal())
IterableExpectation{Vector{Float64}, Vector{Float64}}([-9.706235997359522, -8.680837722732212, -7.825051744352812, -7.055396866960293, -6.339997686869597, -5.662381850082874, -5.012600596486518, -4.384020365898052, -3.7718944231592366, -3.1726346394204032  …  3.1726346394204032, 3.7718944231592366, 4.384020365898052, 5.012600596486518, 5.662381850082874, 6.339997686869597, 7.055396866960293, 7.825051744352812, 8.680837722732212, 9.706235997359522], [1.6408069991319633e-21, 1.5855609449663175e-17, 1.6240801299724176e-14, 4.5734258713261955e-12, 5.178459467189773e-10, 2.88217515404765e-8, 8.90908886862124e-7, 1.6579981630673396e-5, 0.00019651294398482591, 0.001544707339866085  …  0.001544707339866085, 0.00019651294398482591, 1.6579981630673396e-5, 8.90908886862124e-7, 2.88217515404765e-8, 5.178459467189773e-10, 4.5734258713261955e-12, 1.6240801299724176e-14, 1.5855609449663175e-17, 1.6408069991319633e-21])

julia> E3 = expectation(Gamma())
IterableExpectation{Vector{Float64}, Vector{Float64}}([0.04448936583326675, 0.23452610951961822, 0.5768846293018851, 1.0724487538178178, 1.7224087764446459, 2.5283367064257956, 3.4922132730219944, 4.616456769749767, 5.903958504174245, 7.358126733186242  …  44.50920799575494, 49.22439498730864, 54.33372133339691, 59.89250916213402, 65.97537728793506, 72.68762809066271, 80.18744697791352, 88.7353404178924, 98.82954286828398, 111.7513980979377], [0.10921834195236964, 0.21044310793883128, 0.23521322966984406, 0.1959033359728804, 0.12998378628607243, 0.07057862386571843, 0.03176091250917527, 0.011918214834838561, 0.0037388162946115316, 0.0009808033066149562  …  2.1197922901636102e-19, 2.054429673788039e-21, 1.3469825866373934e-23, 5.661294130397274e-26, 1.418560545462764e-28, 1.9133754944542112e-31, 1.1922487600982224e-34, 2.6715112192401224e-38, 1.3386169421062427e-42, 4.5105361938989814e-48])

julia> E = MixtureExpectation([E1, E2, E3], [1/3, 1/3, 1/3])
MixtureExpectation{Vector{IterableExpectation{Vector{Float64}, Vector{Float64}}}, Vector{Float64}}(IterableExpectation{Vector{Float64}, Vector{Float64}}[IterableExpectation{Vector{Float64}, Vector{Float64}}([0.001553257962675192, 0.008165938360126412, 0.019989067515846226, 0.036899976285362845, 0.05871973210397369, 0.08521711880861582, 0.11611128394758691, 0.15107475260334208, 0.18973690850537855, 0.23168792592899007  …  0.7683120740710099, 0.8102630914946214, 0.848925247396658, 0.8838887160524131, 0.9147828811913842, 0.9412802678960264, 0.9631000237146372, 0.9800109324841537, 0.9918340616398735, 0.9984467420373249], [0.003984096248083297, 0.009233234155545493, 0.01439235394166165, 0.019399596284813594, 0.024201336415297054, 0.02874657810880951, 0.03298711494109026, 0.036877987368852574, 0.040377947614710114, 0.04344989360054148  …  0.04344989360054148, 0.040377947614710114, 0.036877987368852574, 0.03298711494109026, 0.02874657810880951, 0.024201336415297054, 0.019399596284813594, 0.01439235394166165, 0.009233234155545493, 0.003984096248083297]), IterableExpectation{Vector{Float64}, Vector{Float64}}([-9.706235997359522, -8.680837722732212, -7.825051744352812, -7.055396866960293, -6.339997686869597, -5.662381850082874, -5.012600596486518, -4.384020365898052, -3.7718944231592366, -3.1726346394204032  …  3.1726346394204032, 3.7718944231592366, 4.384020365898052, 5.012600596486518, 5.662381850082874, 6.339997686869597, 7.055396866960293, 7.825051744352812, 8.680837722732212, 9.706235997359522], [1.6408069991319633e-21, 1.5855609449663175e-17, 1.6240801299724176e-14, 4.5734258713261955e-12, 5.178459467189773e-10, 2.88217515404765e-8, 8.90908886862124e-7, 1.6579981630673396e-5, 0.00019651294398482591, 0.001544707339866085  …  0.001544707339866085, 0.00019651294398482591, 1.6579981630673396e-5, 8.90908886862124e-7, 2.88217515404765e-8, 5.178459467189773e-10, 4.5734258713261955e-12, 1.6240801299724176e-14, 1.5855609449663175e-17, 1.6408069991319633e-21]), IterableExpectation{Vector{Float64}, Vector{Float64}}([0.04448936583326675, 0.23452610951961822, 0.5768846293018851, 1.0724487538178178, 1.7224087764446459, 2.5283367064257956, 3.4922132730219944, 4.616456769749767, 5.903958504174245, 7.358126733186242  …  44.50920799575494, 49.22439498730864, 54.33372133339691, 59.89250916213402, 65.97537728793506, 72.68762809066271, 80.18744697791352, 88.7353404178924, 98.82954286828398, 111.7513980979377], [0.10921834195236964, 0.21044310793883128, 0.23521322966984406, 0.1959033359728804, 0.12998378628607243, 0.07057862386571843, 0.03176091250917527, 0.011918214834838561, 0.0037388162946115316, 0.0009808033066149562  …  2.1197922901636102e-19, 2.054429673788039e-21, 1.3469825866373934e-23, 5.661294130397274e-26, 1.418560545462764e-28, 1.9133754944542112e-31, 1.1922487600982224e-34, 2.6715112192401224e-38, 1.3386169421062427e-42, 4.5105361938989814e-48])], [0.3333333333333333, 0.3333333333333333, 0.3333333333333333])

julia> E(x -> abs(x))
0.7696518563989267

Supported Distributions, Algorithms, Keywords, and Defaults

Here is a list of currently supported distributions, along with keyword arguments and their defaults.

Distribution NameAlgorithm (Julia Type)Keywords and DefaultsRestrictions
Discrete UnivariateFiniteDiscrete <: QuadratureAlgorithmN/ASupport must be finite.
Continuous UnivariateGauss-Legendre (Gaussian <: QuadratureAlgorithm)n = 500Support must be a compact interval $[a, b]$.
Continuous UnivariateQNWDist[1] (QuantileRange <: ...)n = 50, q0 = 0.001, qN = 0.999Distribution must be nondegenerate.
Normal <: Continuous UnivariateGauss-Hermite (...)n = 30...
LogNormal <: ...Gauss-Hermite (...)n = 30...
Beta <: ...Gauss-Jacobi (...)n = 32...
ChiSq <: ...Gauss-Laguerre (...)n = 32...
Uniform <: ...Gauss-Legendre (...)n = 30...
Exponential <: ...Gauss-Laguerre (...)n = 32...
Gamma <: ...Gauss-Laguerre (...)n = 32...
UnivariateTrapezoidal <: ExplicitQuadratureAlgorithmN/AAll nodes must be inside distribution's support.

Some unbounded distributions are currently not supported (e.g., Poisson). Depending on your use case truncating may be a feasible option:

julia> E = expectation(Pareto()) # Throws error
ERROR: UndefVarError: `Pareto` not defined

julia> E = expectation(truncated(Pareto(),0.0,1000.0)) # Truncated Pareto on [0,1000]
ERROR: UndefVarError: `Pareto` not defined

See Distributions.truncated for more. Of course, truncating the distribution also affects its properties.

Mathematical Details and References

The specific quadrature algorithms come from the FastGaussQuadrature.jl library, which is maintained by Alex Townsend of Cornell University. Much of the quadrature code came from the DistQuads.jl library, which is maintained by Patrick K. Mogensen at the University of Copenhagen. In addition, there are some objects contributed by individual users; see docstring for citations.

WARNING: It is important to be aware of the deficiencies of numerical quadrature schemes. For example, it is recommended to be careful when using these methods for the following classes of functions and situations:

  • Discontinuous or nondifferentiable functions (even if the function is a.e.-differentiable)
  • Periodic/oscillatory functions with a high frequency
  • Extremely large numbers of quadrature nodes, which may lead to vanishingly small weights.

Contact

If you would like to get in touch, please do one of the following:

  • Issue requests: Open an issue on the package repository with the tag feature request.
  • Bugs: Same as above, but with the tag bug.
  • Pull Request: We are always open to new functionality. If you have a feature you'd like to add (say, a new distribution or algorithm), once you prepare a PR with the feature and some tests, open it in the usual way.
  • Other: You can reach out to Arnav Sood at misc@arnavsood.com
  • Citation: If this package was helpful in your research work, you may consider citing the package in whatever method is appropriate for your field.
  • 1This is a quadrature scheme written by Spencer Lyon (PhD. NYU) as part of the QuantEcon project. Used with permission.