Variance Estimation

[Martin-Lacroix]

Dear Dr.Feinberg,

I am working on polynomial chaos expansion for a very short time and I recently discovered your Python library Chaospy. I realized some expansions were not delivering the same result as for a Monte-Carlo simulation concerning the variance. A minimal working example is an equation of the type y(t) = exp(-k² t) where k is a random uniform variable in [-1,1] and t is the time.

orderBase = 5
orderQuad = 15
dist = chaospy.Uniform(-1,1)
point = numpy.random.uniform(-1,1,nbrPts)

def y(x,k):

    y = numpy.exp(-k**2*x)
    return y

point,weight = chaospy.generate_quadrature(order=orderQuad,dist=dist,rule="Gaussian")
sample = [y(x,*point) for point in point.T]

poly,norms = chaospy.orth_ttr(order=orderBase,dist=dist,sort='G',retall=True)
model = chaospy.fit_quadrature(poly,point,weight,sample)

mean = chaospy.E(model,dist)
var = chaospy.Var(model,dist)

I do not really understand what is it happening, some other probability distributions like Uniform[0,1] seem to work correctly.

[jonathf]

That is an excellent question. I usually don't have no answer for questions like these.

My best educated guess is that the problem you are describing happens to be a bad fit for polynomial approximation, so the polynomial and quadrature order has to be increased both to the point where the method becomes unstable. But as it work better with [0, 1], something looks like it going on.

I will not give up just yet. Will see if I can solve this.

[Martin-Lacroix]

Thank you, I figured out this problem occurs especially with more complicated function and using a gamma law (Increasing the number of polynomials in this simple example can lead to correct results before observing instabilities).

Unrelated question: May I ask you where did you find the analytical expressions of the three term recurrence coefficients used for the construction of polynomial bases in your code ? I have some difficulties in finding papers providing clear information about them.

[jonathf]

Yes, I figured out that it was problem specific, and that increasing order helps. But I am still curious as to why the convergence rate drops so much as it does for this problem.

A little bit here and there. Most of the standard TTRs are gathered from the paper "The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations" by D. Xiu and G. Karniadakis.