Handling Correlated PDFs

PUQ assumes input parameters are independent, but sometimes this is not true. In that case, you need to have either a Joint PDF, or a large number of samples drawn from that Joint PDF. In this section we’ll show how to take two correlated output PDFs from one PUQ run and use them as inputs to another PUQ run.

These examples are from puq/examples/correlated. The first simulation, which creates obviously correlated outputs is c1_prog.py.

#!/usr/bin/env python

import optparse
from puqutil import dump_hdf5
import time

usage = "usage: %prog --x x --y y"
parser = optparse.OptionParser(usage)
parser.add_option("--x", type=float)
parser.add_option("--y", type=float)

(options, args) = parser.parse_args()
x = options.x
y = options.y

# need to produce two correlated output variables, f and g

f = x
g = x + y
z = f * g

dump_hdf5('f', f, "$f(x,y)=x+y$")
dump_hdf5('g', g, "$g(x,y)=x+y$")
dump_hdf5('z', z, "$z(x,y)=x(x+y)$")

This has three outputs, f, g, and z, which will be created as PDFs by PUQ.

The input script is

from puq import *

def run():
    x = UniformParameter('x', 'x', min=-5, max=5)
    y = UniformParameter('y', 'y', min=-5, max=5)

    host = InteractiveHost()
    uq = Smolyak([x,y], level=2)
    prog = TestProgram('./c1_prog.py')
    return Sweep(uq, host, prog)

Lets run PUQ and create the PDFs:

~/puq/examples/correlated> puq start -f c1.hdf5 c1.py
Saving run to c1.hdf5

Processing <HDF5 dataset "f": shape (13,), type "<f8">
        Surface   = 1.0*x
        RMSE      = 1.07e-12 (1.07e-11 %)

SENSITIVITY:
Var      u*            dev
-----------------------------
x    1.0000e+01    9.7391e-12
y    0.0000e+00    0.0000e+00

Processing <HDF5 dataset "g": shape (13,), type "<f8">
        Surface   = x + y
        RMSE      = 1.52e-12 (7.59e-12 %)

SENSITIVITY:
Var      u*            dev
-----------------------------
y    1.0000e+01    9.7391e-12
x    1.0000e+01    9.7391e-12

Processing <HDF5 dataset "z": shape (13,), type "<f8">
        Surface   = x*(x + y)
        RMSE      = 7.59e-12 (1.52e-11 %)

SENSITIVITY:
Var      u*            dev
-----------------------------
x    5.5178e+01    6.8073e+01
y    2.5000e+01    3.5355e+01

As you can see, PUQ correctly generated each response surface. Let’s look at z.

~/puq/examples/correlated> puq plot -v z c1.hdf5

OK. So far things look good. The second simulation takes f and g as input PDFs and outputs f*g. So it should look the same as z. However, f and g are very correlated so it won’t work as we expect. But before we try that, we need to save f and g to files.

~/puq/examples/correlated> puq analyze c1.hdf5

Save f and g to f.json and g.json.

You can look at c2_prog.py and verify it simply multiplies the two input parameters. Then take a look at c2.py.

from puq import *

def run():
    f = Parameter('f', 'f', pdf=LoadObj('f.json'), use_samples=True)
    g = Parameter('g', 'g', pdf=LoadObj('g.json'), use_samples=True)

    host = InteractiveHost()
    uq = Smolyak([f,g], level=2)
    prog = TestProgram('./c2_prog.py')
    return Sweep(uq, host, prog)

For the first run, change lines 4 and 5 so use_samples is False. This is the default and will cause the input parameters to be treated as independent. Then run puq

~/puq/examples/correlated> puq start -f c2_indep.hdf5 c2.py
Saving run to c2_indep.hdf5

Processing <HDF5 dataset "k": shape (13,), type "<f8">
        Surface   = 1.0*f*g
        RMSE      = 7.92e-14 (7.96e-14 %)

SENSITIVITY:
Var      u*            dev
-----------------------------
f    4.9775e+01    7.0000e+01
g    4.9498e+01    7.0000e+01

The response surface is correct (f*g). Looking at the output, k,

~/puq/examples/correlated> puq plot -v k c2_indep.hdf5

This looks nothing like the plot for z. Now let’s try it with correlated input parameters. Change c2.py back so lines 4 and 5 have the use_samples arg set to True.

from puq import *

def run():
    f = Parameter('f', 'f', pdf=LoadObj('f.json'), use_samples=True)
    g = Parameter('g', 'g', pdf=LoadObj('g.json'), use_samples=True)

    host = InteractiveHost()
    uq = Smolyak([f,g], level=2)
    prog = TestProgram('./c2_prog.py')
    return Sweep(uq, host, prog)

~/puq/examples/correlated> puq start -f c2.hdf5 c2.py
     Saving run to c2.hdf5

     Processing <HDF5 dataset "k": shape (13,), type "<f8">
             Surface   = 1.0*f*g
             RMSE      = 7.92e-14 (7.96e-14 %)

     SENSITIVITY:
     Var      u*            dev
     -----------------------------
     f    4.9775e+01    7.0000e+01
     g    4.9498e+01    7.0000e+01

     ~/puq/examples/correlated> puq plot -v k c2.hdf5

This looks very much like the plot for z we got in the first section. For a better comparison, you can use puq analyze to save z and k to files, then do

~/puq/examples/correlated> puq read -c z.json k.json