{ "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 }