In this example, we do some basic orthogonal expansions of 1-D functions of random variables and then use sampling from the approximation and the true function to evaluate the expansion. In future work we will show how to compute the expansions using quadrature. The source for the markdown is here.
Polynomial Chaos Expansions
Ryan G. McClarren
November 4, 2014
Given the function \(f(x) = e^{-x}\) where \(X \sim N(0,1)\) we can expand \[f(x) \approx \sum_{n=0}^4 c_n \psi_n(x),\] where the \(\psi_n(x)\) are Hermite polynomials which have the weight function \[w(x) = \frac{1}{\sqrt{2 \pi}} e^{-x^2/2}.\] We will use a third-order expansion for demonstration purposes: \[f(x) \approx c_0 + c_1 x + c_2 \frac{1}{\sqrt{2}}(x^2-1) + c_3 \frac{1}{\sqrt{6}}(x^3-3x).\]
We can show that the values of the expansion constants lead to \[f(x) \approx e^{1/2} \left( 1 – x + \frac{1}{2}(x^2-1) – \frac{1}{{6}}(x^3-3x)\right).\]
Now to generate samples from our expansion, we sample \(x\) and then compute \(f(x)\) using the approximation. Of course, in this case we can evaluate the full model, but in practice we would be running a simulation to propagate the sample through the model. With polynomial chaos, once we have an expansion the work is basically done.
Nsamples = 10000;
samples = data.frame(x = rnorm(Nsamples, mean = 0, sd = 1),num = 1:Nsamples)
samples$f = exp(-samples$x)
samples$f.approx = with(samples,exp(-1/2)*(1 - x + 1/2*(x^2-1) - 1/6*(x^3-3*x)))
meltDF <- melt(samples[,c('num','f','f.approx')],id.vars = "num")
ggplot(meltDF, aes(x=value,color=variable)) + geom_density(size=2) + scale_x_continuous(lim=c(0,5)) + scale_color_discrete('Variable')## Warning: Removed 524 rows containing non-finite values (stat_density).
## Warning: Removed 169 rows containing non-finite values (stat_density).
The mean of the approximation is 0.5977 compared to the value of 1.6245 from the actual function. The standard deviation for the approximation is 0.7705 compared to the value of 2.098 from the actual function.
The exponential is hard because it’s constants do not decay. Let’s use a simpler function: \[f(x) = (x-1)^2 + \frac{1}{10}\cos(x),\] with \(x\) distributed the same way as before.
Nsamples = 10000;
samples = data.frame(x = rnorm(Nsamples, mean = 0, sd = 1),num = 1:Nsamples)
samples$f = (samples$x-1)^2 + 0.1*cos(samples$x)
samples$f.approx = with(samples,(2+1/(10*exp(.5)) - 2*x + 1/2*(2-1/(10*exp(.5)))*(x^2-1)))
meltDF <- melt(samples[,c('num','f','f.approx')],id.vars = "num")
ggplot(meltDF, aes(x=value,color=variable)) + geom_density(size=2) + scale_x_continuous(lim=c(0,5)) + scale_color_discrete('Variable')## Warning: Removed 1093 rows containing non-finite values (stat_density).
## Warning: Removed 1095 rows containing non-finite values (stat_density).
The mean of the approximation is 2.0522 compared to the value of 2.0522 from the actual function. The standard deviation for the approximation is 2.3784 compared to the value of 2.3774 from the actual function.