Hi everyone,
I’m currently building a likelihood estimator using BayesFlow v1. The model setup is relatively simple — we aim to estimate the likelihood p(y∣θ), where θ is a two-dimensional parameter vector and y is a scalar (one-dimensional) observation. The structure of the model is outlined below:
def simulator(theta):
x1 = np.random.default_rng().normal(loc=theta[0], scale=0.3)
x2 = np.random.default_rng().normal(loc=theta[1], scale=0.3)
g1x = (x1**2 * x2)/20-1
return g1x
def prior_m1(D=2, n_samples=2000):
return np.random.default_rng().uniform(low=0, high=10, size=(n_samples, D))
theta_samples = prior_m1(D=2, n_samples=20000)
data_samples = np.array([simulator(theta) for theta in theta_samples]).reshape(-1, 1)
I created 20,000 sets of training samples for this problem and trained the likelihood model using the following setup:
SETTINGS_LIK = {
"dense_args": dict(units=128, activation="relu", kernel_regularizer=tf.keras.regularizers.l2(1e-4)),
"num_dense": 2,
"spec_norm": False,
"mc_dropout": False,
"dropout": True,
"residual": False,
"dropout_prob": 0.2,
}
inference_net = InvertibleNetwork(num_params=1, num_coupling_layers=10,coupling_design='affine',
coupling_settings=SETTINGS_LIK,permutation="learnable")
amortizer = AmortizedLikelihood(inference_net, name="simple").
The model was successfully trained. The goal is to **compare the estimated distribution of y from the likelihood model to the true distribution obtained via Monte Carlo simulations, given the same θ. For example:
sims = {}
theta1,theta2 = 4, 7
sims["conditions"] = np.array([[theta1,theta2]]).astype(np.float32) # 5,1.5
sims["observables"] = np.zeros((1,180,1)).astype(np.float32)
sims["observables"][0,:,0] = np.arange(-14,4,0.1)
x_sim_s = amortizer.sample(sims, n_samples=2000)
from scipy.stats import gaussian_kde
from sklearn.mixture import GaussianMixture
N = int(1e4)
x1 = np.random.normal(theta1, 0.3, size=(N,))
x2 = np.random.normal(theta2, 0.3, size=(N,))
gx = (x1**2 * x2)/20-1
kde_true = gaussian_kde(gx)
x_vals = np.linspace(gx.min(), gx.max(), 500)
pdf_vals = kde_true(x_vals)
x_sim_s_flat = x_sim_s.ravel()
kde_predict = gaussian_kde(x_sim_s_flat)
x_kde = np.linspace(np.min(x_sim_s_flat), np.max(x_sim_s_flat), 500)
kde_vals = kde_predict(x_kde)
plt.figure(figsize=(8, 5))
plt.plot(x_vals, pdf_vals, color='k', label='True', linewidth=2)
plt.plot(x_kde, kde_vals, color='red', linestyle='--', label='KDE (Predicted)', linewidth=2)
plt.xlabel('gx', fontsize=14)
plt.ylabel('Density', fontsize=14)
plt.title('PDF Estimate of Composite Variable gx', fontsize=16)
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()
The model performance over test data is acceptable, as you see:
However, when I tried new theta values, the estimated pdf of y deviates far from the true pdf:
I’m currently investigating potential reasons behind this issue. Since the observation y is one-dimensional, I suspect that applying an affine coupling flow in the neural likelihood model may be problematic—specifically, it might not be meaningful to split a 1D input into two parts, in standard affine coupling layers. This could be one of the underlying causes.
I also experimented with neural spline flows and interleaved flow architectures, but none of them were able to train the model effectively.
Does anyone have suggestions or insights on how to handle likelihood modeling for scalar outputs using BayesFlow or similar flow-based approaches?
