MCMC proposal distribution error

Python usage
1

Hello, I am trying to implement an MCMC Random-walk Metropolis-Hastings calibration procedure that includes gaussian model discrepancy. I’ve already worked with the RandomWalkMetropolisHastings class before without any trouble but I am now encountering an error I don’t understand. I am using OT1.23.

I have three calibration parameters (\theta, \sigma^{2}_{\delta}, \xi_{\delta}), the last two correspond to hyperparameters in the discrepancy covariance.

The supports of my distributions are the following :

support = ot.Interval([0.0, 1.0, 0.0], [15e-4, 100, 1e-5])

I use the following uniform proposal :

proposal = ot.ComposedDistribution([ot.Uniform(-step, step) for step in step_sizes])

with step-sizes :

step_sizes = [1e-3, 10, 1e-4]

When I try calling the ot.RandomWalkMetropolisHastings class with this proposal I get the following error :

TypeError: InvalidArgumentException : The range of the proposal distribution does not contain the origin with enough probability

I don’t understand what this error represents, I have tried tweaking with the step-sizes without any success. Since the origin is contained in my proposal’s supports, I really don’t see where the error comes from. My questions are the following :

  • Why do we have an error here?
  • What are the theoretical reasons for having apparently “stagnant” random walks (at least that’s how I understand that it contains the origin) with strong probability for the MCMC to work?
  • (Partially related) If someone has experience with bayesian calibration using model discrepancy: is there any rule of thumb for the proposal distribution in the hyperparameters?

Thank you very much for your help!
Edgar

2

Hi Edgar, could you send a standalone script allowing to reproduce this bug? I have tried but failed to reproduce your error message…

Concerning the two other points :

  • In a random-walk Metropolis-Hastings algorithm, convergence to the equilibrium distribution is automatically ensured for proposal distributions that are symmetrical with respect to the origin, which for continuous variables usually implies that their origin is in the proposal support. It is still possible to work with more exotic proposals, using for instance the UserDefinedMetropolisHastings class, but with less guarantees on the convergence.
  • Bayesian calibration using model discrepancy (Kennedy & O’Hagan style) is technically difficult and prone to identifiability issues. My personal choice is to use uniform priors with carefully chosen bounds; we will try to add a working example for this.

Best,
Merlin

3

Hello Merlin, thank you for your answer and your explanations.
Below is a script with some details. We first have a Gaussian Process (GP) model and some data

#GP metamodel
gp = kriging_result.getMetaModel()

#Size of the data, list of the form [data_times, data_values]
e = len(data)

We define the different covariance matrices,

  • the noise matrix is \Sigma_{\epsilon} = \hat{\sigma}^{2}_{\epsilon}I_{e}
  • the covariance of the GP model is \Sigma_{\text{GP}} = \{\gamma((t^{*}_{i},\theta_{i}),(t^{*}_{j},\theta_{j}))\}_{1\leq i,j \leq e}, where t_{1}^{*},\ldots,t_{e}^{*} are the experiment times and (t^{*}_{i},\theta_{i})_{1\leq i \leq e} are in the GP design of experiment and \gamma is the posterior kernel.
  • the parametric discrepancy kernel \Sigma_{\delta} = \{\sigma_{\delta}\exp(-\xi_{\delta}|t^{*}_{i} - t^{*}_{j}|)\}_{1\leq i,j \leq e} with new parameters \sigma_{\delta}, \xi_{\delta}
#Covariance of the data noise
sigma_data = np.diag([1]*e)*np.var(data_values)

#Covariance of the GP model
sigma_gp = np.empty([e,e])
for i in range(e):
    for j in range(e):
        sigma_gp[i,j] = np.asarray(np.matrix(kriging_result.getCovarianceModel()(doe[i], doe[j])))[0][0]

#Covariance of the discrepancy
def sigma_delta(sigma_delta, xi_delta):
    sigma_delta = np.empty([e,e])
    for i in range(e):
        for j in range(e):
            sigma_delta[i,j] = sigma_delta*np.exp(-xi_delta*np.abs(data_time[i]-data_time[j]))
    return sigma_delta

We proceed in defining the GP function and the posterior log-pdf, where we take non-informative Jeffreys priors on both p(\sigma_{\delta})\propto 1/\sigma_{\delta},\; p(\xi_{\delta}) \propto 1/\xi_{\delta} (maybe this can be changed?):

def gp_func(theta):
    global data_time, gp
    return(np.asarray(gp(np.asarray([[data_time[i]] + [theta[0]] for i in range(len(data_time))])))[:,0])

def log_density(theta):
    global sigma_data, sigma_gp

    diff = data_values - gp_func([theta[0]])

    #assemble the covariance-matrix
    sigma_delta = sigma_delta(theta[1], theta[2])
    sigma_total = sigma_data + sigma_gp + sigma_delta 
    sigma_total_inv = np.linalg.inv(sigma_total)

    #compute the weighted loss
    wlss = np.dot(diff, np.dot(sigma_total_inv, diff))

    log_pdf = np.log(1/theta[1]) + np.log(1/theta[2]) - ss/2 
    return [log_pdf]

#define the Python function for MCMC
log_density = ot.PythonFunction(3, 1, log_density)

We now choose the parameters for the MCMC. I chose the values of the support for the discrepancy parameters more or less arbitrarily, based on the scale of my data.

#Random-walk step-sizes
step_sizes = [1e-3, 10, 1e-4]

#Define the support of the posterior density
support = ot.Interval([0.0, 1.0, 0.0], [15e-4, 10, 1e-5])

#Uniform random-walk proposal
proposal = ot.ComposedDistribution([ot.Uniform(-step, step) for step in step_sizes])

#Random initialization for (\theta, \sigma_\delta, \xi_\delta)
rand_theta0 = ot.ComposedDistribution([ot.Uniform(0, 15e-4), ot.Uniform(1, 100), ot.Uniform(0.0, 1e-5)]).getSample(1)

Finally we can run the OT RWMH class:

RWMHsampler = ot.RandomWalkMetropolisHastings(log_density, support, rand_theta0, proposal)

and I get the aforementioned error …

TypeError: InvalidArgumentException : The range of the proposal distribution does not contain the origin with enough probability