md"""

```math

\begin{align*}

F &= \int_a^b f(x) dx  \\

​

\eta &= \frac{f(\xi)}{p(\xi)}, \int_a^b p(\xi)d\xi = 1 \\

​

\text{여기서 $p$ 는 probability density function(PDF)} \\

E\left[\eta\right] &= E\left[ \frac{f(\xi)}{p(\xi)} \right] = \int_a^b \frac{f(x)}{p(x)}p(x) dx = \int_a^b f(x)dx = I \\

​

E\left[\eta\right] &= F = E\left[ \frac{f(\xi)}{p(\xi)} \right] \approx \frac{1}{N}\sum_{i=0}^{N-1}\frac{f(x)}{p(x)} \\

F &= \int_a^b f(x) dx \approx \langle F^N \rangle=\frac{1}{N}\sum_{i=0}^{N-1}\frac{f(X_i)}{p(X_i)}, X_i \sim p \\

\text{$X_i$는 p 분포에서 샘플링한 값}

​

\end{align*}

```

"""

using Distributions, Plots, Random, LaTeXStrings

md"""

###### $p(x) =  \frac{2}{\pi}$ 인 경우 즉 Uniform distribution

```math

\begin{align*}

f(x) &= \sin(x) \\

p(x) &= \color{red} \frac{2}{π} \\

CDF(x) & = \frac{2}{π}x \\

r_i &= 0 \sim 1 = CDF(x_i) = \frac{2}{π}x_i\\

x_i &= 0 \sim \frac{\pi}{2} =  \color{red}\frac{π}{2}r_i \\

\langle F^N \rangle &= \frac{1}{N}\sum_{i=0}^{N-1}\frac{\sin(x_i)}{p(x_i)} \\

\end{align*} \\

​

```

"""

md"""

###### $p(x) =  \frac{8x}{\pi^2}$ 인 경우 즉 $f(x)$ 와 유사한 분포를 선택

```math

\begin{align*}

f(x) &= \sin(x) \\

p(x) &= \color{red} \frac{8}{\pi^2}x \\

CDF(x) & = \frac{4x^2}{\pi^2} \\

r_i &= 0 \sim 1 = CDF(x_i) = \frac{4x_i^2}{\pi^2}\\

x_i &= 0 \sim \frac{\pi}{2} = \color{red}\frac{\pi}{2}\sqrt{r_i} \\

\langle F^N \rangle &= \frac{1}{N}\sum_{i=0}^{N-1}\frac{\sin(x_i)}{p(x_i)} \\

\end{align*} \\

​

```

"""

Random.seed!(13)

f(x) = sin(x)

p₁(x) = 2/π

p₂(x) = 8x/π^2

N=16

epoch = 100

begin

    const π² = π*π

    const ∑ = sum

    const σ = std

    const μ = mean  

end

begin

    x = 0:0.01:π/2

    plot(x,x->2/π,linewidth=2,label=L"p_1(x)=2/\pi")

    plot!(x,x->8x/π²,linewidth=2,label=L"p_2(x)=8x/\pi^2")

    plot!(x,x->sin(x),linewidth=2,label=L"f(x)=sin(x)",legend=:bottomright)

end

begin   

    arrF₁ = Array{Float64,1}()

    arrF₂ = Array{Float64,1}()  

    for i in 1:epoch

        r = rand(N)

        x₁ = π * r / 2

        x₂ = π * .√(r) / 2

        F₁ = ∑(sin.(x₁)./p₁.(x₁)) / N

        F₂ = ∑(sin.(x₂)./p₂.(x₂)) / N

        push!(arrF₁,F₁)

        push!(arrF₂,F₂)

    end 

    plot([1:epoch],x->1,line = (:dot, 2),label=L"F")

    plot!(arrF₁,label=L"uniform\;sampling:\langle F^N_1 \rangle",linewidth=2,marker =(:hex,3))

    plot!(arrF₂,label=L"importance\;sampling:\langle F^N_2 \rangle",linewidth=2,yticks=(-0.5:0.1:1.5),marker=(:circle,3),legend=:topright ) 

end

μ(arrF₁), μ(arrF₂)

σ(arrF₁),σ(arrF₂)

begin

    mse₁ = ∑( (arrF₁ .- 1).^2 ) / N

    mse₂ = ∑( (arrF₂ .- 1).^2 ) / N

    mse₁,mse₂

end

​

"Mean squre error of (arrF₁, arrF₂) : ($mse₁, $mse₂)"