diff --git a/example/poisson_dirichlet_example.py b/example/poisson_dirichlet_example.py
index d55ddea7..8e70bae7 100644
--- a/example/poisson_dirichlet_example.py
+++ b/example/poisson_dirichlet_example.py
@@ -105,7 +105,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
     misfit_form = PoissonMisfitForm(d, noise_variance)
     misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form, [bc0])
     prior_mean = dlx.fem.Function(Vh_m)
-    prior_mean.x.array[:] = 0.01
+    prior_mean.x.array[:] = np.log(2)
     prior_mean = prior_mean.x
 
     prior = hpx.BiLaplacianPrior(
@@ -242,43 +242,21 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
         prior_samples.append(prior_sample)
         posterior_samples.append(posterior_sample)
 
-    with dlx.io.XDMFFile(
+    with dlx.io.VTXWriter(
         msh.comm,
-        "poisson_Dirichlet_prior_Bilaplacian_samples_prior_np{0:d}.xdmf".format(nproc),
-        "w",
-    ) as file:
-        file.write_mesh(msh)
-        for i, sample in enumerate(prior_samples):
-            time_value = i / (num_samples_generate - 1)
-            file.write_function(sample, time_value)
-
-    with dlx.io.XDMFFile(
+        "poisson_Dirichlet_prior_Bilaplacian_samples_prior_np{0:d}.bp".format(nproc),
+        prior_samples,
+    ) as vtx:
+        vtx.write(0.0)
+
+    with dlx.io.VTXWriter(
         msh.comm,
-        "poisson_Dirichlet_prior_Bilaplacian_samples_posterior_np{0:d}.xdmf".format(
+        "poisson_Dirichlet_prior_Bilaplacian_samples_posterior_np{0:d}.bp".format(
             nproc
         ),
-        "w",
-    ) as file:
-        file.write_mesh(msh)
-        for i, sample in enumerate(posterior_samples):
-            time_value = i / (num_samples_generate - 1)
-            file.write_function(sample, time_value)
-
-    # with dlx.io.VTXWriter(
-    #     msh.comm,
-    #     "poisson_Dirichlet_prior_Bilaplacian_samples_prior_np{0:d}.bp".format(nproc),
-    #     prior_samples,
-    # ) as vtx:
-    #     vtx.write(0.0)
-
-    # with dlx.io.VTXWriter(
-    #     msh.comm,
-    #     "poisson_Dirichlet_prior_Bilaplacian_samples_posterior_np{0:d}.bp".format(
-    #         nproc
-    #     ),
-    #     posterior_samples,
-    # ) as vtx:
-    #     vtx.write(0.0)
+        posterior_samples,
+    ) as vtx:
+        vtx.write(0.0)
 
     eigen_decomposition_results = {"A": Hmisfit, "B": prior, "k": k, "d": d, "U": U}
 
diff --git a/example/poisson_dirichlet_example_reg.py b/example/poisson_dirichlet_example_reg.py
index 2b9e927d..50a8649b 100644
--- a/example/poisson_dirichlet_example_reg.py
+++ b/example/poisson_dirichlet_example_reg.py
@@ -109,7 +109,7 @@ def top_bottom_boundary(x: Sequence[float]) -> Sequence[bool]:
     misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form, [bc0])
 
     prior_mean = dlx.fem.Function(Vh_m)
-    prior_mean.x.array[:] = 0.01
+    prior_mean.x.array[:] = np.log(2)
 
     prior_gamma = prior_param["gamma"]
     prior_delta = prior_param["delta"]
diff --git a/example/poisson_example.py b/example/poisson_example.py
index 0ef0e7b2..905887eb 100644
--- a/example/poisson_example.py
+++ b/example/poisson_example.py
@@ -87,7 +87,7 @@ def run_inversion(
     misfit_form = PoissonMisfitForm(d, noise_variance)
     misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form)
     prior_mean = dlx.fem.Function(Vh_m)
-    prior_mean.x.array[:] = 0.01
+    prior_mean.x.array[:] = np.log(2)
     prior_mean = prior_mean.x
 
     prior = hpx.BiLaplacianPrior(
@@ -254,9 +254,10 @@ def run_inversion(
 if __name__ == "__main__":
     nx = 64
     ny = 64
-    noise_variance = 1e-4
+    noise_variance = 1e-6
     prior_param = {"gamma": 0.07, "delta": 0.7}
     final_results = run_inversion(nx, ny, noise_variance, prior_param)
+
     k, d = (
         final_results["eigen_decomposition_results"]["k"],
         final_results["eigen_decomposition_results"]["d"],
diff --git a/example/poisson_example_reg.py b/example/poisson_example_reg.py
index 6b8034aa..a658b948 100644
--- a/example/poisson_example_reg.py
+++ b/example/poisson_example_reg.py
@@ -85,7 +85,7 @@ def run_inversion(
     misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form)
 
     prior_mean = dlx.fem.Function(Vh_m)
-    prior_mean.x.array[:] = 0.01
+    prior_mean.x.array[:] = np.log(2)
 
     prior_gamma = prior_param["gamma"]
     prior_delta = prior_param["delta"]
@@ -216,7 +216,7 @@ def run_inversion(
 if __name__ == "__main__":
     nx = 64
     ny = 64
-    noise_variance = 1e-4
+    noise_variance = 1e-6
     prior_param = {"gamma": 0.02, "delta": 0.2}
     final_results = run_inversion(nx, ny, noise_variance, prior_param)
     k, d = (