Polynomial representations and random variables

[robertsawko]

Hi,

Many thanks for contributing this package. I've done a few basic examples and I like your API a lot. At the moment I am trying to construct a random variable based on the polynomial representation. Is there a way to do it in ChaosPy already?

For instance, if I have a code which looks like this

a = Uniform(a_low, a_high)
I = Uniform(I_low, I_high)
dist = J(a, I)

P, norms = orth_ttr(order, dist, retall=True)
nodes, weights = generate_quadrature(
       order, dist, sparse=False, rule='G')
solves = [u(t, s[0], s[1]) for s in nodes.T]
u_hat = fit_quadrature(P, nodes, weights, solves, norms=norms)

I would like to make a random variable based on u_hat. If I just use algebraic operations on a and I get an expected behaviour of a transformed RV, but how to perform this with the u_hat representation?

Obviously, I can print u_hat and see it in terms of q0 and q1, but I would like to make an actual substitution and see in terms of a and I and then for instance plot a pdf.

Please let me know if it's there already. Otherwise, I am happy to work on a pull request to implement it.

[jonathf]

I have thought about creating a variable as you suggest, but struggled with generalising the idea.
Chaospy builds on top of Rosenblatt transformations from the [0,1]^D hypercube. And unless the resulting polynomial is seperable, making a dist becomes difficult.

For your use case, your output is a univariate variable, the job is more doable. But I need a couple of components that likely have to be estimated numerically:

• the resulting cumulative distribution function
• practical lower and upper bounds of the RV.

The easy solution, though perhaps not very elegant is to use kernel density estimation on samples generated from your polynomial. This is already implemented:

N = 1000
samples =  u_hat(*dist.sample(N))
dist_poly = cp.SampleDist(samples)

If you have a better approach, let me know, and we can discuss making a PR.

[robertsawko]

I see that constructing the PDF in the general case is going to be quite difficult, yes. I suppose for now, I can do with the KDE approach though. It's still less difficult than estimating the same from the model. Thanks for the comment.

Closing for now.

[flo2k]

Hi,

here is this function called cp.QoI_Dist() that can do the construction of the variable (internally, it uses also a KDE) for you.

qoi_dist = cp.QoI_Dist(u_hat, dist)

Now you can use the pdf() function of qoi_dist, depending if qoi_dist is a single variable or an array. If it is an array you have to use indices like qoi_dist[i].pdf(...).

Best,

Florian

[robertsawko]

Is cp.QoI_Dist() working? This is my example:

from chaospy import (
    Uniform, QoI_Dist, generate_quadrature, fit_quadrature, orth_ttr, J)
from numpy import e


def u(t, a, I):
    return I * e**(-a * t)

t = 10.0
a_low, a_high = 0.0,  0.1
I_low, I_high = 8.0, 10.0
a = Uniform(a_low, a_high)
I = Uniform(I_low, I_high)
dist = J(a, I)
order = 4

P, norms = orth_ttr(order, dist, retall=True)
nodes, weights = generate_quadrature(order + 1, dist, rule='G')
solves = [u(t, s[0], s[1]) for s in nodes.T]
u_hat = fit_quadrature(P, nodes, weights, solves, norms=norms)

qoi_dist = QoI_Dist(u_hat, dist)

qoi_dist does not appear to have a pdf function. Instead it is ndarray of objects. Also mean and std functions don't work when applied on this object.

[jonathf]

Hm... It seems stuck inside numpy array.

You can extract it by doing

qoi_dist = qoi_dist.item()

I'll also make a bugfix in the code for that.

[jonathf]

Should be solved in #24.

Feel free to re-open the issue if that is not the case.