BYO Prior Parameter Ensemble

[gmerritt123]

Hi, I'm trying to understand how to point my pst.pestpp_options['ies_parameter_ensemble'] to a covariance matrix that comes out of some sort of my own contrived geostatistical analysis. More or less trying to understand the workflow and maybe to a lesser extent exactly what the numbers in that covariance matrix mean.

My basic example below is just gonna attempt to use a distance matrix as a surrogate for parameter covariance (very very basic just to flesh out the perceived workflow).

import pandas as pd
import numpy as np
from scipy.spatial import distance_matrix

#say I have a dataframe of my 10 pilot point locations
hru_df = pd.DataFrame(data={'X':np.random.random(10)*100,'Y':np.random.random(10)*100})
hru_df['parID'] = 'hru_'+hru_df.index.astype(str)

#populate an inverse distance matrix as surrogate for covariance?
dm = 1/distance_matrix(hru_df[['X','Y']].to_numpy(),hru_df[['X','Y']].to_numpy())
#basically set the "range" parameter --> anything further than 25 m should not be correlated
#make up a sill value to "strengthen" the correlation
sill = 25
dm = np.where(dm<1/25,dm,0) * sill
print(pd.DataFrame(data=dm)) #pretty print as df to visualize
#put into Cov instance???
jcb = pyemu.mat.mat_handler.Cov(x=dm,names=hru_df['parID'].to_list(),isdiagonal=True) #it IS diagonal eh?
jcb.to_binary("dm.jco")
# then just point my Pst instance (pst) to this???
pst.pestpp_options['ies_parameter_ensemble'] = Pst.ParameterEnsemble.from_binary("dm.jco")
#???
#profit

Am I on the right track with this workflow --> i.e. will it incorporate my imposed correlation between parameters into the ensemble generation/sampling?

[jtwhite79]

I think you want to use parcov instead of ies_parameter_ensemble if you want o BYO prior parameter covariance matrix - IES can read it and draw from it. If you want to draw from your own cov matrix to get an ensemble, I use an eigen factorization of the matrix to form a projection matrix (FE^1/2), then you can simply matrix-vector multiply a vector of standard-normal variates thru this projection matrix to get a realization of correlated parameter values that have the correct variance. This is what ies does under the hood. Note however, this approach really only works for <30K pars - after that you run out of memory (ies will do things by parameter group, assuming each group is an independent block of the prior parameter covariance matrix, but all that handling can be slow).
Also, if you want, pyemu can generate covariance matrix for you using variogram and x-y pairs - its all pure python so could follow it with a debugger to see what it is doing...

[gmerritt123]

Ah yes sorry and thanks for interpreting/clarifying that --> I am waffling between creating parcov and letting ies do the sampling under the hood OR doing my initial sampling beforehand WITH a cov matrix, and then pointing ies to that result (i.e. ies_parameter_ensemble). I'm leaning towards the latter at least for this first run through simply to be able to visualize and understand how different variogram choices and parameter grouping can affect my initial ensembles.

So I think I have this working like this:

v = pyemu.utils.geostats.ExpVario(contribution=20,a=1000) #sill and range args

dm = v.covariance_matrix(x=hru_df['X'],y=hru_df['Y'],names=hru_df['parID'])

Also, if you want, pyemu can generate covariance matrix for you using variogram and x-y pairs - its all pure python so could follow it with a debugger to see what it is doing...

That'd be taking a look at the covariance_matrix method on v eh?

Then to sample a set of ensembles that follow this covariance:

pe = pyemu.ParameterEnsemble.from_gaussian_draw(pst=pst,cov=dm,num_reals=10)

And looking at the produced ensemble values they appear to be spatially correlated, and the range sampled fits with the "95% confidence" par_sigma_range default (I think?). Cool.

So I think that covers me off nicely for building a parcov (and subsequent ensemble) "manually" with one parameter group. The next step would be mixing different param groups with potentially differing variograms... like you mentioned IES will hold them independent of each other, which is okay... and that kinda means I supply a covariance matrix for each parameter group and then loop through them sequentially calling ParameterEnsemble.from_gaussian_draw on pe each time? That should have the effect of independently sampling each parameter group and updating the ensemble dataframe? Is there a better way/am I thinking about this workflow correctly?

Finally, the most complicated question (and kind of an open question I'd appreciate your take/advice on): how to go about incorporating covariance between parameter groups that have different units? Like the key case @userbarryb and I are thinking about is that recharge and surficial K should be correlated. We'd need to somehow build up a covariance matrix that captures the spatial correlation but also a correlation between each location's K parameter and Recharge parameter. Any thoughts/advice on the maths/implementation of something like that would be most helpful. Thanks!

[jtwhite79]

I prefer to generate my own prior ensembles before ies runs just to have some more control of it. The way I usually make a combined prior ensemble (with both geostat covariance and indep parameters) is with pyemu.helpers.geostatistical_draws() - it assumes the sill of the geostatistical structure(s) is 1.0 and that the parameter bounds approximate the prior parameter uncertainty (by default it assumes bounds approximate 4 standard deviations). This lets you have parameters of varying prior uncertainty but also correlated. I wouldnt use grid-scale parameters with that function tho (unless there are less than 20K pars per group) otherwise, memory trouble...

Your last question re correlation between HK and recharge...do you think they are correlated in a prior sense or are you expecting them to be correlated in the responses of the model (like high HK and high recharge can yield similar sim water levels as low HK and low recharge)? If the latter, then I'd let that correlation arise naturally in the inverse problem and not hard-wire it. If you really insist on adding this correlation in the prior, it can get kind of dicey bc to you have to start adding off-diagonal elements - it can be done but Ive had a hard time with this previously bc the matrix quickly becomes singular...