jonathf / chaospy

Chaospy - Toolbox for performing uncertainty quantification.
https://chaospy.readthedocs.io/
MIT License
439 stars 87 forks source link

Dependency using copula #154

Open ahenkes opened 5 years ago

ahenkes commented 5 years ago

I have four normal distributions and want to add dependency to them in a simple way and than do some PCE. Lets say, I have the the correlation matrix from experimental data, my guess would be to use a nataf copula with the respective R matrix. Unfortunately, for more than two random variables the memory gets full. Maybe this is not the canonical way to do it? To stress my point, here is some of my code

def polynomial_chaos_expansion(random_variables):
    import chaospy as chp
    import numpy as np

    standardDev_1 = random_variables[0] * 0.1
    standardDev_2 = random_variables[1] * 0.1
    standardDev_3 = random_variables[2] * 0.1
    standardDev_4 = random_variables[3] * 0.1

    distNormal_1 = chp.Normal(random_variables[0], standardDev_1)
    distNormal_2 = chp.Normal(random_variables[1], standardDev_2)
    distNormal_3 = chp.Normal(random_variables[2], standardDev_3)
    distNormal_4 = chp.Normal(random_variables[3], standardDev_4)

    distribution = chp.J(distNormal_1, distNormal_2, distNormal_3, distNormal_4)

    R = np.zeros(4)
    R[0, 1] = 0.75
    R[1, 0] = R[0, 1]
    nataf = chp.Nataf(distribution, R, ordering=None)

    quad_dim = 1
    pol_dim = 1
    abscissas, weights = chp.generate_quadrature(quad_dim, nataf, rule="F")
    polynomial_expansion = chp.orth_chol(pol_dim, nataf)

    solves = np.zeros((abscissas[0].size, 9))
    counter_pce = 0
    for ab in abscissas.T:
        solves[counter_pce, :] = fun(ab)
        counter_pce = counter_pce + 1

    count_approx = 0
    matpar_expectation = np.zeros((9, 1))
    matpar_deviation = np.zeros((9, 1))

    for loop_gp in range(9):
        approx = \
               chp.fit_quadrature(
                        polynomial_expansion, abscissas, weights, solves[:, loop_gp])
        expected_pce = chp.E(approx, nataf)
        deviation_pce = chp.Std(approx, nataf)
        matpar_expectation[count_approx, 0] = expected_pce
        matpar_deviation[count_approx, 0] = deviation_pce
        count_approx = count_approx + 1

    return matpar_expectation, matpar_deviation
jonathf commented 5 years ago

That is a bug. I will take a closer look to see what is going on.

In addition:

ahenkes commented 5 years ago

Thank you very much for the quick answer!

Indeed, it is a typo, it should be eye(4)! I will check your attached slides and report back!

ahenkes commented 5 years ago

So, I used the proposed procedure from your slides, and it seems to work out, which is great! I now try to rebuild my previous results using this technique.

Thank you very much for your great work!

jonathf commented 5 years ago

Glad to hear.

I'll keep the issue open for the sake of the original bug.