High order gPC estimators derived from non-standard distributions are unstable about chaospy HOT 8 CLOSED

I have just confirmed that the problem occurs also for Uniform distribution.

from chaospy.

It works for me if I use a rv with a mean of 0. I haven't used chaospy much, but generally high-ish order polynomials work better with mean-0, standard-deviation-1 abscissae. If you can afford to, work in such a space, and then translate any results back to the original space as needed by an appropriate (scaling and) shifting.

Unless this capability is supposed to be built-in to chaospy, which would make sense.

from chaospy.

As @tsbertalan said. Try to parameterize the problem in terms of standard distributions when going to higher order.

Polynomial chaos expansions allows you to address the problem trough a proxy variable:

from numpy import std
from chaospy import Normal, generate_quadrature, orth_ttr, fit_quadrature, Std

def f(x):
    return x

order = 6
rv = Normal(10, 0.1)
rv_proxy = Normal(0, 1)

data = f(rv.sample(order+1))

nodes, weights = generate_quadrature(order, rv_proxy, rule='G')
P, norms = orth_ttr(order, rv_proxy, normed=False, retall=True)
evals = f(rv.inv(rv_proxy.fwd(nodes[0])))
u_hat = fit_quadrature(P, nodes, weights, evals, norms=norms)
print(Std(u_hat, rv_proxy))
print(std(data))

This approximation works well as long as the mapping between rv and rv_proxy is smooth.

from chaospy.

I just want to say for now, thanks for your feedback and ideas. Both suggestions look like exactly what I was looking for. I will leave it open for now at least before I get the time to properly try it all out. Sorry for late responses.

from chaospy.

Right, I finally sat down and tested it out. Briefly, proxy method rocks. I focused on a variance and used the code below. I take that that the evals in your example simply transform to unit square through forward Rosenblatt and then transform to desired distribution through inverse Rosenblatt.

So clearly this is the property of the method. Is there a place in the literature where this has been described?

def f(x):
    return x**3


def direct(rv, order=4):
    nodes, weights = generate_quadrature(order, rv, rule='G')
    evals = f(nodes[0])

    Ps, norms = orth_ttr(order, rv, normed=False, retall=True)
    u_hat = fit_quadrature(Ps, nodes, weights, evals, norms=norms)
    return Var(u_hat, rv)


def proxy(rv, order=4):
    rv_proxy = Normal(0, 1)
    nodes, weights = generate_quadrature(order, rv_proxy, rule='G')
    evals = f(rv.inv(rv_proxy.fwd(nodes[0])))

    Ps, norms = orth_ttr(order, rv_proxy, normed=False, retall=True)
    u_hat = fit_quadrature(Ps, nodes, weights, evals, norms=norms)
    return Var(u_hat, rv_proxy)

rv = Normal(10, 0.1)

orders = arange(1, 11)
proxy_variance = array([proxy(rv, o) for o in orders])
direct_variance = array([direct(rv, o) for o in orders])

from chaospy.

That is correct, those are Rosenblatt transformations.

Rosenblatt is usually discussed in the context of transforming from and to stochastically independent variables, but the rules are exactly the same applied to achieve stability in higher orders. There should be something more specific to what you are asking out there, but unfortunately I don't remember what they are.

If you want to read up on the transformation itself, then the guys who made the DAKOTA software library discusses the topic here:
https://www.researchgate.net/profile/Michael_Eldred/publication/242068651_Recent_Advances_in_Non-Intrusive_Polynomial_Chaos_and_Stochastic_Collocation_Methods_for_Uncertainty_Analysis_and_Design/links/53d6abb30cf220632f3dc635.pdf

from chaospy.

I consider the issue as resolved.

If you consider that not to be the case, or if there are other problems related to the same issue, feel free to re-open the issue.

from chaospy.

Yes, Jonathan. Thanks a lot for your help here - the discussion was very useful to me. I should closed it before actually.

from chaospy.