---
title: "Polynomial Chaos Expansions"
author: "Ryan G. McClarren"
date: "November 4, 2014"
output:
  html_document:
    toc: yes
---

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.

```{r, include=FALSE}
library(ggplot2)
library(reshape)
```
```{r}
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')
```

The mean of the approximation is `r mean(samples$f.approx)` compared to the value of `r mean(samples$f)` from the actual function. The standard deviation for the approximation is `r sd(samples$f.approx)` compared to the value of `r sd(samples$f)` 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.

```{r, include=FALSE}
library(ggplot2)
library(reshape)
```
```{r}
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')
```

The mean of the approximation is `r mean(samples$f.approx)` compared to the value of `r mean(samples$f)` from the actual function. The standard deviation for the approximation is `r sd(samples$f.approx)` compared to the value of `r sd(samples$f)` from the actual function.