using Plots

using LaTeXStrings

using QuadGK

f(x) = x.*exp.(x)

integral,err = quadgk(x->f(x),0,1,rtol=1e-8)

using Distributions, Random

Random.seed!(12345)

udist = Uniform(0,1)

begin   

    N1 = 100

    x1 = rand(udist,N1) 

    Ihat1 = sum(f(x1))/N1

end

begin

    N2 = 5000

    x2 = rand(udist,N2)

    Ihat2 = sum(f(x2))/N2

end

md"""

Approximating the expected value of the Beta distribution

​

$\mathbb{E}[x] = \int_{p(x)}p(x)xdx$ \

where $x \sim p(x) = Beta(\alpha_1, \alpha_2)$

​

**Step 1**: we identify $h(x) = x$ \

**Step 2**: the function $g(x)$ is simply the probability density function $p(x)$ due the expression for $p(x)$ above:

​

$p(x)^* = \frac{p(x)}{\int p(x)dx} = p(x)$

​

**Step 3** we can use Julia to easily draw $N$ independent sample $p(x)$ using the function $Beta$. And finally \

**Step 4** we approximate the expectation with the expression

​

$\mathbb{E}_{Beta(\alpha_1,\alpha_2)} \approx \frac{1}{N}\sum_i x_i$

​

​

"""

begin

    alpha1 = 2

    alpha2 = 10

    N = 100000

    betaDist = Beta(alpha1,alpha2)

    x_b = rand(betaDist,N)

end

# Montecarlo expectation

expectMC = mean(x_b)

# Analytic Expression for Beta Mean

expectAnalytic = alpha1 / (alpha1 + alpha2)

md"""

#### Monte Carlo Approximation for Optimization

​

Monte Carlo Approximation can also be used to solve optimization problems of the form:

​

$\hat{x} = \underset{x\in(a,b)}{\operatorname{argmax}} \;g(x)$

​

if g(x) fulfills the same criteria described above (namely that it is a scaled version of a probability distribution), then (as above) we can define the probability function

​

$p(x) = \frac{g(x)}{C}$

​

This allows us to instead solve problem

​

$\hat{x} = \underset{x}{\operatorname{argmax}}\;Cp(x)$

​

if we can sample form $p(x)$, the solution $\hat{x}$ is easily found by drawing samples from $p(x)$ and determining the location of the samples that has the highest density (Note that the solution is not dependent of the value of $C$). The following example demonstrates Monte Carlo optimization:

"""

md"""

#### Example: Monte Carlo Optimization of $g(x) = e^{-\frac{(x-4)^2}{2}}$

​

Say we would like to find the value of $x_{opt}$ which optimizes the function $g(x) = e^{-\frac{(x-4)^2}{2}}$. In other words we want to solve the problem

​

$x_{opt} = \underset{x}{\operatorname{argmax}}\; e^{-\frac{(x-4)^2}{2}}$

​

We could solve for $x_{opt}$ using standard calculus methods, but a clever trick is to use Monte Carlo approximation to solve the problem. First, notice that $g(x)$ is a scaled version of a Normal distribution with mean equal to 4 and unit variance:

​

$g(x) = C \times \frac{1}{\sqrt{2\pi}}e^{-\frac{(x-4)^2}{2}}$

$= C \times \mathcal{N}(4,1)$

​

where $C=\sqrt{2\pi}$. This means we can solve for $x_{opt}$ by drawing samples from the normal distribution and determining where those samples have the highest density. The following chunck of Julia code solves the optimization problem in this way

"""

# Monte Carlo Optimization of exp(x-4)^2

begin

    N₃ = 100000

    x₃ = 0:0.1:6

    C₃ = sqrt(2*pi)

    g₃(x) = exp.(-0.5(x.-4)^2 ) 

end

using StatsBase

begin

    y₃ = pdf.(Normal(4,1),x₃)

    # Calculate Monte Carlo Approximation

    x_normal = rand(Normal(4,1),N₃)

    bins = range(min(x_normal...),stop=max(x_normal...),length=100)

    h = fit(Histogram,x_normal,bins)

    xHat = bins[argmax(h.weights)]

    bar(bins[1:99],h.weights ./ max(h.weights...),fillcolor=:yellow,fillalpha=0.2,linewidth=0.1,label="Normalized Histogram")

    plot!(x₃,g₃.(x₃)/C₃,linecolor=:red,linewidth=3,label=L"g_3(x)")

    plot!(x₃,g₃.(x₃),linecolor=:blue,linewidth=2, label=L"p(x)=\frac{g_3(x)}{C}")   

    plot!([4,4],[0,1],linecolor=:black,linewidth=2,label=L"x_{opt}")

    plot!([xHat,xHat],[0,1],linecolor=:lightgreen,linewidth=2,label=L"\hat{x}") 

end

begin

    hh = fit(Histogram,x_normal,bins);  

    bar(hh)

end