Handling double dependencies in problem

[KillingVector]

Hi Jonathan! First off, thank you for creating this package, it has been very interesting to work with so far. My question is with framing a dependency problem. The idea is that in experiment, we might analyse two variables, but these two variables might be dependent on two other variables. Eg, analysing V and W such that: V(x,y) and W(x,y) x and y are independent in this case, and so, I can easily analyse their impact on some objective function. But I am interested in the behaviours of V and W.

I have tried to implement this as below. I think the issue I am having is that the implementation I have used works fine for the case when one variable is dependent on the other. That is: when the joint is assembled with x and V(x,y) - or some other like combination.

When I assemble the joint as below: cp.J( v, w ), I get the dependency error: "dangling dependencies". I have tried to include three, cp.J( x, v, w ), in order to capture at least one of the initial dependencies (though I am sure this is not a good way to frame it), and it fails with errors for under-defined dependencies.

Implementation

# test starter values
x = cp.Normal(0,1)
y = cp.Normal(1,1)
v = x+y
w = x*y # x-y
order = 4

joint = cp.J( v, w )
dist = [ cp.Normal(), cp.Normal ] # two variables
dist = cp.J( *dist )
P = cp.expansion.stieltjes( order, dist )

node_count = 2 * order + 2
node_sample = dist.sample( node_count, "sobol" )
nodes = joint.inv( dist.fwd( nodes_sample ) ) # fails here with dangling dependencies

evals = [ nodes[0] * nodes[1] for node in nodes.T ]
model = cp.git_regression( P, nodes, evals )

Additional I know that framing dependencies is not necessarily straightforward. I am simultaneously reading and learning in this area as I play around with your code - I am aware there could be understanding that I am lacking, and am still trying to build that understanding. I am hoping you might be able to provide direction in how I would do so with your code.

Please let me know if I can provide any additional information, and thank you for your time. (:

[jonathf]

Hi!

Sorry for the delayed response.

In chaospy handling of dependencies relies on the ability to create Rosenblatt transformation. In essense this means the dependencies has to be decomposable into P(X1), P(X2|X1), P(X3|X1,X2), .... In theory you can alsways create this, but in pratice it might be unfeasable to get there, and beyond the built-in multivariate distributions, chaospy can not automatically create this decomposition.

In your case, having P(V | X,Y) and P(W | X,Y) is not feasable to decompose on its own.

Your best bet, if you intend to go this way can be achieved in your case is to go through raw statistical moments. Doing so isn't really supported directly, but you can get there.

You need to create a function that outputs raw statistical moments for you variables. I.e. every i and j: E(V^i W^j) = E((X+Y)^i (X*Y)^j) = ... = mu_ij

Then make a custom (hack) distribution that produces these moments:

class MyDist(cp.Distribution):

    def __init__(self):
        dependencies, parameters, rotation = cp.declare_dependencies(self, {}, length=2)
        super().__init__(
            dependencies=dependencies,
            parameters=parameters,
            rotation=rotation,
        )

    def get_mom_parameters(self):
        return {}

    def _mom(self, loc):
        i, j = loc
        # calculate mu_ij
        return mu_ij

This can then be used to create orthogonal polynomials using Cholesky:

P = cp.expansion.cholesky(4, MyDist())

For the sampling, you can always do:

xy_nodes = x, y = dist.sample(...)
vw_nodes = np.array([x+y, x*y])

Hope this helps.

[jonathf]

You might want to check out Isserlis' theorem on how to calculate the moments you need.

[KillingVector]

Hi Jonathan, Thank you so much for taking the time to respond. I am working with skew and kurtosis elsewhere, so I at least have some familiarity with higher order moments. I will need to process what's you've responded with, and I do appreciate the detail!

[KillingVector]

Hi Jonathan! The solution I ended up coming up with was to follow the moment_approximation function. It does rely on the user having knowledge of how the V, W variables are related to the X,Y variables. I am still going through the motions of working out how to clean it up, and the approach for handling the relating equations in a clean fashion. There would be some work to do for it to be written in a more similar fashion to your own work, but would you be interested in this dependent distribution? I was planning to create a fork and develop it there.

[jonathf]

Good for you. :)

I think getting to specific with the dependent distribution doesn't make too much sense to merge. Fork sounds like a good way to go. I'll likely take a look if you go that route.