Created an example specifically to test this feature

Example 'only-model' created based on the 'analytical-function' to test Inference without a surrogate.
Currently it's not running without building the surrogate - do tests first and see if able to refracture?

examples/only-model/analytical_function.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Nov 20 14:48:43 2019
+
+@author: farid
+"""
+import numpy as np
+import scipy.stats as stats
+import scipy.stats as st
+import seaborn as sns
+
+
+def analytical_function(xx, t=None):
+    """
+    Analytical Non-Gaussian Function
+
+    Authors: Farid Mohammadi, University of Stuttgart
+             Sergey Oladyshkin, University of Stuttgart
+    Questions/Comments: Please email Farid Mohammadi at:
+       farid.mohammadi@iws.uni-stuttgart.de
+
+    For function details and reference information, see:
+        https://doi.org/10.3390/e21111081
+
+    Parameters
+    ----------
+    xx : array
+        [x1, x2, ..., xn] where xn ~ Uinform(-5, 5).
+    t : array, optional
+        vector of times. The default is None. ( k − 1 ) /9 and k = 1,..., 10
+
+    Returns
+    -------
+    array
+        row vector of time vectors (s, t).
+
+    """
+    nParamSets, nParams = xx.shape
+
+    if t is None:
+        t = np.arange(0, 10, 1.) / 9
+
+    term1 = (xx[:, 0]**2 + xx[:, 1] - 1)**2
+
+    term2 = xx[:, 0]**2
+
+    term3 = 0.1 * xx[:, 0] * np.exp(xx[:, 1])
+
+    term5 = 0
+    if nParams > 2:
+        for i in range(2, nParams):
+            term5 = term5 + xx[:, i]**3/i
+
+    const = term1 + term2 + term3 + 1 + term5
+
+    # Compute time dependent term
+    term4 = np.zeros((nParamSets, len(t)))
+    for idx in range(nParamSets):
+        term4[idx] = -2 * xx[idx, 0] * np.sqrt(0.5*t)
+
+    Output = term4 + np.repeat(const[:, None], len(t), axis=1)
+
+    return {'x_values': t, 'Z': Output[0]}
+
+
+if __name__ == "__main__":
+
+    MCSize = 10000
+    ndim = 10
+    sigma = 2
+
+    # -------------------------------------------------------------------------
+    # ----------------------- Synthetic data generation -----------------------
+    # -------------------------------------------------------------------------
+    t = np.arange(0, 10, 1.) / 9
+
+    MAP = np.zeros((1, ndim))
+    synthethicData = analytical_function(MAP, t=t)
+
+    # -------------------------------------------------------------------------
+    # ---------------------- Generate Prior distribution ----------------------
+    # -------------------------------------------------------------------------
+
+    xx = np.zeros((MCSize, ndim))
+
+    params = (-5, 5)
+
+    for idxDim in range(ndim):
+        lower, upper = params
+        xx[:, idxDim] = stats.uniform(
+            loc=lower, scale=upper-lower).rvs(size=MCSize)
+
+    # -------------------------------------------------------------------------
+    # ------------- BME and Kullback-Leibler Divergence -----------------------
+    # -------------------------------------------------------------------------
+    Outputs = analytical_function(xx, t=t)
+
+    cov_matrix = np.diag(np.repeat(sigma**2, synthethicData.shape[1]))
+
+    Likelihoods = st.multivariate_normal.pdf(
+        Outputs['Z'], mean=synthethicData[1], cov=cov_matrix)
+
+    sns.kdeplot(np.log(Likelihoods[Likelihoods > 0]),
+                shade=True, color="g", label='Ref. Likelihood')
+
+    normLikelihood = Likelihoods / np.nanmax(Likelihoods)
+    # Random numbers between 0 and 1
+    unif = np.random.rand(1, MCSize)[0]
+
+    # Reject the poorly performed prior
+    accepted = normLikelihood >= unif
+
+    # Prior-based estimation of BME
+    logBME = np.log(np.nanmean(Likelihoods))
+    print(f'\nThe Naive MC-Estimation of BME is {logBME:.5f}.')
+
+    # Posterior-based expectation of likelihoods
+    postExpLikelihoods = np.mean(np.log(Likelihoods[accepted]))
+
+    # Calculate Kullback-Leibler Divergence
+    KLD = postExpLikelihoods - logBME
+    print("The Kullback-Leibler divergence estimation is {KLD:.5f}.")
+
+    # -------------------------------------------------------------------------
+    # ----------------- Save the arrays as .npy files -------------------------
+    # -------------------------------------------------------------------------
+    if MCSize > 500000:
+        np.save(f"data/refBME_KLD_{ndim}.npy", (logBME, KLD))
+        np.save(f"data/mean_{ndim}.npy", np.mean(Outputs['Z'], axis=0))
+        np.save(f"data/std_{ndim}.npy", np.std(Outputs['Z'], axis=0))
+    else:
+        np.save(f"data/Prior_{ndim}.npy", xx)
+        np.save(f"data/origModelOutput_{ndim}.npy", Outputs[1:])
+        np.save(f"data/validLikelihoods_{ndim}.npy", Likelihoods)

examples/only-model/bayesvalidrox/__init__.py

+# -*- coding: utf-8 -*-
+__version__ = "0.0.5"
+
+from .pylink.pylink import PyLinkForwardModel
+from .surrogate_models.surrogate_models import MetaModel
+from .surrogate_models.meta_model_engine import MetaModelEngine
+from .surrogate_models.inputs import Input
+from .post_processing.post_processing import PostProcessing
+from .bayes_inference.bayes_inference import BayesInference
+from .bayes_inference.bayes_model_comparison import BayesModelComparison
+from .bayes_inference.discrepancy import Discrepancy
+
+__all__ = [
+    "__version__",
+    "PyLinkForwardModel",
+    "Input",
+    "Discrepancy",
+    "MetaModel",
+    "MetaModelEngine",
+    "PostProcessing",
+    "BayesInference",
+    "BayesModelComparison"
+    ]

examples/only-model/bayesvalidrox/bayes_inference/__init__.py

+# -*- coding: utf-8 -*-
+
+from .bayes_inference import BayesInference
+from .mcmc import MCMC
+
+__all__ = [
+    "BayesInference",
+    "MCMC"
+    ]

examples/only-model/bayesvalidrox/bayes_inference/discrepancy.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import scipy.stats as stats
+from bayesvalidrox.surrogate_models.exp_designs import ExpDesigns
+
+
+class Discrepancy:
+    """
+    Discrepancy class for Bayesian inference method.
+    We define the reference or reality to be equal to what we can model and a
+    descripancy term \\( \\epsilon \\). We consider the followin format:
+
+    $$\\textbf{y}_{\\text{reality}} = \\mathcal{M}(\\theta) + \\epsilon,$$
+
+    where \\( \\epsilon \\in R^{N_{out}} \\) represents the the effects of
+    measurement error and model inaccuracy. For simplicity, it can be defined
+    as an additive Gaussian disrepancy with zeromean and given covariance
+    matrix \\( \\Sigma \\):
+
+    $$\\epsilon \\sim \\mathcal{N}(\\epsilon|0, \\Sigma). $$
+
+    In the context of model inversion or calibration, an observation point
+    \\( \\textbf{y}_i \\in \\mathcal{y} \\) is a realization of a Gaussian
+    distribution with mean value of \\(\\mathcal{M}(\\theta) \\) and covariance
+    matrix of \\( \\Sigma \\).
+
+    $$ p(\\textbf{y}|\\theta) = \\mathcal{N}(\\textbf{y}|\\mathcal{M}
+                                             (\\theta))$$
+
+    The following options are available:
+
+    * Option A: With known redidual covariance matrix \\(\\Sigma\\) for
+    independent measurements.
+
+    * Option B: With unknown redidual covariance matrix \\(\\Sigma\\),
+    paramethrized as \\(\\Sigma(\\theta_{\\epsilon})=\\sigma^2 \\textbf{I}_
+    {N_{out}}\\) with unknown residual variances \\(\\sigma^2\\).
+    This term will be jointly infer with the uncertain input parameters. For
+    the inversion, you need to define a prior marginal via `Input` class. Note
+    that \\(\\sigma^2\\) is only a single scalar multiplier for the diagonal
+    entries of the covariance matrix \\(\\Sigma\\).
+
+    Attributes
+    ----------
+    InputDisc : obj
+        Input object. When the \\(\\sigma^2\\) is expected to be inferred
+        jointly with the parameters (`Option B`).If multiple output groups are
+        defined by `Model.Output.names`, each model output needs to have.
+        a prior marginal using the `Input` class. The default is `''`.
+    disc_type : str
+        Type of the noise definition. `'Gaussian'` is only supported so far.
+    parameters : dict or pandas.DataFrame
+        Known residual variance \\(\\sigma^2\\), i.e. diagonal entry of the
+        covariance matrix of the multivariate normal likelihood in case of
+        `Option A`.
+
+    """
+
+    def __init__(self, InputDisc='', disc_type='Gaussian', parameters=None):
+        self.InputDisc = InputDisc
+        self.disc_type = disc_type
+        self.parameters = parameters
+
+    # -------------------------------------------------------------------------
+    def get_sample(self, n_samples):
+        """
+        Generate samples for the \\(\\sigma^2\\), i.e. the diagonal entries of
+        the variance-covariance matrix in the multivariate normal distribution.
+
+        Parameters
+        ----------
+        n_samples : int
+            Number of samples (parameter sets).
+
+        Returns
+        -------
+        sigma2_prior: array of shape (n_samples, n_params)
+            \\(\\sigma^2\\) samples.
+
+        """
+        self.n_samples = n_samples
+        ExpDesign = ExpDesigns(self.InputDisc)
+        self.sigma2_prior = ExpDesign.generate_ED(
+            n_samples, sampling_method='random', max_pce_deg=1
+            )
+        # Store BoundTuples
+        self.ExpDesign = ExpDesign
+
+        # Naive approach: Fit a gaussian kernel to the provided data
+        self.ExpDesign.JDist = stats.gaussian_kde(ExpDesign.raw_data)
+
+        # Save the names of sigmas
+        if len(self.InputDisc.Marginals) != 0:
+            self.name = []
+            for Marginalidx in range(len(self.InputDisc.Marginals)):
+                self.name.append(self.InputDisc.Marginals[Marginalidx].name)
+
+        return self.sigma2_prior

examples/only-model/bayesvalidrox/bayesvalidrox.mplstyle

+figure.titlesize : 30
+axes.titlesize : 30
+axes.labelsize : 30
+axes.linewidth : 3
+axes.grid : True
+lines.linewidth : 3
+lines.markersize : 10
+xtick.labelsize : 30
+ytick.labelsize : 30
+legend.fontsize : 30
+font.family : serif
+font.serif : Arial
+font.size : 30
+text.usetex : True
+grid.linestyle : -
+figure.figsize : 24, 16

examples/only-model/bayesvalidrox/desktop.ini

+[LocalizedFileNames]
+exploration.py=@exploration.py,0

examples/only-model/bayesvalidrox/post_processing/__init__.py

+# -*- coding: utf-8 -*-
+
+from .post_processing import PostProcessing
+
+__all__ = [
+    "PostProcessing"
+    ]

examples/only-model/bayesvalidrox/pylink/__init__.py

+# -*- coding: utf-8 -*-
+
+from .pylink import PyLinkForwardModel
+
+__all__ = [
+    "PyLinkForwardModel"
+    ]