Skip to content
Snippets Groups Projects
Commit 005ae4a7 authored by Farid Mohammadi's avatar Farid Mohammadi
Browse files

[tests][AnalyticalFunc] add example notebook.

parent d2ebdfd7
No related branches found
No related tags found
No related merge requests found
$X_1$,$X_2$
9.993116948359226370e-01,3.945526636751626182e-01
9.961447327121788486e-01,3.667186570442192428e-01
9.947646363203489495e-01,3.647464239325610147e-01
9.957567531711847275e-01,3.731759386892447572e-01
9.938740724688709394e-01,3.496238556702422340e-01
9.934295067779727040e-01,3.489816158549155500e-01
1.000000000000000000e+00,2.733114696357644147e-01
9.929163528111150105e-01,3.544359760567539208e-01
9.927050241098982486e-01,3.584158152985525980e-01
9.915077464334921542e-01,3.702678640580663316e-01
9.366935781467614630e-01,1.296011206606640376e-01
9.366778238387639677e-01,1.296333541870885453e-01
9.366737110183974702e-01,1.296417560822205883e-01
9.366715749203039509e-01,1.296461022957972054e-01
9.366705838329740796e-01,1.296480821342860246e-01
9.366707796164559507e-01,1.296475551494160605e-01
9.366843916343090548e-01,1.296184596592687976e-01
9.366660378655075059e-01,1.296580126549548295e-01
9.366679119152595545e-01,1.296539068017917851e-01
9.366689440746930329e-01,1.296517202892714382e-01
......@@ -11,7 +11,7 @@ import scipy.stats as st
import seaborn as sns
def AnalyticalFunction(xx, t=None):
def analytical_function(xx, t=None):
"""
Analytical Non-Gaussian Function
......@@ -68,56 +68,56 @@ if __name__ == "__main__":
MCSize = 10000 #1000000
ndim = 10
sigma = 2
# -------------------------------------------------------------------------
# ----------------------- Synthetic data generation -----------------------
# -------------------------------------------------------------------------
t = np.arange(0, 10, 1.) / 9
MAP = np.zeros((1, ndim))
synthethicData = AnalyticalFunction(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 = AnalyticalFunction(xx, t=t)
cov_matrix = np.diag(np.repeat(sigma**2, synthethicData.shape[1]))
Likelihoods = st.multivariate_normal.pdf(Outputs[1:], 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('\nThe Naive MC-Estimation of BME is %.5f.'%(logBME))
# 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 %.5f."%KLD)
# -------------------------------------------------------------------------
# ----------------- Save the arrays as .npy files -------------------------
# -------------------------------------------------------------------------
......
This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment