{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Translation of the first example in the talk \"[Emily Gorcenski - Polynomial Chaos: A technique for modeling uncertainty](https://www.youtube.com/watch?v=Z-Qio-n6yPc)\". The idea is to approximate a general, potentially empirical, probability distribution using a suitable polynomial expansion (Polynomial/Wiener Chaos). With [Julia](https://github.com/JuliaLang/julia) and [PolyChaos](https://timueh.github.io/PolyChaos.jl/stable) package it simplified beautifully. For full explanation see related [blog post](https://john-hearn.info/articles/polynomial-chaos)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"#]add PolyChaos, Distributions, Plots"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"using PolyChaos, Distributions, Plots"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First define our Hermite polynomials, $H_i(x)$\n",
"(called Gauss in PolyChaos because Hermite is reserved for the physicists Hermite polynomials which are slightly different)."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"op_gauss = GaussOrthoPoly(20)\n",
"H = [x -> evaluate(i, x, op_gauss) for i in 0:5]\n",
"plot(-1.5:0.01:1.5, H, leg=false)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"#7 (generic function with 1 method)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"inv_cdf(dist) = u -> quantile(dist, u)\n",
"\n",
"h = inv_cdf(Exponential())\n",
"l = inv_cdf(Normal())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"PolyChaos calculates our inner (scalar) product for us, conveniently. \n",
"\n",
"The integrand is just the Galerkin projection with the transformed functions, $h(u)$ and $l(u)$."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"integrand (generic function with 1 method)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sp = computeSP2(op_gauss)\n",
"integrand(i) = u -> h(u)*H[i](l(u))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Perform the integration up to $p$ polynomials, we can use PolyChaos for this too although any integration scheme works."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Float64,1}:\n",
" 0.9999993692484951\n",
" 0.9031921803407967\n",
" 0.29780023616476886\n",
" 0.03334230632910607\n",
" -0.0024169819236717414"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p = 5\n",
"int_op = Uniform01OrthoPoly(1000, addQuadrature=true)\n",
"ki = [integrate(integrand(i), int_op) / sp[i] for i in 1:p]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we have the $k_i$s to do the transform we can reconstitute the approximated distribution by adding up a load of transformed Gaussian random variables, $\\zeta_i$."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ζ = randn(5000)\n",
"Σ = sum(ki[i] * H[i](ζ) for i in 1:p)\n",
"histogram(Σ, normed=true, leg=false)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"kernelspec": {
"display_name": "Julia 1.4.1",
"language": "julia",
"name": "julia-1.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.4.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}