PCE with least squares minimization

pce-least-squares.jl

using Distributions

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

p = 4

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],
]

x = rand(Uniform(-1, 1), 100, 2)

function f(x::AbstractVector)
    return prod(x)
end

M = map(f, eachrow(x))

A = mapreduce(row -> [prod(P.(row, ind)) for ind in indices], hcat, eachrow(x))'

y = inv(transpose(A) * A) * transpose(A) * M

μ = y[1]
global σ² = 0.0

for i in 2:length(y)
    global σ² +=
        (y[i] * sqrt(1 / (2 * indices[i][1] + 1)) * sqrt(1 / (2 * indices[i][2] + 1)))^2
end

M2 = mapreduce(row -> sum(y .* [prod(P.(row, ind)) for ind in indices]), hcat, eachrow(x))

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

ϵ = mean((M .- M2) .^ 2)
println("MSE: $ϵ")