Hi @Eilidh_Radcliff,
You are currently calculating the PCE coefficients using the expectation y_\alpha=\mathbb{E}[\Psi_\alpha(\boldsymbol{X})\mathcal{M}(\boldsymbol{X})] (refer to Eq. (1.25) in the UQLab user manual), which you are computing via Monte Carlo integration. While this method is valid, it is not the most efficient approach. For example, UQLab approximates this expectation using quadrature which significantly reduces the number of samples needed. You can see the impact this has in the example below.
However, it is unclear from your example what exactly you are doing with UQLab. By default, UQLab calculates the coefficients using regression or LARS to be more specific. To make a meaningful comparison between your implementation and UQLab’s, it is also essential to ensure both use the same data points.
Unless for educational purposes, I would generally not recommend implementing PCE from scratch. It may be better to rely on existing implementations that incorporate state-of-the-art techniques. They are often challenging to develop and optimize but also deliver in many cases significantly better results.
I hope this brings some clarity!
Best regards,
Styfen
Click me
import numpy as np
from numpy.polynomial.hermite import hermgauss
from scipy.stats import norm, qmc
def Y(x):
return 5*x**3 - 3*x**2 + 2*x + 4
def H0(x):
return 1
def H1(x):
return x
def H2(x):
return (x**2 - 1) / np.sqrt(2)
def H3(x):
return (x**3 - 3*x) / np.sqrt(6)
# -------------------------------------------------------------------
# Use Gaussian quadrature
# -------------------------------------------------------------------
N = 4 # Number of quadrature points
nodes, weights = hermgauss(N)
# Transform to standard normal
x_samples_normal = np.sqrt(2) * nodes
Y_samples_normal = Y(x_samples_normal)
H0_samples = H0(x_samples_normal)
H1_samples = H1(x_samples_normal)
H2_samples = H2(x_samples_normal)
H3_samples = H3(x_samples_normal)
# Compute coefficients using quadrature
a0 = np.sum(weights * Y_samples_normal * H0_samples) / np.sqrt(np.pi)
a1 = np.sum(weights * Y_samples_normal * H1_samples) / np.sqrt(np.pi)
a2 = np.sum(weights * Y_samples_normal * H2_samples) / np.sqrt(np.pi)
a3 = np.sum(weights * Y_samples_normal * H3_samples) / np.sqrt(np.pi)
print("a0 = {:.2f}, a1 = {:.2f}, a2 = {:.2f}, a3 = {:.2f}".format(a0, a1, a2, a3))
# -------------------------------------------------------------------
# Use MC integration
# -------------------------------------------------------------------
N = 1_000_000 # Number of samples
sampler = qmc.LatinHypercube(d=1)
sample = sampler.random(n=N)
x_samples_normal = norm.ppf(sample) # Transform to standard normal
Y_samples_normal = Y(x_samples_normal)
H0_samples = H0(x_samples_normal)
H1_samples = H1(x_samples_normal)
H2_samples = H2(x_samples_normal)
H3_samples = H3(x_samples_normal)
a0 = np.mean(Y_samples_normal * H0_samples)
a1 = np.mean(Y_samples_normal * H1_samples)
a2 = np.mean(Y_samples_normal * H2_samples)
a3 = np.mean(Y_samples_normal * H3_samples)
print("a0 = {:.2f}, a1 = {:.2f}, a2 = {:.2f}, a3 = {:.2f}".format(a0, a1, a2, a3))