PCE Quadrature

pce-quadrature.jl

using Distributions
using FastGaussQuadrature

function P(x::Float64, d::Int)
    if d == 0
        return 1.0
    end

    if d == 1
        return x
    end

    return ((2d - 1) * x * P(x, d - 1) - (d - 1) * P(x, d - 2)) / float(d)
end

indices = [
    [0, 0],
    [1, 0],
    [0, 1],
    [2, 0],
    [1, 1],
    [0, 2],
    [3, 0],
    [2, 1],
    [1, 2],
    [0, 3],
    [4, 0],
    [3, 1],
    [2, 2],
    [1, 3],
    [0, 4],
]
function f(x::AbstractVector)
    return π * (x[1] - 1) * sin(π * x[1]) * (1 - x[2]^2)
end

function quadrature(p::Int)
    y = zeros(15)

    x, w = gausslegendre(p + 1)

    for i in 1:(p + 1)
        for j in 1:(p + 1)
            y +=
                (w[i] / 2) *
                (w[j] / 2) *
                f([x[i], x[j]]) *
                [P(x[i], ind[1]) * P(x[j], ind[2]) for ind in indices]
        end
    end
    return y
end

y = quadrature(6)

μ = y[1]
σ² = sum(y[2:end] .^ 2)

println("Mean: $μ")
println("Variance: $σ²")