diff --git a/.gitattributes b/.gitattributes
index 3e5cf230..23dcd25c 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -2,3 +2,4 @@
 examples/**/*.npz filter=lfs diff=lfs merge=lfs -text
 examples/**/*.npy filter=lfs diff=lfs merge=lfs -text
 examples/**/*.dat filter=lfs diff=lfs merge=lfs -text
+
diff --git a/.gitignore b/.gitignore
index 35e613a9..d19426db 100644
--- a/.gitignore
+++ b/.gitignore
@@ -171,6 +171,16 @@ runs_sc25/*
 settings.json
 *.out
 
+# data files that can easily be generated
+examples/brainiac/**/*.npy
+examples/brainiac/**/*.npz
+
+examples/g_ar1/**/*.npy
+examples/g_ar1/**/*.npz
+
+examples/p_ar1/**/*.npy
+examples/p_ar1/**/*.npz
+
 # Profiler outputs
 *.nsys-rep
 *.qdstrm
diff --git a/examples/g_ar1/generate_data.py b/examples/g_ar1/generate_data.py
new file mode 100644
index 00000000..e9bb0a6b
--- /dev/null
+++ b/examples/g_ar1/generate_data.py
@@ -0,0 +1,101 @@
+import os
+import sys
+
+import numpy as np
+import scipy.sparse as sp
+from scipy.sparse.linalg import spsolve, spsolve_triangular
+from scipy.sparse import csc_matrix
+from scipy.linalg import cholesky
+
+BASE_DIR = os.path.dirname(os.path.abspath(__file__))
+
+if __name__ == "__main__":
+
+    np.random.seed(5)
+    n = 1000
+
+    ## define priors
+    s2 = 5  
+    tau = 1 / s2
+    phi = 0.9  
+    # noise obs
+    obs_noise_prec = 100
+    theta_original = [
+        phi,
+        tau,
+        obs_noise_prec,
+    ]
+
+    denom = s2 * (1 - phi**2)
+
+    diag = [(1 + phi**2) / denom] * n
+    diag[0] = diag[-1] = 1 / denom
+    off_diag = [-phi / denom] * (n - 1)
+
+    Q = sp.diags([diag, off_diag, off_diag], [0, -1, 1])
+    
+    # Compute sparse Cholesky factorization: Q = L @ L.T
+    # For tridiagonal matrix, we can use dense Cholesky on small blocks or scipy
+    Q_csc = Q.tocsc()
+    
+    print("Q shape:", Q.shape, "Q nnz:", Q.nnz)
+    print("Q sparsity:", 100 * Q.nnz / (Q.shape[0] * Q.shape[1]), "%")
+    print(Q.toarray()[:6, :6])
+    
+    # Method 1: Use dense Cholesky (for moderate sizes this is still efficient)
+    Q_dense = Q.toarray()
+    L_dense = cholesky(Q_dense, lower=True)
+    L = csc_matrix(L_dense)
+    
+    print("L nnz:", L.nnz, "L sparsity:", 100 * L.nnz / (L.shape[0] * L.shape[1]), "%")
+    
+    # Efficient sampling: generate z ~ N(0,I), then solve L @ u = z
+    z = np.random.normal(0, 1, size=n)
+    
+    # Solve L @ u = z using sparse triangular solver
+    u = spsolve_triangular(L, z, lower=True)
+    
+    # Verify the sampling worked correctly
+    print("Sample u statistics - mean:", np.mean(u), "std:", np.std(u), ". Should be around sqrt(s2) =", np.sqrt(s2))
+    
+    intercept = 2
+
+    x = np.concatenate((u, [intercept]))
+    print("x: ", x[:10])
+    
+    np.save("reference_outputs/x_original.npy", x)
+    np.save("inputs_ar1/x.npy", u)
+    np.save("reference_outputs/theta_original.npy", theta_original)
+
+    a_ar1 = sp.eye(n)
+    sp.save_npz("inputs_ar1/a.npz", a_ar1)
+
+    a_regression = sp.csr_matrix(np.ones((n, 1)))
+    sp.save_npz("inputs_regression/a.npz", a_regression)
+
+    eta = a_ar1 @ u + intercept
+
+    print("eta: ", eta[:6])
+    np.save("inputs_ar1/x_original.npy", eta)
+
+    noise = np.random.normal(0, np.sqrt(1 / obs_noise_prec), size=eta.shape)
+    print("noise: ", noise[:10])
+    y = eta + noise
+    np.save("y.npy", y)
+
+    print("y: ", y[:10])
+
+    Qprior = sp.block_diag([Q, sp.csr_matrix([[0.001]])])
+
+    a = sp.hstack([a_ar1, a_regression]) # a_ar1 #
+    Qcond = Qprior + obs_noise_prec * a.T @ a
+    print("Qcond: \n", Qcond.toarray()[:6,:6])
+
+    b = obs_noise_prec * a.T @ y
+    print("b: ", b[:10])
+    #x_est = np.linalg.solve(Qcond.toarray(), b)
+    x_est = spsolve(csc_matrix(Qcond), b)
+    print("norm(x - x_est): ", np.linalg.norm(x - x_est))
+
+    print("norm(eta - eta_est): ", np.linalg.norm(a @ x - a @ x_est))
+    print("normalized norm(eta - eta_est): ", np.linalg.norm(a @ x - a @ x_est) / np.linalg.norm(a @ x))
\ No newline at end of file
diff --git a/examples/g_ar1/run.py b/examples/g_ar1/run.py
new file mode 100644
index 00000000..2bddfe30
--- /dev/null
+++ b/examples/g_ar1/run.py
@@ -0,0 +1,153 @@
+import sys
+import os
+
+parent_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
+sys.path.append(parent_dir)
+
+import numpy as np
+
+from dalia import xp, sp, backend_flags
+from dalia.configs import likelihood_config, dalia_config, submodels_config
+from dalia.core.model import Model
+from dalia.core.dalia import DALIA
+from dalia.submodels import AR1SubModel, RegressionSubModel
+from dalia.utils import get_host, print_msg, scaled_logit  # , extract_diagonal
+from examples_utils.parser_utils import parse_args
+
+BASE_DIR = os.path.dirname(os.path.abspath(__file__))
+
+if __name__ == "__main__":
+
+    np.random.seed(3)
+
+    # load reference output
+    theta_original = np.load("reference_outputs/theta_original.npy")
+    print(
+        "theta original: ",
+        theta_original,
+    )
+
+    theta_initial = theta_original #[0.6, 1.0, 3.0]
+    print("theta initial: ", theta_initial)
+
+    x_original = np.load("reference_outputs/x_original.npy")
+    print("x original: ", x_original[:10])
+    print("dim(x original): ", x_original.shape)
+
+    ar1_dict = {
+        "type": "ar1",
+        "input_dir": f"{BASE_DIR}/inputs_ar1",
+        "phi": theta_initial[0],  # has to be between 0 and 1
+        "tau": theta_initial[1],  # assume to already be in log-scale
+        "ph_phi": {"type": "beta", "alpha": 5.0, "beta": 1.0},
+        "ph_tau": {"type": "gaussian", "mean": 0.0, "precision": 0.5},
+    }
+    ar1 = AR1SubModel(
+        config=submodels_config.parse_config(ar1_dict),
+    )
+
+    # Configurations of the regression submodel
+    regression_dict = {
+        "type": "regression",
+        "input_dir": f"{BASE_DIR}/inputs_regression",
+        "n_fixed_effects": 1,
+        "fixed_effects_prior_precision": 0.001,
+    }
+    regression = RegressionSubModel(
+        config=submodels_config.parse_config(regression_dict),
+    )
+
+    likelihood_dict = {
+        "type": "gaussian",
+        "prec_o": xp.log(theta_initial[2]),
+        # "prior_hyperparameters": {
+        #     "type": "penalized_complexity",
+        #     "alpha": 0.01,
+        #     "u": 5,
+        # },
+        "prior_hyperparameters": {"type": "gaussian", "mean": xp.log(theta_original[2]), "precision": 0.05},
+    }
+
+    model = Model(
+        submodels=[ar1, regression], #
+        likelihood_config=likelihood_config.parse_config(likelihood_dict),
+    )
+    print_msg(model)
+
+    Qprior = model.construct_Q_prior()
+    print("Qprior: \n", Qprior.toarray()[:6, :6])
+    Qinv = np.linalg.inv(Qprior.toarray())
+    geom_mean = np.exp(np.mean(np.log(Qinv.diagonal())))
+    print("Geometric mean of Qinv diagonal: ", geom_mean)
+
+    # in gaussian case x = 0, thus eta = 0
+    x_i = np.zeros(model.n_latent_parameters)
+    eta = model.a @ x_i
+    Qcond = model.construct_Q_conditional(eta=eta)
+    print("Qcond: \n", Qcond.toarray()[:6, :6])
+
+    b = model.construct_information_vector(eta=eta, x_i=x_i)
+    print("b: ", b[:10])
+
+    x_est = np.linalg.solve(Qcond.toarray(), b)
+    #print("x est: ", x_est)
+    print("norm(x_original - x_est): ", np.linalg.norm(x_original - x_est))
+
+    # Configurations of DALIA
+    dalia_dict = {
+        "solver": {"type": "dense"},
+        "minimize": {
+            "max_iter": 100,
+            "gtol": 1e-3,
+            "disp": True,
+            "maxcor": len(model.theta_external),
+        },
+        "f_reduction_tol": 1e-3,
+        "theta_reduction_tol": 1e-4,
+        "inner_iteration_max_iter": 50,
+        "eps_inner_iteration": 1e-3,
+        "eps_gradient_f": 1e-3,
+        "simulation_dir": ".",
+        "verbosity": 0,
+    }
+
+    dalia = DALIA(
+        model=model,
+        config=dalia_config.parse_config(dalia_dict),
+    )
+
+    print("initial model theta: ", model.theta_external)
+
+    print("\nCalling DALIA.run()")
+    results = dalia.run()
+
+    print_msg("\n--- Results ---")
+
+    theta = results["theta"]
+    print("theta:          ", np.round(theta, 4))
+    print("theta original: ", theta_original)
+
+    print_msg("Covariance of theta:\n", results["cov_theta"])
+    print_msg(
+        "Mean of the fixed effects:\n",
+        results["x"][-model.submodels[-1].n_fixed_effects :],
+    )
+
+    # print("x:          ", results["x"])
+    # print("x_original: ", x_original)
+
+    # print("eta: ", model.a @ x_original)
+    # print("eta est: ", model.a @ results["x"])
+
+    print_msg("\n--- Comparisons ---")
+    print("norm(eta - eta_est): ", np.linalg.norm(model.a @ x_original - model.a @ results["x"]))
+    print("normalized norm(eta - eta_est): ", np.linalg.norm(model.a @ x_original - model.a @ results["x"]) / np.linalg.norm(model.a @ x_original))
+
+    # Compare marginal variances of latent parameters
+    var_latent_params = results["marginal_variances_latent"]
+    Qconditional = dalia.model.construct_Q_conditional(eta=model.a @ model.x)
+    Qinv_ref = xp.linalg.inv(Qconditional.toarray())
+    print_msg(
+        "Norm (marg var latent - ref):    ",
+        f"{np.linalg.norm(var_latent_params - xp.diag(Qinv_ref)):.4e}",
+    )
\ No newline at end of file
diff --git a/examples/gr/preprocessing.py b/examples/gr/generate_data.py
similarity index 88%
rename from examples/gr/preprocessing.py
rename to examples/gr/generate_data.py
index 51330182..b5561e26 100644
--- a/examples/gr/preprocessing.py
+++ b/examples/gr/generate_data.py
@@ -28,12 +28,12 @@
     x = L_Sigma_prior @ z
     a = sparse.random(n_observations, n_latent_parameters, density=0.5)
 
-    theta_observations = np.log(3)
+    theta_observations = 3.0
     print(f"theta_observations: {theta_observations}")
     theta_likelihood: dict = {"theta_observations": theta_observations}
 
     # generate x from a gaussian distribution of dimensions n_latent_parameters with mean 0 and precision exp(theta_observations)
-    variance = 1 / np.exp(theta_observations)
+    variance = 1 / theta_observations
     eta = a @ x
     y = np.random.normal(eta, scale=np.sqrt(variance), size=n_observations)
     print(f"x: {x}")
@@ -53,7 +53,7 @@
     sparse.save_npz(f"{path}/inputs/a.npz", a)
 
     # save original latent parameters
-    np.save(f"{path}/inputs/x_original.npy", x)
+    np.save(f"{path}/reference_outputs/x_ref.npy", x)
 
     # save original hyperparameter theta
-    np.save(f"{path}/inputs/theta_original.npy", theta_observations)
+    np.save(f"{path}/reference_outputs/theta_ref.npy", theta_observations)
diff --git a/examples/gr/inputs/y.npy b/examples/gr/inputs/y.npy
new file mode 100644
index 00000000..a5755615
--- /dev/null
+++ b/examples/gr/inputs/y.npy
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4f5b972c648de0244b4fe3af2aac7e87210706d8ac879007011e061c2330c0ba
+size 8128
diff --git a/examples/gr/reference_outputs/theta_ref.npy b/examples/gr/reference_outputs/theta_ref.npy
index ccb2de10..01fbab03 100644
--- a/examples/gr/reference_outputs/theta_ref.npy
+++ b/examples/gr/reference_outputs/theta_ref.npy
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:722e1a2b1465e290a283836da486e044a051dc5c7d0b420a1fc0addfd4411c1f
+oid sha256:7b9d65f6717a0911c1cdf1cc680cb167f2ef1395be2db8ca2065cb2bd8a0a8f4
 size 136
diff --git a/examples/gr/reference_outputs/x_ref.npy b/examples/gr/reference_outputs/x_ref.npy
index fbf38892..87f3421e 100644
--- a/examples/gr/reference_outputs/x_ref.npy
+++ b/examples/gr/reference_outputs/x_ref.npy
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:86bdf16672dbf80add11a285c0fe68f2532436b2d232a86992144017d02fcd61
+oid sha256:4ddd3a0b4a76c50a77c276c3c7a3a37b506a633b3e31564e8fd3bf859d3e2184
 size 176
diff --git a/examples/gr/run.py b/examples/gr/run.py
index d12fdaad..78e613e0 100644
--- a/examples/gr/run.py
+++ b/examples/gr/run.py
@@ -35,8 +35,9 @@
     # Likelihood
     likelihood_dict = {
         "type": "gaussian",
-        "prec_o": 1.5,
-        "prior_hyperparameters": {"type": "gaussian", "mean": 3.5, "precision": 0.5},
+        "prec_o": 1.0,
+        "prior_hyperparameters": {"type": "gamma", "alpha": 2.0, "beta": 2.0},
+        #"prior_hyperparameters": {"type": "gaussian", "mean": 1.0, "precision": 0.5},
     }
     # Creation of the first model by combining the Regression submodel and the likelihood
     model = Model(
@@ -69,7 +70,8 @@
     results = dalia.run()
 
     print_msg("\n--- Results ---")
-    print_msg("Theta values:\n", results["theta"])
+    print_msg("Theta values external:\n", results["theta"])
+    print_msg("Theta values internal:\n", results["theta_internal"])
     print_msg("Covariance of theta:\n", results["cov_theta"])
     print_msg(
         "Mean of the fixed effects:\n",
@@ -79,6 +81,7 @@
     print_msg("\n--- Comparisons ---")
     # Compare hyperparameters
     theta_ref = xp.load(f"{BASE_DIR}/reference_outputs/theta_ref.npy")
+    print_msg("Reference theta:", theta_ref)
     print_msg(
         "Norm (theta - theta_ref):        ",
         f"{np.linalg.norm(results['theta'] - get_host(theta_ref)):.4e}",
@@ -101,7 +104,7 @@
     )
 
     # Compare marginal variances of observations
-    var_obs = dalia.get_marginal_variances_observations(theta=theta_ref, x_star=x_ref)
+    var_obs = dalia.get_marginal_variances_observations(theta_external=theta_ref, x_star=x_ref)
     var_obs_ref = extract_diagonal(model.a @ Qinv_ref @ model.a.T)
     print_msg(
         "Norm (var_obs - var_obs_ref):    ",
diff --git a/examples/gs_small/run.py b/examples/gs_small/run.py
index e94f20d1..d4dadd1e 100644
--- a/examples/gs_small/run.py
+++ b/examples/gs_small/run.py
@@ -58,6 +58,8 @@
         likelihood_config=likelihood_config.parse_config(likelihood_dict),
     )
 
+    print_msg(model)
+
     # Configurations of DALIA
     dalia_dict = {
         "solver": {"type": "dense"},
@@ -65,7 +67,7 @@
             "max_iter": args.max_iter,
             "gtol": 1e-3,
             "disp": True,
-            "maxcor": len(model.theta),
+            "maxcor": len(model.theta_external),
         },
         "f_reduction_tol": 1e-3,
         "theta_reduction_tol": 1e-4,
diff --git a/examples/gst_coreg2_small/run.py b/examples/gst_coreg2_small/run.py
index 55328b44..1dba8a7f 100644
--- a/examples/gst_coreg2_small/run.py
+++ b/examples/gst_coreg2_small/run.py
@@ -206,7 +206,7 @@
             "max_iter": args.max_iter,
             "gtol": 1e-3,
             "disp": True,
-            "maxcor": len(coreg_model.theta),
+            "maxcor": len(coreg_model.theta_external),
         },
         "f_reduction_tol": 1e-3,
         "theta_reduction_tol": 1e-4,
diff --git a/examples/gst_coreg3_small/run.py b/examples/gst_coreg3_small/run.py
index 8406d013..fa0fd90d 100644
--- a/examples/gst_coreg3_small/run.py
+++ b/examples/gst_coreg3_small/run.py
@@ -266,7 +266,7 @@
             "max_iter": args.max_iter,
             "gtol": 1e-3,
             "disp": True,
-            "maxcor": len(coreg_model.theta),
+            "maxcor": len(coreg_model.theta_external),
         },
         "f_reduction_tol": 1e-3,
         "theta_reduction_tol": 1e-4,
diff --git a/examples/gst_small/run.py b/examples/gst_small/run.py
index 2c36f753..24899cec 100644
--- a/examples/gst_small/run.py
+++ b/examples/gst_small/run.py
@@ -52,12 +52,13 @@
     likelihood_dict = {
         "type": "gaussian",
         "prec_o": 4,
-        # "prior_hyperparameters": {"type": "gaussian", "mean": 1.4, "precision": 0.5},
-        "prior_hyperparameters": {
-            "type": "penalized_complexity",
-            "alpha": 0.01,
-            "u": 4,
-        },
+        "prior_hyperparameters": {"type": "gamma", "alpha": 2.0, "beta": 2.0},
+        #"prior_hyperparameters": {"type": "gaussian", "mean": 1.4, "precision": 0.5},
+        # "prior_hyperparameters": {
+        #     "type": "penalized_complexity",
+        #     "alpha": 0.01,
+        #     "u": 4,
+        # },
     }
 
     # Creation of the model by combining the submodels and the likelihood
@@ -74,7 +75,7 @@
             "max_iter": args.max_iter,
             "gtol": 1e-3,
             "disp": True,
-            "maxcor": len(model.theta),
+            "maxcor": len(model.theta_external),
         },
         "f_reduction_tol": 1e-3,
         "theta_reduction_tol": 1e-4,
@@ -88,11 +89,18 @@
         config=dalia_config.parse_config(dalia_dict),
     )
 
-    results = dalia.run()
+    results = dalia.minimize()
 
     print_msg("\n--- Results ---")
+    theta_ref = np.load(f"{BASE_DIR}/reference_outputs/theta_ref.npy")
+    theta_external = np.array([0.08423457, 2.52313066, 1.46267965, np.exp(1.36076756)])
+
     print_msg("Theta values:\n", results["theta"])
-    print_msg("Covariance of theta:\n", results["cov_theta"])
+    print_msg("Theta values internal:\n", results["theta_internal"])
+    
+    #cov_theta = dalia.compute_covariance_hp(results["theta"])
+    cov_theta = dalia.compute_covariance_hp(theta_external)
+    print_msg("Covariance of theta:\n",  cov_theta["internal"])
     print_msg(
         "Mean of the fixed effects:\n",
         results["x"][-model.submodels[-1].n_fixed_effects :],
@@ -100,7 +108,6 @@
 
     print_msg("\n--- Comparisons ---")
     # Compare hyperparameters
-    theta_ref = np.load(f"{BASE_DIR}/reference_outputs/theta_ref.npy")
     print_msg(
         "Norm (theta - theta_ref):        ",
         f"{np.linalg.norm(results['theta'] - get_host(theta_ref)):.4e}",
diff --git a/examples/p_ar1/generate_data.py b/examples/p_ar1/generate_data.py
new file mode 100644
index 00000000..c1eef39d
--- /dev/null
+++ b/examples/p_ar1/generate_data.py
@@ -0,0 +1,90 @@
+import os
+import sys
+
+import numpy as np
+import scipy.sparse as sp
+
+from scipy.stats import multivariate_normal, poisson
+
+BASE_DIR = os.path.dirname(os.path.abspath(__file__))
+
+if __name__ == "__main__":
+
+    n = 1000
+
+    ## define priors
+    s2 = 1  # 0.7
+    tau = 1 / s2
+    ### note: phi between -1 and 1 for discrete timesteps 
+    # (doesn't make sense for negative in cts case)
+
+    ## rescale beta prior 2 theta - 1
+    phi = 0.9  # 0.9
+    theta_original = [phi, tau]
+
+    denom = s2 * (1 - phi**2)
+
+    diag = [(1 + phi**2) / denom] * n
+    diag[0] = diag[-1] = 1 / denom
+    off_diag = [-phi / denom] * (n - 1)
+
+    Q = sp.diags([diag, off_diag, off_diag], [0, -1, 1])
+    L = np.linalg.cholesky(Q.toarray())
+    Cov = np.linalg.inv(Q.toarray())
+
+    geom_mean = np.exp(np.mean(np.log(Cov.diagonal())))
+    print("Geometric mean of Qinv diagonal: ", geom_mean)
+
+    print(Q.toarray()[:6, :6])
+    print(np.linalg.inv(Q.toarray())[:6, :6])
+    print(np.round(Q.toarray() @ np.linalg.inv(Q.toarray()), 6)[:6, :6])
+    # exit()
+
+    mv = multivariate_normal(mean=np.zeros(n), cov=Cov, seed=3)
+
+    intercept = 2
+    u = mv.rvs()
+    print("u: ", u[:10])
+    eta = u + intercept
+    x = np.concatenate((u, [intercept]))
+    np.save("reference_outputs/x_original.npy", x)
+    x_initial = u + np.random.normal(0, 0.3, size=len(u))
+    np.save("inputs_ar1/x.npy", u)
+    np.save("reference_outputs/theta_original.npy", theta_original)
+
+    a_ar1 = sp.eye(n)
+    sp.save_npz("inputs_ar1/a.npz", a_ar1)
+
+    a_regression = sp.csr_matrix(np.ones((n, 1)))
+    sp.save_npz("inputs_regression/a.npz", a_regression)
+
+    print("eta: ", eta[:10])
+    np.save("inputs_ar1/x_original.npy", eta)
+
+    # sample with repitition
+    E = np.random.choice([1, 2, 3], size=n, replace=True)
+    # E = [1] * n
+    np.save("e.npy", E)
+
+    y = poisson.rvs(E * np.exp(eta), random_state=3)
+    np.save("y.npy", y)
+
+    print("y[:10]: ", y[:10])
+
+    Qprior = sp.block_diag([Q, sp.csr_matrix([[0.001]])])
+    # print("Qprior : \n", Qprior.toarray())
+
+    # a = sp.hstack([a_ar1, a_regression])
+    # Qcond = Qprior + obs_noise_prec * a.T @ a
+    # print("Qcond: \n", Qcond.toarray())
+
+    # b = obs_noise_prec * a.T @ y
+    # # -xp.exp(theta) * (eta - y)
+    # # beta_initial + np.linalg.solve(
+    # #     Qconditional.toarray(), information_vector
+    # # )
+    # x_est = np.linalg.solve(Qcond.toarray(), b)
+    # print("x_est: ", x_est)
+
+    # print("eta est : ", a @ x_est)
+    # print("eta :     ", a @ x)
diff --git a/examples/p_ar1/run.py b/examples/p_ar1/run.py
new file mode 100644
index 00000000..22edf823
--- /dev/null
+++ b/examples/p_ar1/run.py
@@ -0,0 +1,128 @@
+import sys
+import os
+
+parent_dir = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
+sys.path.append(parent_dir)
+
+import numpy as np
+
+from dalia import xp, sp, backend_flags
+from dalia.configs import likelihood_config, dalia_config, submodels_config
+from dalia.core.model import Model
+from dalia.core.dalia import DALIA
+from dalia.submodels import AR1SubModel, RegressionSubModel
+from dalia.utils import get_host, print_msg, scaled_logit  # , extract_diagonal
+from examples_utils.parser_utils import parse_args
+
+
+BASE_DIR = os.path.dirname(os.path.abspath(__file__))
+
+if __name__ == "__main__":
+
+    n = 1000
+
+    # load reference output
+    theta_original = np.load("reference_outputs/theta_original.npy")
+    print("theta original: ", theta_original)
+
+    x_original = np.load("reference_outputs/x_original.npy")
+    print("x original: ", x_original[:10])
+    print("dim(x original): ", x_original.shape)
+
+    ar1_dict = {
+        "type": "ar1",
+        "input_dir": f"{BASE_DIR}/inputs_ar1",
+        "n_latent_parameters": n,
+        "phi": theta_original[0],  # has to be between 0 and 1
+        "tau": xp.log(theta_original[1]),  # assume to already be in log-scale
+        "ph_phi": {"type": "beta", "alpha": 5.0, "beta": 1.0},
+        "ph_tau": {"type": "gaussian", "mean": 0.0, "precision": 0.5},
+    }
+    ar1 = AR1SubModel(
+        config=submodels_config.parse_config(ar1_dict),
+    )
+
+    # Configurations of the regression submodel
+    regression_dict = {
+        "type": "regression",
+        "input_dir": f"{BASE_DIR}/inputs_regression",
+        "n_fixed_effects": 1,
+        "fixed_effects_prior_precision": 0.001,
+    }
+    regression = RegressionSubModel(
+        config=submodels_config.parse_config(regression_dict),
+    )
+
+    likelihood_dict = {
+        "type": "poisson",
+        "input_dir": f"{BASE_DIR}",
+    }
+
+    model = Model(
+        submodels=[ar1, regression],
+        likelihood_config=likelihood_config.parse_config(likelihood_dict),
+    )
+    print_msg(model)
+
+    Qprior = model.construct_Q_prior()
+    print("Qprior: \n", Qprior.toarray())
+    Qinv = np.linalg.inv(Qprior.toarray())
+    geom_mean = np.exp(np.mean(np.log(Qinv.diagonal())))
+    print("Geometric mean of Qinv diagonal: ", geom_mean)
+
+    eta = model.a @ model.x
+    Qcond = model.construct_Q_conditional(eta=eta)
+
+    # L = np.linalg.cholesky(Qcond.toarray())
+
+    # plt.spy(Qcond, markersize=2)
+    # plt.title("Sparsity pattern of Qcond")
+    # plt.show()
+
+    # Configurations of DALIA
+    dalia_dict = {
+        "solver": {"type": "dense"},
+        "minimize": {
+            "max_iter": 100,
+            "gtol": 1e-3,
+            "disp": True,
+            "maxcor": len(model.theta_external),
+        },
+        "f_reduction_tol": 1e-3,
+        "theta_reduction_tol": 1e-4,
+        "inner_iteration_max_iter": 50,
+        "eps_inner_iteration": 1e-3,
+        "eps_gradient_f": 1e-3,
+        "simulation_dir": ".",
+    }
+
+    dalia = DALIA(
+        model=model,
+        config=dalia_config.parse_config(dalia_dict),
+    )
+
+    print("theta external: ", model.theta_external)
+    # print("x : ", model.x)
+    f_value = dalia._evaluate_f(model.theta_external)
+    print("after evaluate f. x: ", model.x)
+
+    results = dalia.minimize()
+
+    theta_unscaled = results["theta"]
+    print("theta unscaled: ", theta_unscaled)
+    theta_original_log = xp.array(
+        [
+            theta_original[0],
+            xp.log(theta_original[1]),
+        ]
+    )
+    print("theta original log: ", theta_original_log)
+
+    print("norm(x_original - x): ", np.linalg.norm(x_original - results["x"]))
+    print("normalized norm: ", np.linalg.norm(x_original - results["x"]) / np.linalg.norm(x_original))
+
+    print("norm(eta_original - eta_est): ", np.linalg.norm(model.a @ x_original - model.a @ results["x"]))  
+    print(
+        "normalized norm: ",
+        np.linalg.norm(model.a @ x_original - model.a @ results["x"]) / np.linalg.norm(model.a @ x_original),
+    )
diff --git a/src/dalia/configs/dalia_config.py b/src/dalia/configs/dalia_config.py
index 0dbbbec0..8d5e24fc 100644
--- a/src/dalia/configs/dalia_config.py
+++ b/src/dalia/configs/dalia_config.py
@@ -55,6 +55,10 @@ class DaliaConfig(BaseModel):
     simulation_dir: Path = Path("./dalia/")
     output_dir: Path = Path.joinpath(simulation_dir, "output/")
 
+    # --- Verbosity level ------------------------------------------------------
+    verbosity: int = 0  # 0: minimal, 1: more info
+
+
 
 def parse_config(config: dict | str) -> DaliaConfig:
     if isinstance(config, str):
diff --git a/src/dalia/configs/priorhyperparameters_config.py b/src/dalia/configs/priorhyperparameters_config.py
index ddbea90c..7004ba8d 100644
--- a/src/dalia/configs/priorhyperparameters_config.py
+++ b/src/dalia/configs/priorhyperparameters_config.py
@@ -13,7 +13,9 @@
 class PriorHyperparametersConfig(BaseModel):
     model_config = ConfigDict(extra="forbid", arbitrary_types_allowed=True)
 
-    type: Literal["gaussian", "penalized_complexity", "beta", "gaussian_mvn"] = None
+    type: Literal[
+        "gaussian", "penalized_complexity", "beta", "gaussian_mvn", "gamma"
+    ] = None
 
 
 class GaussianPriorHyperparametersConfig(PriorHyperparametersConfig):
@@ -44,6 +46,11 @@ class BetaPriorHyperparametersConfig(PriorHyperparametersConfig):
     beta: float = None
 
 
+class GammaPriorHyperparametersConfig(PriorHyperparametersConfig):
+    alpha: float = None
+    beta: float = None
+
+
 def parse_config(config: dict) -> PriorHyperparametersConfig:
     prior_type = config.get("type")
     if prior_type == "gaussian":
@@ -54,5 +61,7 @@ def parse_config(config: dict) -> PriorHyperparametersConfig:
         return PenalizedComplexityPriorHyperparametersConfig(**config)
     elif prior_type == "beta":
         return BetaPriorHyperparametersConfig(**config)
+    elif prior_type == "gamma":
+        return GammaPriorHyperparametersConfig(**config)
     else:
         raise ValueError(f"Unknown prior hyperparameters config type: {prior_type}")
diff --git a/src/dalia/configs/submodels_config.py b/src/dalia/configs/submodels_config.py
index f31a9914..8101565b 100644
--- a/src/dalia/configs/submodels_config.py
+++ b/src/dalia/configs/submodels_config.py
@@ -24,7 +24,7 @@ class SubModelConfig(BaseModel, ABC):
 
     # Input folder for this specific submodel
     input_dir: str = None
-    type: Literal["spatio_temporal", "spatial", "regression", "brainiac"] = None
+    type: Literal["spatio_temporal", "spatial", "regression", "brainiac", "ar1"] = None
 
     @abstractmethod
     def read_hyperparameters(self) -> tuple[ArrayLike, list]:
@@ -39,6 +39,30 @@ def read_hyperparameters(self):
         return xp.array([]), []
 
 
+class AR1SubModelConfig(SubModelConfig):
+
+    ## prior on phi
+    phi: float = None  # AR(1) coefficient
+    phi_scaled: float = None
+    ph_phi: PriorHyperparametersConfig = None
+    ## check that phi is between -1 and 1 (use pc prior)
+    # check inla.doc("pc.cor1")
+
+    ## prior (in log-scale or not?) on tau (marginal precision) or s2 marginal variance
+    tau: float = None  # Marginal variance
+    ph_tau: PriorHyperparametersConfig = None
+
+    def read_hyperparameters(self):
+
+        # input of phi is in (0,1), rescale to -/+ INF
+        self.phi_scaled = scaled_logit(self.phi, direction="forward")
+        theta = xp.array([self.phi_scaled, self.tau])
+        theta_internal = xp.array([self.phi_scaled, self.tau])
+        theta_keys = ["phi", "tau"]
+
+        return theta, theta_keys
+
+
 class SpatioTemporalSubModelConfig(SubModelConfig):
     spatial_domain_dimension: PositiveInt = 2
 
@@ -122,6 +146,10 @@ def parse_config(config: dict | str) -> SubModelConfig:
         config["ph_h2"] = parse_priorhyperparameters_config(config["ph_h2"])
         config["ph_alpha"] = parse_priorhyperparameters_config(config["ph_alpha"])
         return BrainiacSubModelConfig(**config)
+    elif type == "ar1":
+        config["ph_tau"] = parse_priorhyperparameters_config(config["ph_tau"])
+        config["ph_phi"] = parse_priorhyperparameters_config(config["ph_phi"])
+        return AR1SubModelConfig(**config)
     # Add more elif branches for other submodel types
     else:
         raise ValueError(f"Unknown submodel type: {type}")
diff --git a/src/dalia/core/dalia.py b/src/dalia/core/dalia.py
index ac06436f..e502751a 100644
--- a/src/dalia/core/dalia.py
+++ b/src/dalia/core/dalia.py
@@ -26,6 +26,8 @@
     smartsplit,
     synchronize,
     synchronize_gpu,
+    check_vector_consistency,
+    compute_outer_covariance_matrix,
 )
 
 if backend_flags["mpi_avail"]:
@@ -73,6 +75,8 @@ def __init__(
         self.eps_gradient_f = self.config.eps_gradient_f
         self.eps_hessian_f = self.config.eps_hessian_f
 
+        self.verbosity = self.config.verbosity
+
         # --- Configure HPC
         set_device(comm_rank, comm_size)
 
@@ -198,14 +202,14 @@ def __init__(
             dtype=xp.float64,
         )
         self.theta_mat = xp.zeros(
-            (self.model.theta.size, self.n_f_evaluations), dtype=xp.float64
+            (self.model.theta_internal.size, self.n_f_evaluations), dtype=xp.float64
         )
-        self.theta_optimizer = xp.zeros_like(self.model.theta)
-        self.theta_optimizer[:] = self.model.theta
+        self.theta_optimizer = xp.zeros_like(self.model.theta_internal)
+        self.theta_optimizer[:] = self.model.theta_internal
 
         # --- Metrics
         self.f_values: ArrayLike = []
-        self.theta_values: ArrayLike = []
+        self.theta_values_internal: ArrayLike = []
         self.objective_function_time: ArrayLike = []
         self.solver_time: ArrayLike = []
         self.construction_time: ArrayLike = []
@@ -303,14 +307,34 @@ def run(self) -> dict:
 
         theta_star = get_device(minimization_result["theta"])
         x_star = get_device(minimization_result["x"])
+        print("Finished the optimization procedure.")
+
+        # need to update theta_star and x_star to be the same across all ranks
+        theta_star[:] = self.comm_world.bcast(theta_star, root=0)
+        x_star[:] = self.comm_world.bcast(x_star, root=0)
+
+        # need to update theta_star and x_star to be the same across all ranks
+        theta_star[:] = self.comm_world.bcast(theta_star, root=0)
+        x_star[:] = self.comm_world.bcast(x_star, root=0)
 
         # compute covariance of the hyperparameters theta at the mode
-        cov_theta = self.compute_covariance_hp(theta_star)
+        print("theta_star: ", theta_star)
+        cov_theta_dict = self.compute_covariance_hp(theta_star)
+        print("Computed covariance of the hyperparameters at the mode.")
+
+        # need to update theta_star and x_star to be the same across all ranks
+        theta_star[:] = self.comm_world.bcast(theta_star, root=0)
+        x_star[:] = self.comm_world.bcast(x_star, root=0)
+
+        # need to update theta_star and x_star to be the same across all ranks
+        theta_star[:] = self.comm_world.bcast(theta_star, root=0)
+        x_star[:] = self.comm_world.bcast(x_star, root=0)
 
         # compute marginal variances of the latent parameters
         marginal_variances_latent = self.get_marginal_variances_latent_parameters(
             theta_star, x_star
         )
+        print("Computed marginal variances of the latent parameters.")
 
         # compute marginal variances of the observations
         # TODO: only run by default when dense multiplcation issue is fixed, see issue #78
@@ -321,13 +345,14 @@ def run(self) -> dict:
         # construct new dictionary with the results
         results = {
             "theta": minimization_result["theta"],
-            "theta_interpret": minimization_result["theta_interpret"],
+            "theta_internal": minimization_result["theta_internal"],
             "x": minimization_result["x"],
             "f": minimization_result["f"],
             "grad_f": minimization_result["grad_f"],
             "f_values": minimization_result["f_values"],
             "theta_values": minimization_result["theta_values"],
-            "cov_theta": cov_theta,
+            "cov_theta_internal": cov_theta_dict["internal"],
+            "cov_theta": cov_theta_dict["external"],
             "marginal_variances_latent": marginal_variances_latent,
             # "marginal_variances_observations": get_host(
             #     marginal_variances_observations
@@ -349,13 +374,19 @@ def minimize(self) -> optimize.OptimizeResult:
             Result of the optimization procedure.
         """
 
-        if len(self.model.theta) == 0:
+        # ensure that all ranks are initialized to the same theta
+        check_vector_consistency(
+            self.model.theta_external,
+            comm=self.comm_world,
+        )
+
+        if len(self.model.theta_external) == 0:
             # Only run the inner iteration
             print_msg("No hyperparameters, just running inner iteration.")
-            self.f_value = self._evaluate_f(self.model.theta)
+            self.f_value = self._evaluate_f(self.model.theta_external)
             self.minimization_result: dict = {
-                "theta": self.model.theta,
-                "theta_interpret": self.model.get_theta_interpret(),
+                "theta_internal": self.model.theta_internal,
+                "theta": self.model.theta_external,
                 "x": self.model.x,  # [self.model.inverse_permutation_latent_variables],
                 "f": self.f_value,
             }
@@ -386,12 +417,12 @@ def callback(intermediate_result: optimize.OptimizeResult):
                     f"Iteration: {self.accepted_iter:2d} (took: {self.objective_function_time[-1]:.2f}) | "
                     f"Theta: [{theta_str}] | "
                     f"Function Value: {fun_i: .6f} | "
-                    f"Gradient: [{gradient_str}] | ",
+                    #f"Gradient: [{gradient_str}] | ",
                     f"Norm(Grad): [{xp.linalg.norm(self.gradient_f): .6f}]",
                     flush=True,
                 )
 
-                self.theta_values.append(theta_i)
+                self.theta_values_internal.append(theta_i)
                 self.f_values.append(fun_i)
 
                 # check if f_values have been decreasing over last iterations
@@ -409,9 +440,9 @@ def callback(intermediate_result: optimize.OptimizeResult):
                         )
 
                         self.minimization_result = {
-                            "theta": get_host(self.model.theta),
-                            "theta_interpret": get_host(
-                                self.model.get_theta_interpret()
+                            "theta_internal": get_host(self.model.theta_internal),
+                            "theta": get_host(
+                                self.model.theta_external
                             ),
                             "x": get_host(
                                 self.model.x
@@ -422,7 +453,8 @@ def callback(intermediate_result: optimize.OptimizeResult):
                             "f": fun_i,
                             "grad_f": self.gradient_f,
                             "f_values": self.f_values,
-                            "theta_values": self.theta_values,
+                            ### these values are in internal scale (!!)
+                            "theta_values": self.theta_values_internal,
                         }
 
                         raise OptimizationConvergedEarlyExit()
@@ -430,13 +462,13 @@ def callback(intermediate_result: optimize.OptimizeResult):
                 if self.accepted_iter > self.config.theta_reduction_lag:
                     if (
                         xp.linalg.norm(
-                            self.theta_values[-self.config.theta_reduction_lag]
+                            self.theta_values_internal[-self.config.theta_reduction_lag]
                             - theta_i
                         )
                         < self.config.theta_reduction_tol
                     ):
                         norm_diff = xp.linalg.norm(
-                            self.theta_values[
+                            self.theta_values_internal[
                                 self.accepted_iter - self.config.theta_reduction_lag
                             ]
                             - theta_i
@@ -450,10 +482,9 @@ def callback(intermediate_result: optimize.OptimizeResult):
                         )
 
                         self.minimization_result = {
-                            "theta": get_host(self.model.theta),
-                            "theta_interpret": get_host(
-                                self.model.get_theta_interpret()
-                            ),
+                            "theta_internal": get_host(self.model.theta_internal),
+                            "theta": get_host(
+                                self.model.theta_external),
                             "x": get_host(
                                 self.model.x
                                 # self.model.x[
@@ -463,7 +494,7 @@ def callback(intermediate_result: optimize.OptimizeResult):
                             "f": fun_i,
                             "grad_f": self.gradient_f,
                             "f_values": self.f_values,
-                            "theta_values": self.theta_values,
+                            "theta_values": self.theta_values_internal,
                         }
 
                         raise OptimizationConvergedEarlyExit()
@@ -513,15 +544,15 @@ def callback(intermediate_result: optimize.OptimizeResult):
                 )
 
             self.minimization_result: dict = {
-                "theta": scipy_result.x,
-                "theta_interpret": self.model.get_theta_interpret(),
+                "theta_internal": get_host(self.model.theta_internal), #  scipy_result.x, #
+                "theta": get_host(self.model.theta_external),
                 "x": get_host(
                     self.model.x,  # [self.model.inverse_permutation_latent_variables]
                 ),
                 "f": scipy_result.fun,
                 "grad_f": self.gradient_f,
                 "f_values": self.f_values,
-                "theta_values": self.theta_values,
+                "theta_values": self.theta_values_internal,
             }
 
         return self.minimization_result
@@ -607,7 +638,7 @@ def _objective_function(
             self.t_construction_qprior + self.t_construction_qconditional
         )
 
-        if self.iter > 0:
+        if self.iter > 0 and self.verbosity > 0:
             print(
                 f"rank {comm_rank} | objfunc_time: {self.objective_function_time[1:]} | solver_time: {self.solver_time[1:]} | construction_time: {self.construction_time[1:]}",
                 flush=True,
@@ -644,7 +675,8 @@ def _evaluate_f(
 
         tic = time.time()
 
-        self.model.theta[:] = theta_i
+        #self.model.theta_internal[:] = theta_i
+        self.model.theta_internal = theta_i
         f_theta = xp.zeros(1, dtype=xp.float64)
 
         # --- Optimize x and evaluate the conditional of the latent parameters
@@ -664,6 +696,7 @@ def _evaluate_f(
             task_mapping = [i % n_qeval_comm for i in range(2)]
 
             if task_mapping[0] == self.color_qeval:
+
                 # Done by processes "even"
                 tic = time.perf_counter()
                 synchronize_gpu()
@@ -698,6 +731,7 @@ def _evaluate_f(
                 log_prior_hyperparameters: float = (
                     self.model.evaluate_log_prior_hyperparameters()
                 )
+
                 likelihood: float = float(self.model.evaluate_likelihood(eta=eta))
                 prior_latent_parameters: float = (
                     self._evaluate_prior_latent_parameters()
@@ -716,6 +750,7 @@ def _evaluate_f(
                     comm=self.comm_feval,
                 )
                 synchronize(comm=self.comm_qeval)
+            
         else:
             self.model.construct_Q_prior()
 
@@ -764,7 +799,7 @@ def _evaluate_f(
 
         return f_theta[0]
 
-    def compute_covariance_hp(self, theta_i: NDArray) -> NDArray:
+    def compute_covariance_hp(self, theta_external: NDArray) -> NDArray:
         """compute the covariance matrix of the hyperparameters theta.
 
         Parameters
@@ -777,24 +812,45 @@ def compute_covariance_hp(self, theta_i: NDArray) -> NDArray:
         cov_theta : NDArray[dim_theta, dim_theta]
             Covariance matrix of the hyperparameters theta.
         """
-        self.model.theta[:] = theta_i
+
+        # ensure that all ranks are initialized to the same theta
+        check_vector_consistency(
+            theta_external,
+            comm=self.comm_world,
+        )
+
+        # self.model.rescale_hyperparameters_to_internal(theta_interpret, direction="forward")
+        # ensure that all ranks are initialized to the same theta
+        check_vector_consistency(
+            theta_external,
+            comm=self.comm_world,
+        )
         print_msg(
-            f"Computing covariance of hyperparameters theta at {theta_i}.",
+            f"Computing covariance of hyperparameters at theta_external {theta_external}.",
             flush=True,
         )
+        
+        self.model.theta_external = theta_external
 
-        hess_theta = self._evaluate_hessian_f(theta_i)
-        # print_msg(
-        #     f"hessian_f: \n {hess_theta}",
-        #     flush=True,
-        # )
-        cov_theta = xp.linalg.inv(hess_theta)
+        hess_theta_internal = self._evaluate_hessian_f(self.model.theta_internal)
+        print_msg(
+            f"hessian_f: \n {hess_theta_internal}",
+            flush=True,
+        )
+        self.cov_theta_internal = xp.linalg.inv(hess_theta_internal)
+        
+        # rescale to external scale
+        cov_theta_external = compute_outer_covariance_matrix(self.model.theta_internal, self.cov_theta_internal, self.model.rescale_hyperparameters_to_internal)
 
-        return cov_theta
+        dict_cov = {"internal": self.cov_theta_internal, "external": cov_theta_external}
+        print("Cov Internal: ", dict_cov["internal"])
+        print("Cov External: ", dict_cov["external"])
+
+        return dict_cov
 
     def _evaluate_hessian_f(
         self,
-        theta_i: NDArray,
+        theta_internal: NDArray,
     ) -> NDArray:
         """Approximate the hessian of the function f(theta) = log(p(theta|y)).
 
@@ -811,10 +867,9 @@ def _evaluate_hessian_f(
         -----
         Compute finite difference approximation of the hessian of f at theta_i.
         """
-
+        
         ## TODO: this is the quick fix ...
-        theta_internal = theta_i.copy()
-        self.model.theta[:] = theta_i
+        # self.model.theta[:] = theta_i
         dim_theta = self.model.n_hyperparameters
 
         # pre-allocate storage for the hessian & f_values
@@ -846,8 +901,10 @@ def _evaluate_hessian_f(
         counter = 0
         # compute f(theta)
         if self.color_feval == task_mapping[0]:
+            theta_i = theta_internal.copy()
             f_theta = self._evaluate_f(theta_i)
             f_ii_loc[1, :] = f_theta
+
         counter += 1
 
         for k in range(loop_dim):
@@ -858,49 +915,62 @@ def _evaluate_hessian_f(
             if i == j:
                 if self.color_feval == task_mapping[counter]:
                     # theta+eps_i
-                    theta_i = theta_internal.copy()
-                    f_ii_loc[0, i] = self._evaluate_f(theta_i + eps_mat[i, :])
+                    #theta_i = theta_internal.copy()
+                    #f_ii_loc[0, i] = self._evaluate_f(theta_i + eps_mat[i, :])
+                    theta_i = theta_internal + eps_mat[i, :]
+                    result = self._evaluate_f(theta_i)
+                    f_ii_loc[0, i] = result
                 counter += 1
 
                 if self.color_feval == task_mapping[counter]:
                     # theta-eps_i
-                    theta_i = theta_internal.copy()
-                    f_ii_loc[2, i] = self._evaluate_f(theta_i - eps_mat[i, :])
-
+                    # theta_i = theta_internal.copy()
+                    # f_ii_loc[2, i] = self._evaluate_f(theta_i - eps_mat[i, :])
+                    theta_i = theta_internal - eps_mat[i, :]
+                    result = self._evaluate_f(theta_i)
+                    f_ii_loc[2, i] = result
                 counter += 1
 
             # as hessian is symmetric we only have to compute the upper triangle
             elif i < j:
                 # theta+eps_i+eps_j
                 if self.color_feval == task_mapping[counter]:
-                    theta_i = theta_internal.copy()
-                    f_ij_loc[0, k] = self._evaluate_f(
-                        theta_i + eps_mat[i, :] + eps_mat[j, :]
-                    )
+                    # theta_i = theta_internal.copy()
+                    # f_ij_loc[0, k] = self._evaluate_f(
+                    #     theta_i + eps_mat[i, :] + eps_mat[j, :]
+                    # )
+                    theta_i = theta_internal + eps_mat[i, :] + eps_mat[j, :]
+                    f_ij_loc[0, k] = self._evaluate_f(theta_i)
                 counter += 1
 
                 # theta+eps_i-eps_j
                 if self.color_feval == task_mapping[counter]:
-                    theta_i = theta_internal.copy()
-                    f_ij_loc[1, k] = self._evaluate_f(
-                        theta_i + eps_mat[i, :] - eps_mat[j, :]
-                    )
+                    # theta_i = theta_internal.copy()
+                    # f_ij_loc[1, k] = self._evaluate_f(
+                    #     theta_i + eps_mat[i, :] - eps_mat[j, :]
+                    # )
+                    theta_i = theta_internal + eps_mat[i, :] - eps_mat[j, :]
+                    f_ij_loc[1, k] = self._evaluate_f(theta_i)
                 counter += 1
 
                 # theta-eps_i+eps_j
                 if self.color_feval == task_mapping[counter]:
-                    theta_i = theta_internal.copy()
-                    f_ij_loc[2, k] = self._evaluate_f(
-                        theta_i - eps_mat[i, :] + eps_mat[j, :]
-                    )
+                    # theta_i = theta_internal.copy()
+                    # f_ij_loc[2, k] = self._evaluate_f(
+                    #     theta_i - eps_mat[i, :] + eps_mat[j, :]
+                    # )
+                    theta_i = theta_internal - eps_mat[i, :] + eps_mat[j, :]
+                    f_ij_loc[2, k] = self._evaluate_f(theta_i)
                 counter += 1
 
                 # theta-eps_i-eps_j
                 if self.color_feval == task_mapping[counter]:
-                    theta_i = theta_internal.copy()
-                    f_ij_loc[3, k] = self._evaluate_f(
-                        theta_i - eps_mat[i, :] - eps_mat[j, :]
-                    )
+                    # theta_i = theta_internal.copy()
+                    # f_ij_loc[3, k] = self._evaluate_f(
+                    #     theta_i - eps_mat[i, :] - eps_mat[j, :]
+                    # )
+                    theta_i = theta_internal - eps_mat[i, :] - eps_mat[j, :]
+                    f_ij_loc[3, k] = self._evaluate_f(theta_i)
                 counter += 1
 
         allreduce(
@@ -943,23 +1013,24 @@ def _evaluate_hessian_f(
         return hess
 
     def _compute_covariance_latent_parameters(
-        self, theta: NDArray, x_star: NDArray
+        self, theta_external: NDArray, x_star: NDArray
     ) -> None:
         """Compute the marginal distribution of the latent parameters x.
 
         Parameters
         ----------
-        theta_i : NDArray
+        theta_internal : NDArray
             Hyperparameters theta.
         x_star : NDArray
-            Latent parameters x(theta_i).
+            Latent parameters x(theta_internal).
 
         Returns
         -------
         marginal_latent_parameters : NDArray
             Marginal distribution of the latent parameters x.
         """
-        self.model.theta[:] = theta
+        print("Computing covariance of latent parameters at theta external:", theta_external)
+        self.model.theta_external = xp.array(theta_external)
         self.model.x[:] = x_star
 
         eta = self.model.a @ self.model.x
@@ -969,19 +1040,28 @@ def _compute_covariance_latent_parameters(
         self.solver.selected_inversion(sparsity="bta")
 
     def get_marginal_variances_latent_parameters(
-        self, theta: NDArray = None, x_star: NDArray = None
+        self, theta_external: NDArray = None, x_star: NDArray = None
     ) -> NDArray:
+
         # TODO: this should be only called by rank 0?
-        if theta is None and x_star is None:
+        if theta_external is None and x_star is None:
             print(
                 "Computing marginal variances for currently stored latent parameters. "
             )
             x_star = self.model.x
-            theta = self.model.theta
+            theta = self.model.theta_internal
+        elif theta_external is not None and x_star is not None:
+            ## assume theta to be in "external" scale
+            self.model.theta_external = xp.atleast_1d(theta_external)
+            theta = self.model.theta_internal
+
         elif theta is None or x_star is None:
             raise ValueError(
                 "BOTH or NEITHER theta and x_star must be provided to compute the marginal variances."
             )
+        
+        check_vector_consistency(theta, comm=self.comm_world)
+        check_vector_consistency(x_star, comm=self.comm_world)
 
         # check order x_star ... -> potentially need to reorder marginal variances
         self._compute_covariance_latent_parameters(theta, x_star)
@@ -995,7 +1075,7 @@ def get_marginal_variances_latent_parameters(
         return marginal_variances
 
     def get_marginal_variances_observations(
-        self, theta: NDArray, x_star: NDArray
+        self, theta_external: NDArray = None, x_star: NDArray = None
     ) -> NDArray:
         """Extract the marginal variances of the observations.
 
@@ -1018,21 +1098,27 @@ def get_marginal_variances_observations(
             Marginal variances of the observations.
         """
 
+        # TODO: implement this for non-Gaussian likelihoods
+        check_vector_consistency(theta_external, comm=self.comm_world)
+        check_vector_consistency(x_star, comm=self.comm_world)
+
         if self.model.is_likelihood_gaussian():
             # TODO: this should be only called by rank 0?
-            if theta is None and x_star is None:
+            if theta_external is None and x_star is None:
                 print(
                     "Computing marginal variances for currently stored latent parameters. "
                 )
                 x_star = self.model.x
-                theta = self.model.theta
-            elif theta is None or x_star is None:
+                theta_external = self.model.theta_external
+            
+            if theta_external is None or x_star is None:
                 raise ValueError(
                     "BOTH or NEITHER theta and x_star must be provided to compute the marginal variances."
                 )
 
+
                 # check order x_star ... -> potentially need to reorder marginal variances
-            self._compute_covariance_latent_parameters(theta, x_star)
+            self._compute_covariance_latent_parameters(theta_external, x_star)
 
             # now only extract diagonal elements corresponding to marginal variances of the latent parameters
             variances_latent = self.solver._structured_to_spmatrix(
@@ -1077,7 +1163,7 @@ def _inner_iteration(
             if counter > self.inner_iteration_max_iter:
                 print_msg(
                     "Theta value at failing of the inner_iteration: ",
-                    self.model.theta,
+                    self.model.theta_internal,
                     flush=True,
                 )
                 raise ValueError(
diff --git a/src/dalia/core/model.py b/src/dalia/core/model.py
index 91608d62..ab927f31 100644
--- a/src/dalia/core/model.py
+++ b/src/dalia/core/model.py
@@ -14,6 +14,7 @@
     GaussianMVNPriorHyperparametersConfig,
     GaussianPriorHyperparametersConfig,
     PenalizedComplexityPriorHyperparametersConfig,
+    GammaPriorHyperparametersConfig,
 )
 from dalia.core.likelihood import Likelihood
 from dalia.core.prior_hyperparameters import PriorHyperparameters
@@ -24,12 +25,14 @@
     GaussianMVNPriorHyperparameters,
     GaussianPriorHyperparameters,
     PenalizedComplexityPriorHyperparameters,
+    GammaPriorHyperparameters,
 )
 from dalia.submodels import (
     BrainiacSubModel,
     RegressionSubModel,
     SpatialSubModel,
     SpatioTemporalSubModel,
+    AR1SubModel,
 )
 from dalia.utils import add_str_header, boxify, scaled_logit
 
@@ -56,8 +59,11 @@ def __init__(
 
         self.n_fixed_effects: int = 0
 
+        self._theta_external: ArrayLike = []
+        self._theta_internal: ArrayLike = []
+
         # For each submodel...
-        theta: ArrayLike = []
+        theta_external: ArrayLike = []
         theta_keys: ArrayLike = []
         self.hyperparameters_idx: ArrayLike = [0]
         self.prior_hyperparameters: list[PriorHyperparameters] = []
@@ -154,6 +160,36 @@ def __init__(
             elif isinstance(submodel, RegressionSubModel):
                 self.n_fixed_effects += submodel.n_fixed_effects
 
+            elif isinstance(submodel, AR1SubModel):
+
+                if isinstance(submodel.config.ph_phi, BetaPriorHyperparametersConfig):
+                    self.prior_hyperparameters.append(
+                        BetaPriorHyperparameters(
+                            config=submodel.config.ph_phi,
+                        )
+                    )
+                elif isinstance(
+                    submodel.config.ph_phi,
+                    PenalizedComplexityPriorHyperparametersConfig,
+                ):
+                    self.prior_hyperparameters.append(
+                        PenalizedComplexityPriorHyperparameters(
+                            config=submodel.config.ph_phi,
+                            hyperparameter_type="phi",
+                        )
+                    )
+
+                if isinstance(
+                    submodel.config.ph_tau, GaussianPriorHyperparametersConfig
+                ):
+                    self.prior_hyperparameters.append(
+                        GaussianPriorHyperparameters(
+                            config=submodel.config.ph_tau,
+                        )
+                    )
+                else:
+                    raise ValueError("Unknown prior hyperparameter type for ph_tau")
+
             elif isinstance(submodel, BrainiacSubModel):
                 # h2 hyperparameters
                 if isinstance(submodel.config.ph_h2, BetaPriorHyperparametersConfig):
@@ -172,33 +208,28 @@ def __init__(
                             config=submodel.config.ph_alpha,
                         )
                     )
+                if isinstance(
+                    submodel.config.ph_alpha,
+                ):
+                    self.prior_hyperparameters.append(
+                        PenalizedComplexityPriorHyperparameters(
+                            config=submodel.config.ph_alpha,
+                            hyperparameter_type="alpha",
+                        )
+                    )
             else:
                 raise ValueError("Unknown submodel type")
 
             # ...and read their hyperparameters
             theta_submodel, theta_keys_submodel = submodel.config.read_hyperparameters()
 
-            theta.append(theta_submodel)
+            theta_external.append(theta_submodel)
             theta_keys += theta_keys_submodel
 
             self.hyperparameters_idx.append(
                 self.hyperparameters_idx[-1] + len(theta_submodel)
             )
 
-        # Add the likelihood hyperparameters
-        (
-            lh_hyperparameters,
-            lh_hyperparameters_keys,
-        ) = likelihood_config.read_hyperparameters()
-
-        theta.append(lh_hyperparameters)
-        self.theta: NDArray = xp.concatenate(theta)
-
-        theta_keys += lh_hyperparameters_keys
-        self.theta_keys: NDArray = theta_keys
-
-        self.n_hyperparameters = self.theta.size
-
         # --- Initialize the latent parameters and the design matrix
         self.n_latent_parameters: int = 0
         self.latent_parameters_idx: list[int] = [0]
@@ -305,6 +336,24 @@ def __init__(
                         hyperparameter_type="prec_o",
                     )
                 )
+            elif isinstance(
+                likelihood_config.prior_hyperparameters,
+                BetaPriorHyperparametersConfig,
+            ):
+                self.prior_hyperparameters.append(
+                    BetaPriorHyperparameters(
+                        config=likelihood_config.prior_hyperparameters,
+                    )
+                )
+            elif isinstance(
+                likelihood_config.prior_hyperparameters,
+                GammaPriorHyperparametersConfig,
+            ):
+                self.prior_hyperparameters.append(
+                    GammaPriorHyperparameters(
+                        config=likelihood_config.prior_hyperparameters,
+                    )
+                )
         elif likelihood_config.type == "poisson":
             self.likelihood: Likelihood = PoissonLikelihood(
                 n_observations=self.n_observations,
@@ -317,12 +366,60 @@ def __init__(
             )
 
         self.likelihood_config: LikelihoodConfig = likelihood_config
+        
+        # Add the likelihood hyperparameters
+        (
+            lh_hyperparameters,
+            lh_hyperparameters_keys,
+        ) = likelihood_config.read_hyperparameters()
+
+        theta_external.append(lh_hyperparameters)
+        self.theta_external = np.concatenate(theta_external)
+        # self.theta = np.concatenate(theta_external)
+
+        print("Initial hyperparameters (external scale): ", self.theta_external)
+        print("Initial hyperparameters (internal scale): ", self.theta_internal)
+
+        theta_keys += lh_hyperparameters_keys
+        self.theta_keys: NDArray = theta_keys
+
+        self.n_hyperparameters = self.theta_external.size
 
         # --- Recurrent variables
         self.Q_prior = None
         self.Q_prior_data_mapping = [0]
         self.Q_conditional = None
         self.Q_conditional_data_mapping = [0]
+     
+    ########################################################################       
+    @property
+    def theta_external(self):
+        """External/user/interpretable scale theta."""
+        # the copy is important to make sure that in place operations still trigger updating
+        return self._theta_external.copy()
+    
+    @theta_external.setter
+    def theta_external(self, value):
+        """Set external theta and automatically update internal."""
+        
+        self._theta_external = np.array(value) 
+        self._theta_internal = self.rescale_hyperparameters_to_internal(
+            self._theta_external, direction="forward"
+        )
+    
+    @property
+    def theta_internal(self):
+        """Internal/BFGS scale theta."""
+        return self._theta_internal.copy()
+    
+    @theta_internal.setter
+    def theta_internal(self, value):
+        """Set internal theta and automatically update external."""
+        self._theta_internal = np.array(value) 
+        self._theta_external = self.rescale_hyperparameters_to_internal(
+                self._theta_internal, direction="backward"
+        )
+    ########################################################################
 
     def construct_Q_prior(self) -> sp.sparse.spmatrix:
         kwargs = {}
@@ -334,22 +431,32 @@ def construct_Q_prior(self) -> sp.sparse.spmatrix:
             cols = []
             data = []
 
+            ## TODO: improve the if / elif statements
             for i, submodel in enumerate(self.submodels):
                 if isinstance(submodel, SpatioTemporalSubModel):
                     for hp_idx in range(
                         self.hyperparameters_idx[i], self.hyperparameters_idx[i + 1]
                     ):
-                        kwargs[self.theta_keys[hp_idx]] = float(self.theta[hp_idx])
+                        kwargs[self.theta_keys[hp_idx]] = float(self.theta_external[hp_idx])
+                        #kwargs[self.theta_keys[hp_idx]] = float(theta_interpret[hp_idx])
                 elif isinstance(submodel, SpatialSubModel):
                     for hp_idx in range(
                         self.hyperparameters_idx[i], self.hyperparameters_idx[i + 1]
                     ):
-                        kwargs[self.theta_keys[hp_idx]] = float(self.theta[hp_idx])
+                        kwargs[self.theta_keys[hp_idx]] = float(self.theta_external[hp_idx])
+                        #kwargs[self.theta_keys[hp_idx]] = float(theta_interpret[hp_idx])
                 elif isinstance(submodel, BrainiacSubModel):
                     for hp_idx in range(
                         self.hyperparameters_idx[i], self.hyperparameters_idx[i + 1]
                     ):
-                        kwargs[self.theta_keys[hp_idx]] = float(self.theta[hp_idx])
+                        kwargs[self.theta_keys[hp_idx]] = float(self.theta_external[hp_idx])
+                        #kwargs[self.theta_keys[hp_idx]] = float(theta_interpret[hp_idx])
+                elif isinstance(submodel, AR1SubModel):
+                    for hp_idx in range(
+                        self.hyperparameters_idx[i], self.hyperparameters_idx[i + 1]
+                    ):
+                        kwargs[self.theta_keys[hp_idx]] = float(self.theta_external[hp_idx])
+                        #kwargs[self.theta_keys[hp_idx]] = float(theta_interpret[hp_idx])
                 elif isinstance(submodel, RegressionSubModel):
                     ...
 
@@ -382,17 +489,26 @@ def construct_Q_prior(self) -> sp.sparse.spmatrix:
                     for hp_idx in range(
                         self.hyperparameters_idx[i], self.hyperparameters_idx[i + 1]
                     ):
-                        kwargs[self.theta_keys[hp_idx]] = float(self.theta[hp_idx])
+                        kwargs[self.theta_keys[hp_idx]] = float(self.theta_external[hp_idx])
+                        #kwargs[self.theta_keys[hp_idx]] = float(theta_interpret[hp_idx])
                 elif isinstance(submodel, SpatialSubModel):
                     for hp_idx in range(
                         self.hyperparameters_idx[i], self.hyperparameters_idx[i + 1]
                     ):
-                        kwargs[self.theta_keys[hp_idx]] = float(self.theta[hp_idx])
+                        kwargs[self.theta_keys[hp_idx]] = float(self.theta_external[hp_idx])
+                        #kwargs[self.theta_keys[hp_idx]] = float(theta_interpret[hp_idx])
                 elif isinstance(submodel, BrainiacSubModel):
                     for hp_idx in range(
                         self.hyperparameters_idx[i], self.hyperparameters_idx[i + 1]
                     ):
-                        kwargs[self.theta_keys[hp_idx]] = float(self.theta[hp_idx])
+                        kwargs[self.theta_keys[hp_idx]] = float(self.theta_external[hp_idx])
+                        #kwargs[self.theta_keys[hp_idx]] = float(theta_interpret[hp_idx])
+                elif isinstance(submodel, AR1SubModel):
+                    for hp_idx in range(
+                        self.hyperparameters_idx[i], self.hyperparameters_idx[i + 1]
+                    ):
+                        kwargs[self.theta_keys[hp_idx]] = float(self.theta_external[hp_idx])
+                        #kwargs[self.theta_keys[hp_idx]] = float(theta_interpret[hp_idx])
 
                 submodel_Q_prior = submodel.construct_Q_prior(**kwargs)
 
@@ -418,7 +534,7 @@ def construct_Q_conditional(
         if self.likelihood_config.type == "gaussian":
             kwargs = {
                 "eta": eta,
-                "theta": float(self.theta[-1]),
+                "theta": float(self.theta_external[-1]),
             }
         else:
             kwargs = {
@@ -427,7 +543,7 @@ def construct_Q_conditional(
 
         if isinstance(self.submodels[0], BrainiacSubModel):
             # Brainiac specific rule
-            kwargs["h2"] = float(self.theta[0])
+            kwargs["h2"] = float(self.theta_external[0])
             d_matrix = self.submodels[0].evaluate_d_matrix(**kwargs)
         else:
             # General rules
@@ -461,7 +577,7 @@ def construct_information_vector(
         """Construct the information vector."""
 
         if isinstance(self.submodels[0], BrainiacSubModel):
-            kwargs = {"h2": float(self.theta[0])}
+            kwargs = {"h2": float(self.theta_external[0])}
             gradient_likelihood = self.submodels[0].evaluate_gradient_likelihood(
                 eta=eta, y=self.y, **kwargs
             )
@@ -470,7 +586,7 @@ def construct_information_vector(
             gradient_likelihood = self.likelihood.evaluate_gradient_likelihood(
                 eta=eta,
                 y=self.y,
-                theta=self.theta[self.hyperparameters_idx[-1] :],
+                theta=self.theta_external[self.hyperparameters_idx[-1] :],
             )
 
         information_vector: NDArray = (
@@ -487,23 +603,38 @@ def evaluate_log_prior_hyperparameters(self) -> float:
         """Evaluate the log prior hyperparameters."""
         log_prior = 0.0
 
-        # if BFGS and model scale differ: rescale -- generalize
-        if isinstance(self.submodels[0], BrainiacSubModel):
-            #
-            theta_interpret = self.theta.copy()
-            theta_interpret[0] = scaled_logit(self.theta[0], direction="backward")
-            log_prior += self.prior_hyperparameters[0].evaluate_log_prior(
-                theta_interpret[0]
-            )
+        #theta_interpret = self.theta.copy()
 
-            log_prior += self.prior_hyperparameters[1].evaluate_log_prior(
-                theta_interpret[1:]
-            )
-        else:
-            theta_interpret = self.theta
+        for i, prior_hyperparameter in enumerate(self.prior_hyperparameters):
+
+            # if isinstance(prior_hyperparameter, BetaPriorHyperparameters):
+            #     theta_interpret[i] = scaled_logit(
+            #         theta_interpret[i], direction="backward"
+            #     )
+            # theta_interpret[i] = prior_hyperparameter.rescale_hyperparameters_to_internal(
+            #     theta_interpret[i], direction="backward"
+            # )
+
+            log_prior += prior_hyperparameter.evaluate_log_prior(self.theta_external[i])
+            #log_prior += prior_hyperparameter.evaluate_log_prior(theta_interpret[i])
 
-            for i, prior_hyperparameter in enumerate(self.prior_hyperparameters):
-                log_prior += prior_hyperparameter.evaluate_log_prior(theta_interpret[i])
+        # if BFGS and model scale differ: rescale -- generalize
+        # if isinstance(self.submodels[0], BrainiacSubModel):
+        #     #
+        #     theta_interpret = self.theta.copy()
+        #     theta_interpret[0] = scaled_logit(self.theta[0], direction="backward")
+        #     log_prior += self.prior_hyperparameters[0].evaluate_log_prior(
+        #         theta_interpret[0]
+        #     )
+
+        #     log_prior += self.prior_hyperparameters[1].evaluate_log_prior(
+        #         theta_interpret[1:]
+        #     )
+        # else:
+        #     theta_interpret = self.theta
+
+        #     for i, prior_hyperparameter in enumerate(self.prior_hyperparameters):
+        #         log_prior += prior_hyperparameter.evaluate_log_prior(theta_interpret[i])
 
         return log_prior
 
@@ -511,13 +642,13 @@ def get_theta_likelihood(self) -> NDArray:
         """Return the likelihood hyperparameters."""
 
         if isinstance(self.submodels[0], BrainiacSubModel):
-            theta_likelihood = 1 - scaled_logit(self.theta[0], direction="backward")
+            theta_likelihood = 1 - self.theta_external[0]
         else:
-            theta_likelihood = self.theta[self.hyperparameters_idx[-1] :]
+            theta_likelihood = self.theta_external[self.hyperparameters_idx[-1] :]
 
         return theta_likelihood
 
-    def get_theta_interpret(self) -> NDArray:
+    def get_theta_interpret(self, direction: str) -> NDArray:
         theta_interpret = xp.zeros_like(self.theta)
 
         for i, submodel in enumerate(self.submodels):
@@ -528,20 +659,35 @@ def get_theta_interpret(self) -> NDArray:
             ] = submodel.rescale_hyperparameters_to_interpret(
                 self.theta[
                     self.hyperparameters_idx[i] : self.hyperparameters_idx[i + 1]
-                ]
+                ], direction=direction
             )
 
         return theta_interpret
+    
+    def rescale_hyperparameters_to_internal(self, theta, direction):
+
+        # need to iterate over theta and its prior hyperparameters
+        theta_internal = xp.copy(theta)
+
+        for i, prior_hyperparameter in enumerate(self.prior_hyperparameters):
+            ## how to handle priors that have multiple hyperparameters?
+            if isinstance(prior_hyperparameter, GaussianMVNPriorHyperparameters):
+                pass  # no rescaling implemented
+            else:
+                theta_internal[i] = prior_hyperparameter.rescale_hyperparameters_to_internal(theta[i], direction=direction)
+
+        return theta_internal
 
     def evaluate_likelihood(self, eta: NDArray, **kwargs) -> float:
         """Evaluate the likelihood."""
 
         if isinstance(self.submodels[0], BrainiacSubModel):
-            kwargs["h2"] = float(self.theta[0])
+            #kwargs["h2"] = float(self.theta[0])
+            kwargs["h2"] = float(self.theta_external[0])
             likelihood = self.submodels[0].evaluate_likelihood(eta, self.y, **kwargs)
         else:
             likelihood = self.likelihood.evaluate_likelihood(
-                eta, self.y, theta=self.theta[self.hyperparameters_idx[-1] :]
+                eta, self.y, theta=self.theta_external[self.hyperparameters_idx[-1] :]
             )
 
         return likelihood
diff --git a/src/dalia/core/prior_hyperparameters.py b/src/dalia/core/prior_hyperparameters.py
index 88c1ba4e..bd067fb5 100644
--- a/src/dalia/core/prior_hyperparameters.py
+++ b/src/dalia/core/prior_hyperparameters.py
@@ -16,6 +16,20 @@ def __init__(
 
         self.config: PriorHyperparametersConfig = config
 
+    @abstractmethod
+    def rescale_hyperparameters_to_internal(self, theta: float, direction: str) -> float:
+        """Rescale hyperparameters to and from internal scale.
+
+        Args:
+            theta: Hyperparameter
+            direction: "forward" or "backward"
+        Returns:
+            Rescaled hyperparameter.
+
+        """
+        
+        return theta
+
     @abstractmethod
     def evaluate_log_prior(self, theta: float) -> float:
         """Evaluate the log prior hyperparameters."""
diff --git a/src/dalia/core/submodel.py b/src/dalia/core/submodel.py
index 7cefe5e3..0c8123da 100644
--- a/src/dalia/core/submodel.py
+++ b/src/dalia/core/submodel.py
@@ -52,13 +52,6 @@ def __init__(
         except FileNotFoundError:
             self.x_initial: NDArray = xp.zeros((self.a.shape[1]), dtype=float)
 
-    def rescale_hyperparameters_to_interpret(self, theta: NDArray) -> NDArray:
-        """Rescale hyperparameters to interpret them.  Does nothing unless implemented in specific submodel.
-
-        Note: It doesnt include the hyperparameters from the likelihood. If they need rescaling too. Needs to be done in likelihood.
-        """
-
-        return theta
 
     @abstractmethod
     def construct_Q_prior(self, **kwargs) -> sp.sparse.coo_matrix:
diff --git a/src/dalia/likelihoods/gaussian.py b/src/dalia/likelihoods/gaussian.py
index ee69bd16..6ed8a864 100644
--- a/src/dalia/likelihoods/gaussian.py
+++ b/src/dalia/likelihoods/gaussian.py
@@ -29,7 +29,7 @@ def evaluate_likelihood(
 
         Evaluate Gaussian log-likelihood for a given set of observations, latent parameters, and design matrix, where
         the observations are assumed to be identically and independently distributed given eta (=A*x). Leading to:
-        log (p(y|eta)) = -0.5 * n * log(2 * pi) - 0.5 * n * theta_observations - 0.5 * exp(theta_observations) * (y - eta)^T * (y - eta)
+        log (p(y|eta)) = -0.5 * n * log(2 * pi) - 0.5 * n * log(theta) - 0.5 * theta * (y - eta)^T * (y - eta)
         where the constant in front of the likelihood is omitted.
 
         Parameters
@@ -40,7 +40,7 @@ def evaluate_likelihood(
             Vector of the observations.
         kwargs :
             theta : float
-                Specific parameter for the likelihood calculation.
+                precision parameter for the likelihood calculation. theta > 0.
 
         Returns
         -------
@@ -52,11 +52,16 @@ def evaluate_likelihood(
         if theta is None:
             raise ValueError("theta must be provided to evaluate gaussian likelihood.")
 
+        if theta <= 0:
+            raise ValueError(f"theta must be positive, got {theta}")
+
         yEta = eta - y
         # print("xp.exp(theta) in lh:", xp.exp(theta))
 
         likelihood: float = (
-            0.5 * theta * self.n_observations - 0.5 * xp.exp(theta) * yEta.T @ yEta
+            -0.5 * self.n_observations * xp.log(2 * xp.pi)
+            + 0.5 * xp.log(theta) * self.n_observations
+            - 0.5 * theta * yEta.T @ yEta
         )
 
         return likelihood
@@ -77,7 +82,7 @@ def evaluate_gradient_likelihood(
             Vector of the observations.
         kwargs :
             theta : float
-                Specific parameter for the likelihood calculation.
+                precision parameter for the likelihood calculation. theta > 0.
 
         Returns
         -------
@@ -91,7 +96,10 @@ def evaluate_gradient_likelihood(
                 "theta must be provided to evaluate gradient of gaussian likelihood."
             )
 
-        gradient_likelihood: NDArray = -xp.exp(theta) * (eta - y)
+        if theta <= 0:
+            raise ValueError(f"theta must be positive, got {theta}")
+
+        gradient_likelihood: NDArray = -theta * (eta - y)
 
         return gradient_likelihood
 
@@ -121,11 +129,11 @@ def evaluate_hessian_likelihood(
             raise ValueError(
                 "theta must be provided to evaluate gradient of gaussian likelihood."
             )
+        if theta <= 0:
+            raise ValueError(f"theta must be positive, got {theta}")
 
         # print("hessian lh: xp.exp(theta)", xp.exp(theta))
 
-        hessian_likelihood: ArrayLike = -xp.exp(theta) * sp.sparse.eye(
-            self.n_observations
-        )
+        hessian_likelihood: ArrayLike = -theta * sp.sparse.eye(self.n_observations)
 
         return hessian_likelihood
diff --git a/src/dalia/prior_hyperparameters/__init__.py b/src/dalia/prior_hyperparameters/__init__.py
index b561764e..f1b2e9d7 100644
--- a/src/dalia/prior_hyperparameters/__init__.py
+++ b/src/dalia/prior_hyperparameters/__init__.py
@@ -6,10 +6,12 @@
     PenalizedComplexityPriorHyperparameters,
 )
 from dalia.prior_hyperparameters.beta import BetaPriorHyperparameters
+from dalia.prior_hyperparameters.gamma import GammaPriorHyperparameters
 
 __all__ = [
     "GaussianPriorHyperparameters",
     "GaussianMVNPriorHyperparameters",
     "PenalizedComplexityPriorHyperparameters",
     "BetaPriorHyperparameters",
+    "GammaPriorHyperparameters",
 ]
diff --git a/src/dalia/prior_hyperparameters/beta.py b/src/dalia/prior_hyperparameters/beta.py
index aafa91e9..6d43373b 100644
--- a/src/dalia/prior_hyperparameters/beta.py
+++ b/src/dalia/prior_hyperparameters/beta.py
@@ -6,6 +6,7 @@
     GaussianPriorHyperparametersConfig,
 )
 from dalia.core.prior_hyperparameters import PriorHyperparameters
+from dalia.utils.link_functions import scaled_logit
 
 
 class BetaPriorHyperparameters(PriorHyperparameters):
@@ -21,6 +22,20 @@ def __init__(
         self.alpha: float = config.alpha
         self.beta: float = config.beta
 
+
+    def rescale_hyperparameters_to_internal(self, theta, direction):
+
+        ### TODO: on longer term make scaled_logit default but let it be configurable in config
+        ## beta prior is defined on [0,1], while BFGS works on (-inf, inf)
+        if direction == "forward":
+            theta_scaled = scaled_logit(theta, direction="forward")
+        elif direction == "backward":
+            theta_scaled = scaled_logit(theta, direction="backward")
+        else:
+            raise ValueError(f"Unknown direction: {direction}")
+
+        return theta_scaled
+
     def evaluate_log_prior(self, theta: float, **kwargs) -> float:
         """Evaluate the log prior hyperparameters."""
 
diff --git a/src/dalia/prior_hyperparameters/gamma.py b/src/dalia/prior_hyperparameters/gamma.py
new file mode 100644
index 00000000..65757091
--- /dev/null
+++ b/src/dalia/prior_hyperparameters/gamma.py
@@ -0,0 +1,321 @@
+# Copyright 2024-2025 DALIA authors. All rights reserved.
+from dalia import NDArray
+from scipy.sparse import spmatrix
+from dalia import sp, xp
+
+import numpy as np
+
+from dalia.configs.priorhyperparameters_config import (
+    GammaPriorHyperparametersConfig,
+)
+from dalia.core.prior_hyperparameters import PriorHyperparameters
+
+
+class GammaPriorHyperparameters(PriorHyperparameters):
+    """Gamma prior hyperparameters.
+
+    p(theta) = (beta^alpha / Gamma(alpha)) * theta^(alpha - 1) * exp(-beta * theta)
+
+    and in log scale:
+    log p(theta) = alpha * log(beta) - log(Gamma(alpha)) + (alpha - 1) * log(theta) - beta * theta
+
+    where theta is typically a positive parameter such as a precision or rate.
+
+    Parameters
+    ----------
+    config : GammaPriorHyperparametersConfig
+        Configuration object containing alpha and beta parameters.
+
+    Attributes
+    ----------
+    alpha : float. alpha > 0
+        Shape parameter of the Gamma distribution.
+    beta : float. beta > 0
+        Rate parameter of the Gamma distribution.
+    normalizing_constant : float
+        Precomputed normalizing constant for log probability evaluation.
+    """
+
+    def __init__(
+        self,
+        config: GammaPriorHyperparametersConfig,
+    ) -> None:
+        """
+        Initialize the Gamma prior hyperparameters.
+
+        Parameters
+        ----------
+        config : GammaPriorHyperparametersConfig
+            Configuration containing alpha (shape) and beta (rate) parameters.
+
+        Raises
+        ------
+        ValueError
+            If alpha or beta are not positive.
+        """
+        super().__init__(config)
+
+        self.alpha: float = config.alpha
+        self.beta: float = config.beta
+
+        # Validate alpha and beta are positive
+        if self.alpha <= 0:
+            raise ValueError(f"Alpha must be positive, got {self.alpha}")
+        if self.beta <= 0:
+            raise ValueError(f"Beta must be positive, got {self.beta}")
+
+        self.normalizing_constant: float = self.alpha * xp.log(self.beta) - float(
+            sp.special.gammaln(self.alpha)
+        )
+
+    def rescale_hyperparameters_to_internal(self, theta, direction):
+        """
+        Transform between external and internal parameter representations.
+
+        The Gamma distribution is defined for positive values, but optimization
+        often works better in unconstrained space. This method transforms
+        between theta (positive) and log(theta) (unconstrained).
+
+        Parameters
+        ----------
+        theta : float or NDArray
+            Parameter value(s) to transform.
+        direction : str
+            Transformation direction:
+            - "forward": theta -> log(theta) (external to internal)
+            - "backward": log(theta) -> theta (internal to external)
+
+        Returns
+        -------
+        float or NDArray
+            Transformed parameter value(s).
+
+        Raises
+        ------
+        ValueError
+            If direction is not "forward" or "backward".
+        """
+        if direction == "forward":
+            theta_scaled = xp.log(theta)
+        elif direction == "backward":
+            theta_scaled = xp.exp(theta)
+        elif direction == "forward_jacobian":
+            theta_scaled = 1 / theta  # d(log(theta))/d(theta) = 1/theta
+        elif direction == "backward_jacobian":
+            theta_scaled = theta  # d(exp(theta))/d(theta) = exp(theta) = theta 
+        else:
+            raise ValueError(f"Unknown direction: {direction}")
+
+        return theta_scaled
+
+    def evaluate_log_prior(self, theta: float, **kwargs) -> float:
+        """
+        Evaluate the log prior probability density.
+
+        Computes the log probability density of the Gamma distribution
+        at the given theta value in its external/user representation.
+
+        Parameters
+        ----------
+        theta : float
+            Parameter value at which to evaluate the log prior.
+            Must be positive (external/user representation).
+        **kwargs
+            Additional keyword arguments (unused).
+
+        Returns
+        -------
+        float
+            Log prior probability density at theta.
+
+        Notes
+        -----
+        The computation follows:
+            log p(θ) = C + (α - 1) * log(θ) - β * θ
+        where C is the normalizing constant, α is the shape parameter,
+        and β is the rate parameter.
+
+        Raises
+        ------
+        ValueError
+            If theta is not positive (implicitly through log computation).
+        """
+
+        log_prior = (
+            self.normalizing_constant
+            + (self.alpha - 1) * xp.log(theta)
+            - self.beta * theta
+        )
+
+        return log_prior
+
+if __name__ == "__main__":
+    """
+    Test Gaussian quadrature for functions using rescale_hyperparameters_to_internal() 
+    from gamma prior hyperparameters.
+    
+    We start with normally distributed random variables in internal space (unconstrained)
+    that get reparametrized to external space (positive) using the gamma prior's 
+    rescaling function.
+    """
+    
+    from dalia.utils.gaussian_quadrature import compute_variance_gauss_hermite
+    
+    print("=" * 80)
+    print("Testing Gaussian Quadrature with Gamma Prior Rescaling")
+    print("=" * 80)
+    
+    # Create a gamma prior configuration
+    config = GammaPriorHyperparametersConfig(alpha=2.0, beta=1.0)
+    gamma_prior = GammaPriorHyperparameters(config=config)
+    
+    # Define parameters for the normal distribution in internal space
+    # These represent log(theta) where theta > 0 is the gamma-distributed parameter
+    mean_internal = 0.5  # Mean of log(theta)
+    variance_internal = 0.25  # Variance of log(theta)
+    
+    print(f"Internal space (log-scale) parameters:")
+    print(f"  Mean: {mean_internal}")
+    print(f"  Variance: {variance_internal}")
+    print(f"  Standard deviation: {np.sqrt(variance_internal)}")
+    print()
+    
+    # Test 1: Compute statistics using Gaussian quadrature
+    print("1. Computing statistics using Gaussian quadrature:")
+    
+    # Use the rescaling function as the transform
+    def transform_func(x, direction):
+        return gamma_prior.rescale_hyperparameters_to_internal(x, direction)
+    
+    # Compute statistics using different numbers of quadrature points
+    for n_points in [10, 20, 30, 50]:
+        result = compute_variance_gauss_hermite(
+            mean_internal, 
+            variance_internal, 
+            transform_func, 
+            n_points=n_points
+        )
+        
+        print(f"  n_points = {n_points:2d}: Mean = {result['mean']:.6f}, "
+              f"Std = {result['std']:.6f}, Var = {result['variance']:.6f}")
+    
+    print()
+    
+    # Test 2: Compare with analytical solution
+    print("2. Comparison with analytical log-normal distribution:")
+    print("   For log(Y) ~ N(μ, σ²), we have:")
+    print("   E[Y] = exp(μ + σ²/2)")
+    print("   Var[Y] = (exp(σ²) - 1) * exp(2μ + σ²)")
+    
+    # Analytical moments for log-normal distribution
+    mu = mean_internal
+    sigma2 = variance_internal
+    
+    analytical_mean = np.exp(mu + sigma2/2)
+    analytical_variance = (np.exp(sigma2) - 1) * np.exp(2*mu + sigma2)
+    analytical_std = np.sqrt(analytical_variance)
+    
+    print(f"   Analytical mean: {analytical_mean:.6f}")
+    print(f"   Analytical std:  {analytical_std:.6f}")
+    print(f"   Analytical var:  {analytical_variance:.6f}")
+    print()
+    
+    # Compare with quadrature result (using 50 points)
+    quad_result = compute_variance_gauss_hermite(
+        mean_internal, variance_internal, transform_func, n_points=50
+    )
+    
+    print("3. Comparison of quadrature vs analytical:")
+    print(f"   Mean difference: {abs(quad_result['mean'] - analytical_mean):.2e}")
+    print(f"   Std difference:  {abs(quad_result['std'] - analytical_std):.2e}")
+    print(f"   Var difference:  {abs(quad_result['variance'] - analytical_variance):.2e}")
+    
+    # Relative errors
+    mean_rel_error = abs(quad_result['mean'] - analytical_mean) / analytical_mean
+    std_rel_error = abs(quad_result['std'] - analytical_std) / analytical_std
+    var_rel_error = abs(quad_result['variance'] - analytical_variance) / analytical_variance
+    
+    print(f"   Mean rel. error: {mean_rel_error:.2e}")
+    print(f"   Std rel. error:  {std_rel_error:.2e}")
+    print(f"   Var rel. error:  {var_rel_error:.2e}")
+    print()
+    
+    # Test 3: Test with different internal parameters
+    print("4. Testing with different internal parameters:")
+    
+    test_cases = [
+        {"mean": 0.0, "var": 0.1, "name": "Small variance"},
+        {"mean": 1.0, "var": 0.5, "name": "Medium variance"},
+        {"mean": -0.5, "var": 1.0, "name": "Large variance"},
+        {"mean": 2.0, "var": 0.01, "name": "Large mean, small variance"}
+    ]
+    
+    for case in test_cases:
+        mu_test = case["mean"]
+        var_test = case["var"]
+        
+        # Quadrature result
+        quad_result = compute_variance_gauss_hermite(
+            mu_test, var_test, transform_func, n_points=30
+        )
+        
+        # Analytical result
+        anal_mean = np.exp(mu_test + var_test/2)
+        anal_var = (np.exp(var_test) - 1) * np.exp(2*mu_test + var_test)
+        
+        rel_mean_error = abs(quad_result['mean'] - anal_mean) / anal_mean
+        rel_var_error = abs(quad_result['variance'] - anal_var) / anal_var
+        
+        print(f"   {case['name']:25s}: Mean rel. err = {rel_mean_error:.2e}, "
+              f"Var rel. err = {rel_var_error:.2e}")
+    
+    print()
+    
+    # Test 4: Test the rescaling function directions
+    print("5. Testing rescaling function directions:")
+    
+    # Test some values
+    test_values = [0.1, 0.5, 1.0, 2.0, 5.0]
+    
+    print("   Testing forward (external -> internal) and backward (internal -> external):")
+    for theta in test_values:
+        # Forward: theta -> log(theta)
+        log_theta = gamma_prior.rescale_hyperparameters_to_internal(theta, "forward")
+        
+        # Backward: log(theta) -> theta  
+        theta_recovered = gamma_prior.rescale_hyperparameters_to_internal(log_theta, "backward")
+        
+        error = abs(theta - theta_recovered)
+        print(f"   θ = {theta:4.1f} -> log(θ) = {log_theta:6.3f} -> θ = {theta_recovered:6.3f}, "
+              f"error = {error:.2e}")
+    
+    print()
+    
+    # Test 5: Convergence study
+    print("6. Convergence study (increasing number of quadrature points):")
+    
+    mu_conv = 0.3
+    var_conv = 0.4
+    analytical_mean_conv = np.exp(mu_conv + var_conv/2)
+    
+    n_points_list = [5, 10, 15, 20, 25, 30, 40, 50, 75, 100]
+    
+    print("   n_points   |   Mean      |  Rel. Error")
+    print("   -----------|-------------|------------")
+    
+    for n in n_points_list:
+        result = compute_variance_gauss_hermite(
+            mu_conv, var_conv, transform_func, n_points=n
+        )
+        rel_error = abs(result['mean'] - analytical_mean_conv) / analytical_mean_conv
+        
+        print(f"   {n:8d}   | {result['mean']:10.6f} | {rel_error:.3e}")
+    
+    print()
+    print("=" * 80)
+    print("Test completed successfully!")
+    print("All tests show that Gaussian quadrature accurately approximates")
+    print("the moments of log-normal distributions obtained through gamma")  
+    print("prior rescaling transformations.")
+    print("=" * 80)
+
diff --git a/src/dalia/prior_hyperparameters/gaussian.py b/src/dalia/prior_hyperparameters/gaussian.py
index d18215ed..f03ef1ab 100644
--- a/src/dalia/prior_hyperparameters/gaussian.py
+++ b/src/dalia/prior_hyperparameters/gaussian.py
@@ -1,6 +1,7 @@
 # Copyright 2024-2025 DALIA authors. All rights reserved.
 from dalia import NDArray
 from scipy.sparse import spmatrix
+from dalia import sp, xp
 
 import numpy as np
 
@@ -11,19 +12,100 @@
 
 
 class GaussianPriorHyperparameters(PriorHyperparameters):
-    """Gaussian prior hyperparameters."""
+    """
+    Univariate Gaussian prior hyperparameters.
+
+    This class implements prior hyperparameters following a univariate normal
+    (Gaussian) distribution with specified mean and precision (inverse variance).
+
+    Parameters
+    ----------
+    config : GaussianPriorHyperparametersConfig
+        Configuration object containing mean and precision parameters.
+
+    Attributes
+    ----------
+    mean : float
+        Mean of the Gaussian distribution.
+    precision : float
+        Precision (inverse variance) of the Gaussian distribution.
+    normalizing_constant : float
+        Precomputed normalizing constant for log probability evaluation.
+    """
 
     def __init__(
         self,
         config: GaussianPriorHyperparametersConfig,
     ) -> None:
-        """Initializes the Gaussian prior hyperparameters."""
+        """
+        Initialize the Gaussian prior hyperparameters.
+
+        Parameters
+        ----------
+        config : GaussianPriorHyperparametersConfig
+            Configuration containing mean and precision parameters.
+
+        Raises
+        ------
+        ValueError
+            If the precision is not positive.
+        """
         super().__init__(config)
 
         self.mean: float = config.mean
         self.precision: float = config.precision
 
+        # Validate precision is positive
+        if self.precision <= 0:
+            raise ValueError(f"Precision must be positive, got {self.precision}")
+
+        self.normalizing_constant = -0.5 * xp.log(2 * xp.pi) + 0.5 * xp.log(
+            self.precision
+        )
+
+    def rescale_hyperparameters_to_internal(self, theta, direction):
+        """
+        Rescale hyperparameters between internal and external representations.
+
+        Parameters
+        ----------
+        theta : NDArray or float
+            Hyperparameter values to rescale.
+        direction : str
+            Direction of rescaling ('forward' or 'backward', i.e. from interpretable/user/external to internal or vice versa).
+
+        Returns
+        -------
+        NDArray or float
+            Rescaled hyperparameter values. In this case it is the identity, therefore unchanged.
+        """
+        return super().rescale_hyperparameters_to_internal(theta, direction)
+
     def evaluate_log_prior(self, theta: float, **kwargs) -> float:
-        """Evaluate the log prior hyperparameters."""
+        """
+        Evaluate the log prior probability density.
+
+        Computes the log probability density of the univariate normal
+        distribution at the given theta value.
+
+        Parameters
+        ----------
+        theta : float
+            Parameter value at which to evaluate the log prior.
+        **kwargs
+            Additional keyword arguments (unused).
+
+        Returns
+        -------
+        float
+            Log prior probability density at theta.
 
-        return -0.5 * self.precision * (theta - self.mean) ** 2
+        Notes
+        -----
+        The computation follows:
+            log p(θ) = C - 0.5 * τ * (θ - μ)²
+        where C is the normalizing constant, τ is precision, and μ is the mean.
+        """
+        return (
+            self.normalizing_constant - 0.5 * self.precision * (theta - self.mean) ** 2
+        )
diff --git a/src/dalia/prior_hyperparameters/gaussian_mvn.py b/src/dalia/prior_hyperparameters/gaussian_mvn.py
index 8310ef10..6bdc9170 100644
--- a/src/dalia/prior_hyperparameters/gaussian_mvn.py
+++ b/src/dalia/prior_hyperparameters/gaussian_mvn.py
@@ -12,13 +12,44 @@
 
 
 class GaussianMVNPriorHyperparameters(PriorHyperparameters):
-    """Gaussian MVN prior hyperparameters."""
+    """
+    Gaussian multivariate normal (MVN) prior hyperparameters.
+
+    This class implements prior hyperparameters following a multivariate normal
+    distribution with specified mean and precision matrix.
+
+    Parameters
+    ----------
+    config : GaussianMVNPriorHyperparametersConfig
+        Configuration object containing mean and precision matrix.
+
+    Attributes
+    ----------
+    mean : NDArray
+        Mean vector of the multivariate normal distribution.
+    precision : spmatrix
+        Precision matrix (inverse covariance) of the distribution.
+    normalizing_constant : float
+        Precomputed normalizing constant for log probability evaluation.
+    """
 
     def __init__(
         self,
         config: GaussianMVNPriorHyperparametersConfig,
     ) -> None:
-        """Initializes the Gaussian MVN prior hyperparameters."""
+        """
+        Initialize the Gaussian MVN prior hyperparameters.
+
+        Parameters
+        ----------
+        config : GaussianMVNPriorHyperparametersConfig
+            Configuration containing mean vector and precision matrix.
+
+        Raises
+        ------
+        ValueError
+            If the precision matrix is not positive definite.
+        """
         super().__init__(config)
 
         self.mean: NDArray = config.mean
@@ -31,8 +62,58 @@ def __init__(
             self.mean: NDArray = xp.asarray(self.mean)
             self.precision: sp.sparse.spmatrix = sp.sparse.csc_matrix(self.precision)
 
-    def evaluate_log_prior(self, theta: float, **kwargs) -> float:
-        """Evaluate the log prior hyperparameters."""
+        sign, logabsdet = np.linalg.slogdet(self.precision.toarray())
+        if sign != 1:
+            raise ValueError("Precision matrix must be positive definite.")
+
+        self.normalizing_constant = (
+            -0.5 * self.mean.shape[0] * xp.log(2 * xp.pi) + 0.5 * logabsdet
+        )
+        print("Normalizing constant: ", self.normalizing_constant)
+
+    def rescale_hyperparameters_to_internal(self, theta, direction):
+        """
+        Rescale hyperparameters between internal and external/user representations.
+
+        Parameters
+        ----------
+        theta : NDArray
+            Hyperparameter values to rescale.
+        direction : str
+            Direction of rescaling ('forward' or 'backward').
+
+        Returns
+        -------
+        NDArray
+            Rescaled hyperparameter values, which is the identity in this case, therefore unchanged.
+        """
+        return super().rescale_hyperparameters_to_internal(theta, direction)
+
+    def evaluate_log_prior(self, theta: NDArray, **kwargs) -> float:
+        """
+        Evaluate the log prior probability density.
+
+        Computes the log probability density of the multivariate normal
+        distribution at the given theta values.
+
+        Parameters
+        ----------
+        theta : NDArray
+            Parameter values at which to evaluate the log prior.
+            Must have the same shape as the mean vector.
+        **kwargs
+            Additional keyword arguments (unused).
+
+        Returns
+        -------
+        float
+            Log prior probability density at theta.
+
+        Raises
+        ------
+        ValueError
+            If theta and mean have incompatible shapes.
+        """
 
         # TODO: add check in config or somewhere else that dim(theta) and dim(mean) match
         if self.mean.shape != theta.shape:
@@ -41,7 +122,11 @@ def evaluate_log_prior(self, theta: float, **kwargs) -> float:
             )
 
         if isinstance(self.mean, float):
-            return -0.5 * self.precision * (theta - self.mean) ** 2
+            return (
+                self.normalizing_constant
+                - 0.5 * self.precision * (theta - self.mean) ** 2
+            )
         else:
-            # neglect constant as the precision is fixed
-            return -0.5 * (theta - self.mean).T @ self.precision @ (theta - self.mean)
+            return self.normalizing_constant - 0.5 * (
+                theta - self.mean
+            ).T @ self.precision @ (theta - self.mean)
diff --git a/src/dalia/prior_hyperparameters/penalized_complexity.py b/src/dalia/prior_hyperparameters/penalized_complexity.py
index 2c9788d2..5bb7baab 100644
--- a/src/dalia/prior_hyperparameters/penalized_complexity.py
+++ b/src/dalia/prior_hyperparameters/penalized_complexity.py
@@ -44,6 +44,9 @@ def __init__(
 
         # print("lambda_theta: ", self.lambda_theta)
 
+    def rescale_hyperparameters_to_internal(self, theta, direction):
+        return super().rescale_hyperparameters_to_internal(theta, direction)
+
     def evaluate_log_prior(self, theta: float, **kwargs) -> float:
         """Evaluate the prior hyperparameters."""
         log_prior: float = 0.0
diff --git a/src/dalia/solvers/sparse_solver.py b/src/dalia/solvers/sparse_solver.py
index 7403cb60..e1b64822 100644
--- a/src/dalia/solvers/sparse_solver.py
+++ b/src/dalia/solvers/sparse_solver.py
@@ -26,6 +26,7 @@ def cholesky(self, A: sp.sparse.spmatrix, **kwargs) -> None:
         if (LU.U.diagonal() > 0).all():  # Check the matrix A is positive definite.
             self.L = LU.L.dot(sp.sparse.diags(LU.U.diagonal() ** 0.5))
         else:
+            print("min(diag(L)): ", xp.min(LU.U.diagonal()))
             raise ValueError("The matrix is not positive definite")
 
     def solve(
@@ -56,13 +57,27 @@ def logdet(
 
         return 2 * xp.sum(xp.log(self.L.diagonal()))
 
-    def selected_inversion(self, **kwargs):
-        # Placeholder for the selected inversion method.
-        return super().selected_inversion(**kwargs)
+    def selected_inversion(self, **kwargs) -> None:
+        # convert to dense
+        L_dense = self.L.toarray()
+        L_inv = xp.eye(self.L.shape[0])
+
+        L_inv[:] = sp.linalg.solve_triangular(
+            L_dense, L_inv, lower=True, overwrite_b=True
+        )
+        self.A_inv = L_inv.T @ L_inv
+
+        return self.A_inv
+
+    def _structured_to_spmatrix(self, A: sp.sparse.spmatrix, **kwargs) -> None:
+        B = A.tocoo()
+        B.data = self.A_inv[B.row, B.col]
+
+        return B
 
     def get_solver_memory(self) -> int:
         """Return the memory used by the solver in number of bytes"""
         if self.L is None:
             return 0
 
-        return self.L.data.nbytes + self.L.indptr.nbytes + self.L.indices.nbytes
\ No newline at end of file
+        return self.L.data.nbytes + self.L.indptr.nbytes + self.L.indices.nbytes
diff --git a/src/dalia/submodels/__init__.py b/src/dalia/submodels/__init__.py
index 2e2b78fb..c19f8b81 100644
--- a/src/dalia/submodels/__init__.py
+++ b/src/dalia/submodels/__init__.py
@@ -4,5 +4,12 @@
 from dalia.submodels.spatial import SpatialSubModel
 from dalia.submodels.spatio_temporal import SpatioTemporalSubModel
 from dalia.submodels.brainiac import BrainiacSubModel
+from dalia.submodels.ar1 import AR1SubModel
 
-__all__ = ["RegressionSubModel", "SpatialSubModel", "SpatioTemporalSubModel", "BrainiacSubModel"]
+__all__ = [
+    "RegressionSubModel",
+    "SpatialSubModel",
+    "SpatioTemporalSubModel",
+    "BrainiacSubModel",
+    "AR1SubModel",
+]
diff --git a/src/dalia/submodels/ar1.py b/src/dalia/submodels/ar1.py
new file mode 100644
index 00000000..198f28b4
--- /dev/null
+++ b/src/dalia/submodels/ar1.py
@@ -0,0 +1,77 @@
+# Copyright 2024-2025 DALIA authors. All rights reserved.
+from tabulate import tabulate
+
+import numpy as np
+
+from dalia import sp, xp
+from dalia.configs.submodels_config import AR1SubModelConfig
+from dalia.core.submodel import SubModel
+from dalia.utils import add_str_header, scaled_logit
+
+
+class AR1SubModel(SubModel):
+    """Fit an AR(1) model."""
+
+    def __init__(
+        self,
+        config: AR1SubModelConfig,
+    ) -> None:
+        """Initializes the model."""
+        super().__init__(config)
+
+        # check that dimensions match
+
+
+    def construct_Q_prior(self, **kwargs) -> sp.sparse.coo_matrix:
+        """Construct the prior precision matrix."""
+
+        ## TODO: link this to prior hyperparameters class
+        ## i.e. rescale phi and tau according to prior hyperparameters class
+        phi = kwargs.get("phi")
+        #phi = scaled_logit(phi_scaled, direction="backward")
+        #print("phi: ", phi)
+
+        tau = kwargs.get("tau")
+        #exp_tau = xp.exp(tau)
+        #print("tau: ", tau)
+
+        s2 = 1 / tau
+        denom = s2 * (1 - phi**2)
+
+        diag = [(1 + phi**2) / denom] * self.n_latent_parameters
+        diag[0] = diag[-1] = 1 / denom
+        off_diag = [-phi / denom] * (self.n_latent_parameters - 1)
+
+        Q_prior = sp.sparse.diags([off_diag, diag, off_diag], [-1, 0, 1])
+        
+        # need this -> otherwise there might be a sorting issue
+        Q_prior = Q_prior.tocsr()
+        Q_prior.sort_indices()
+
+        return Q_prior.tocoo()
+
+    def __str__(self) -> str:
+        """String representation of the submodel."""
+        str_representation = ""
+
+        # --- Make the Submodel table ---
+        values = [
+            ["Submodel Type", self.submodel_type],
+            ["Number of Latent Parameters", self.n_latent_parameters],
+            ["Phi", f"{self.config.phi:.3f}"],
+            ["tau", f"{self.config.tau:.3f}"],
+        ]
+        submodel_table = tabulate(
+            values,
+            tablefmt="fancy_grid",
+            colalign=("left", "center"),
+        )
+
+        # Add the header title
+        submodel_table = add_str_header(
+            title=self.submodel_type.replace("_", " ").title(),
+            table=submodel_table,
+        )
+        str_representation += submodel_table
+
+        return str_representation
diff --git a/src/dalia/utils/__init__.py b/src/dalia/utils/__init__.py
index b1a25fe8..34063a93 100644
--- a/src/dalia/utils/__init__.py
+++ b/src/dalia/utils/__init__.py
@@ -12,6 +12,9 @@
 )
 from dalia.utils.host import get_host_configuration
 from dalia.utils.link_functions import cloglog, scaled_logit, sigmoid
+from dalia.utils.correlation import compute_outer_covariance_matrix
+from dalia.utils.gaussian_quadrature import compute_variance_gauss_hermite
+from dalia.utils.bivariate_gaussian_quadrature import compute_bivariate_expectation
 from dalia.utils.multiprocessing import (
     allreduce,
     allgather,
@@ -21,6 +24,7 @@
     smartsplit,
     synchronize,
     synchronize_gpu,
+    check_vector_consistency,
     DummyCommunicator,
 )
 from dalia.utils.spmatrix_utils import bdiag_tiling, extract_diagonal, memory_footprint
@@ -36,6 +40,9 @@
     "sigmoid",
     "cloglog",
     "scaled_logit",
+    "compute_outer_covariance_matrix",
+    "compute_variance_gauss_hermite",
+    "compute_bivariate_expectation",
     "print_msg",
     "synchronize",
     "synchronize_gpu",
@@ -44,6 +51,7 @@
     "allreduce",
     "allgather",
     "bcast",
+    "check_vector_consistency",
     "bdiag_tiling",
     "extract_diagonal",
     "memory_footprint",
diff --git a/src/dalia/utils/bivariate_gaussian_quadrature.py b/src/dalia/utils/bivariate_gaussian_quadrature.py
new file mode 100644
index 00000000..4150737d
--- /dev/null
+++ b/src/dalia/utils/bivariate_gaussian_quadrature.py
@@ -0,0 +1,236 @@
+# Copyright 2024-2025 DALIA authors. All rights reserved.
+
+import numpy as np
+from dalia import xp
+from scipy.special import roots_hermite
+
+from dalia.utils.gaussian_quadrature import compute_variance_gauss_hermite
+
+def compute_bivariate_expectation(func1, func2, rho, n_points=20):
+    """
+    Compute E[f(Z₁, Z₂)] where (Z₁, Z₂) ~ N(0, Σ) using bivariate Gauss-Hermite quadrature.
+    
+    E[f(Z₁, Z₂)] = ∬ f(z₁, z₂) φ_ρ(z₁, z₂) dz₁ dz₂
+    
+    where φ_ρ(z₁, z₂) is the bivariate standard normal density with correlation ρ.
+    
+    Parameters
+    ----------
+    func : callable
+        Function f(z₁, z₂) to compute expectation of
+    rho : float, optional
+        Correlation coefficient between Z₁ and Z₂ (default: 0.0)
+    n_points : int, optional
+        Number of quadrature points per dimension (default: 20)
+        
+    Returns
+    -------
+    float
+        Expected value E[f(Z₁, Z₂)]
+    """
+    # Get Gauss-Hermite quadrature points and weights
+    nodes, weights = roots_hermite(n_points)
+    
+    # Transform nodes from Hermite polynomial roots to standard normal
+    z_nodes = xp.sqrt(2) * nodes
+    adjusted_weights = weights / xp.sqrt(xp.pi)
+    
+    # For bivariate case with correlation ρ, we need to transform to correlated variables
+    # If (U₁, U₂) are independent N(0,1), then:
+    # Z₁ = U₁
+    # Z₂ = ρU₁ + √(1-ρ²)U₂
+    # gives (Z₁, Z₂) ~ N(0, [[1, ρ], [ρ, 1]])
+    
+    if abs(rho) > 1:
+        raise ValueError("Correlation coefficient rho must be in [-1, 1]")
+    
+    sqrt_one_minus_rho_sq = xp.sqrt(1 - rho**2) if abs(rho) < 1 else 0.0
+    
+    expectation = 0.0
+    
+    # Double loop over all combinations of quadrature points
+    for i, (u1, w1) in enumerate(zip(z_nodes, adjusted_weights)):
+        for j, (u2, w2) in enumerate(zip(z_nodes, adjusted_weights)):
+            # Transform to correlated variables
+            z1 = u1
+            z2 = rho * u1 + sqrt_one_minus_rho_sq * u2
+            
+            # Evaluate function at transformed points
+            f_val = func1(z1) * func2(z2)
+            
+            # Combined weight (product of univariate weights)
+            weight = w1 * w2
+            
+            expectation += weight * f_val
+    
+    return expectation
+
+
+def test_quadrature_accuracy():
+    """Test quadrature accuracy against known analytical results."""
+    
+    print("=" * 80)
+    print("Testing Bivariate Gaussian Quadrature")
+    print("=" * 80)
+        
+    # Test 2: Bivariate expectations with independence (ρ = 0)
+    print("2. Testing bivariate expectations with independence (ρ = 0):")
+    
+    # E[1] = 1
+    biv_expectation_1 = compute_bivariate_expectation(lambda z: 1.0, lambda z: 1.0, rho=0.0, n_points=10)
+    print(f"   E[1] = {biv_expectation_1:.8f} (should be 1.0, error: {abs(1.0 - biv_expectation_1):.2e})")
+    
+    # E[Z₁Z₂] = 0 when independent
+    biv_expectation_z1z2 = compute_bivariate_expectation(lambda z: z, lambda z: z, rho=0.0, n_points=20)
+    print(f"   E[Z₁Z₂] = {biv_expectation_z1z2:.8f} (should be 0.0, error: {abs(biv_expectation_z1z2):.2e})")
+    
+    # E[Z₁² + Z₂²] = E[Z₁²] + E[Z₂²] = 1 + 1 = 2 (using independence)
+    biv_expectation_sq = compute_bivariate_expectation(lambda z: z**2, lambda z: z**2, rho=0.0, n_points=20)
+    print(f"   E[Z₁² + Z₂²] = {biv_expectation_sq:.8f} (should be 2.0, error: {abs(2.0 - biv_expectation_sq):.2e})")
+    print()
+    
+    # Test 3: Product expectations with independence
+    print("3. Testing product expectations E[f₁(Z₁)f₂(Z₂)] with independence:")
+    
+    # E[Z₁ * Z₂] = E[Z₁] * E[Z₂] = 0 * 0 = 0 when independent
+    prod_exp_z1z2 = compute_bivariate_expectation(lambda z: z, lambda z: z, rho=0.0, n_points=20)
+    print(f"   E[Z₁ * Z₂] = {prod_exp_z1z2:.8f} (should be 0.0, error: {abs(prod_exp_z1z2):.2e})")
+    
+    # E[Z₁² * 1] = E[Z₁²] * E[1] = 1 * 1 = 1
+    prod_exp_z1sq_1 = compute_bivariate_expectation(lambda z: z**2, lambda z: 1.0, rho=0.0, n_points=20)
+    print(f"   E[Z₁² * 1] = {prod_exp_z1sq_1:.8f} (should be 1.0, error: {abs(1.0 - prod_exp_z1sq_1):.2e})")
+    
+    # E[exp(Z₁) * exp(Z₂)] = E[exp(Z₁)] * E[exp(Z₂)] = exp(0.5) * exp(0.5) = exp(1.0)
+    prod_exp_exp = compute_bivariate_expectation(lambda z: xp.exp(z), lambda z: xp.exp(z), rho=0.0, n_points=30)
+    analytical_prod_exp = xp.exp(1.0)
+    print(f"   E[exp(Z₁) * exp(Z₂)] = {prod_exp_exp:.8f} (should be {analytical_prod_exp:.8f}, error: {abs(analytical_prod_exp - prod_exp_exp):.2e})")
+    print()
+    
+    # Test 4: Bivariate expectations with correlation
+    print("4. Testing bivariate expectations with correlation:")
+    
+    correlations = [-0.9, -0.5, -0.3, 0.0, 0.3, 0.5, 0.8, 0.9]
+    
+    print("   Correlation | E[Z₁Z₂]    | Analytical | Error")
+    print("   ------------|------------|------------|----------")
+    
+    for rho in correlations:
+        # E[Z₁Z₂] = ρ for bivariate normal
+        biv_exp_corr = compute_bivariate_expectation(lambda z: z, lambda z: z, rho=rho, n_points=25)
+        error = abs(rho - biv_exp_corr)
+        print(f"   {rho:10.1f} | {biv_exp_corr:10.6f} | {rho:10.6f} | {error:.2e}")
+    print()
+    
+    # Test 4b: Extreme correlation values (±1)
+    print("4b. Testing extreme correlation values (ρ = ±1):")
+    print("    Note: Perfect correlation means Z₂ = ±Z₁")
+    
+    extreme_correlations = [-1.0, 1.0]
+    
+    print("   Correlation | E[Z₁Z₂]    | Analytical | Error     | Note")
+    print("   ------------|------------|------------|----------|------------------")
+    
+    for rho in extreme_correlations:
+        # For extreme correlations, use more quadrature points for accuracy
+        biv_exp_corr = compute_bivariate_expectation(lambda z: z, lambda z: z, rho=rho, n_points=40)
+        error = abs(rho - biv_exp_corr)
+        note = "Perfect positive" if rho == 1.0 else "Perfect negative"
+        print(f"   {rho:10.1f} | {biv_exp_corr:10.6f} | {rho:10.6f} | {error:.2e} | {note}")
+        
+        # Additional test: For ρ = ±1, E[Z₁²Z₂²] should equal E[Z₁⁴] = 3
+        z1_sq_z2_sq = compute_bivariate_expectation(lambda z: z**2, lambda z: z**2, rho=rho, n_points=40)
+        expected_z1_4 = 3.0  # Fourth moment of standard normal
+        error_z4 = abs(z1_sq_z2_sq - expected_z1_4)
+        print(f"   {rho:10.1f} | E[Z₁²Z₂²]={z1_sq_z2_sq:6.4f} | E[Z₁⁴]={expected_z1_4:6.1f} | {error_z4:.2e} | Should equal E[Z₁⁴]")
+    
+    print()
+    
+    # Test 5: Product expectations with correlation  
+    print("5. Testing product expectations with correlation:")
+    print("   For E[f₁(Z₁)f₂(Z₂)], correlation affects the result when f₁ and f₂ are nonlinear")
+    
+    def f1(z):
+        return z**2
+    
+    def f2(z):  
+        return z**3
+    
+    print("   E[Z₁² * Z₂³] for different correlations:")
+    print("   Correlation | E[Z₁²Z₂³]")
+    print("   ------------|----------")
+    
+    for rho in [-0.9, -0.5, 0.0, 0.5, 0.9]:
+        prod_exp_nonlinear = compute_bivariate_expectation(f1, f2, rho=rho, n_points=30)
+        print(f"   {rho:10.1f} | {prod_exp_nonlinear:10.6f}")
+    print()
+    
+    # Test 6: Convergence study
+    print("6. Convergence study (increasing number of quadrature points):")
+    
+    # For this test, we approximate E[exp(0.1*Z₁ + 0.2*Z₂)] using E[exp(0.1*Z₁) * exp(0.2*Z₂)]
+    # This is exact since exp(a+b) = exp(a)*exp(b)
+    def func1_test(z):
+        return xp.exp(0.1 * z)
+    
+    def func2_test(z):
+        return xp.exp(0.2 * z)
+    
+    # Analytical result for E[exp(aZ₁ + bZ₂)] with (Z₁,Z₂) ~ N(0, Σ)
+    # For a=0.1, b=0.2, ρ=0.5: E[exp(aZ₁ + bZ₂)] = exp(0.5 * (a² + b² + 2abρ))
+    a, b, rho_test = 0.1, 0.2, 0.5
+    analytical_mgf = xp.exp(0.5 * (a**2 + b**2 + 2*a*b*rho_test))
+    
+    print(f"   Testing E[exp(0.1*Z₁) * exp(0.2*Z₂)] with ρ = {rho_test}")
+    print(f"   Analytical result: {analytical_mgf:.8f}")
+    print()
+    print("   n_points | Numerical    | Error")
+    print("   ---------|--------------|----------")
+    
+    for n in [5, 10, 15, 20, 25, 30]:
+        numerical_mgf = compute_bivariate_expectation(func1_test, func2_test, rho=rho_test, n_points=n)
+        error = abs(analytical_mgf - numerical_mgf)
+        print(f"   {n:7d}  | {numerical_mgf:12.8f} | {error:.2e}")
+    
+    print()
+    
+    # Test 7: Practical example with transformations
+    print("7. Practical example: Log-normal random variables")
+    print("   Let Y₁ = exp(Z₁), Y₂ = exp(Z₂) where (Z₁, Z₂) have correlation ρ")
+    print("   Computing E[Y₁ * Y₂] = E[exp(Z₁ + Z₂)]")
+    
+    rho_examples = [-0.8, -0.3, 0.0, 0.5, 0.9]
+    
+    print("   Correlation | E[Y₁Y₂] Numerical | E[Y₁Y₂] Analytical | Error")
+    print("   ------------|------------------|-------------------|----------")
+    
+    for rho in rho_examples:
+        # Numerical computation
+        numerical_lognormal = compute_bivariate_expectation(
+            lambda z: xp.exp(z), lambda z: xp.exp(z), 
+            rho=rho, n_points=30
+        )
+        
+        # Analytical: E[exp(Z₁ + Z₂)] = exp(E[Z₁ + Z₂] + 0.5*Var[Z₁ + Z₂])
+        # E[Z₁ + Z₂] = 0, Var[Z₁ + Z₂] = Var[Z₁] + Var[Z₂] + 2*Cov[Z₁,Z₂] = 1 + 1 + 2*ρ = 2 + 2*ρ
+        analytical_lognormal = xp.exp(0.5 * (2 + 2*rho))
+        
+        error = abs(numerical_lognormal - analytical_lognormal)
+        
+        print(f"   {rho:10.1f} | {numerical_lognormal:16.8f} | {analytical_lognormal:17.8f} | {error:.2e}")
+    
+    print()
+    print("=" * 80)
+    print("Bivariate Gaussian Quadrature Test Summary:")
+    print("✓ Univariate expectations computed accurately")
+    print("✓ Bivariate expectations with independence verified")
+    print("✓ Product expectations working correctly")
+    print("✓ Correlation effects properly captured") 
+    print("✓ Convergence behavior as expected")
+    print("✓ Practical log-normal example validated")
+    print()
+    print("All bivariate Gaussian quadrature functions are working correctly!")
+    print("=" * 80)
+
+
+if __name__ == "__main__":
+    test_quadrature_accuracy()
\ No newline at end of file
diff --git a/src/dalia/utils/gaussian_quadrature.py b/src/dalia/utils/gaussian_quadrature.py
new file mode 100644
index 00000000..985497cb
--- /dev/null
+++ b/src/dalia/utils/gaussian_quadrature.py
@@ -0,0 +1,382 @@
+# Copyright 2024-2025 DALIA authors. All rights reserved.
+
+from dalia import xp
+
+from scipy.special import roots_hermite
+
+
+# Dummy classes to avoid circular imports during testing
+class DummyConfig:
+    def __init__(self, alpha=2.0, beta=1.0):
+        self.alpha = alpha
+        self.beta = beta
+
+class DummyGammaPriorHyperparameters:
+    def __init__(self, config):
+        self.config = config
+    
+    def rescale_hyperparameters_to_internal(self, theta, direction):
+        """Log transformation: external (positive) <-> internal (unconstrained)"""
+        if direction == "forward":
+            return xp.log(theta)
+        elif direction == "backward":
+            return xp.exp(theta)
+        else:
+            raise ValueError(f"Unknown direction: {direction}")
+
+def compute_variance_gauss_hermite(mean_internal, variance_internal, transform, n_points=20):
+    """
+    Compute variance of transformed distribution using Gauss-Hermite quadrature
+    
+    For a distribution Y where log(Y) ~ N(μ, σ²), we want to compute:
+    Var(Y) = E[Y²] - (E[Y])²
+    
+    Using Gauss-Hermite quadrature by transforming to standard normal form.
+    
+    Theory:
+    If Z ~ N(0,1), then X = μ + σZ ~ N(μ, σ²)
+    So Y = φ⁻¹(X) = φ⁻¹(μ + σZ) = exp(μ + σZ)
+    
+    E[Y] = E[exp(μ + σZ)] = exp(μ) E[exp(σZ)]
+    E[Y²] = E[exp(2μ + 2σZ)] = exp(2μ) E[exp(2σZ)]
+    
+    Where E[exp(aZ)] can be computed using Gauss-Hermite quadrature.
+    """
+    
+    # Get Gauss-Hermite quadrature points and weights
+    nodes, weights = roots_hermite(n_points)
+    
+    # Transform nodes from Hermite polynomial roots to standard normal
+    # Hermite nodes are for exp(-x²), we want exp(-x²/2)/√(2π)
+    # So we scale by √2: z = √2 * node
+    z_nodes = xp.sqrt(2) * nodes
+    
+    # Compute transformed values Y = φ⁻¹(μ + σZ) for each node
+    internal_values = mean_internal + xp.sqrt(variance_internal) * z_nodes
+    y_values = transform(internal_values, direction="backward")
+    
+    # Gauss-Hermite weights need to be adjusted for standard normal
+    # Original: ∫ f(x) exp(-x²) dx ≈ Σ w_i f(x_i)
+    # For standard normal: ∫ f(z) (1/√(2π)) exp(-z²/2) dz
+    # After substitution x = z/√2: ∫ f(√2 x) (1/√π) exp(-x²) dx
+    adjusted_weights = weights / xp.sqrt(xp.pi)
+    
+    # Compute first and second moments
+    mean_y = xp.sum(adjusted_weights * y_values)
+    second_moment_y = xp.sum(adjusted_weights * y_values**2)
+
+    # Variance = E[Y²] - (E[Y])²
+    variance_y = second_moment_y - mean_y**2
+    
+    return {
+        'mean': mean_y,
+        'second_moment': second_moment_y,
+        'variance': variance_y,
+        'std': xp.sqrt(variance_y)
+    }
+
+
+def test_gaussian_quadrature():
+    """
+    Comprehensive test suite for the Gaussian quadrature function.
+    
+    Tests multiple transformation functions and compares results with 
+    analytical solutions where available.
+    """
+    import numpy as np
+    
+    print("=" * 80)
+    print("COMPREHENSIVE GAUSSIAN QUADRATURE TESTS")
+    print("=" * 80)
+    
+    # Test parameters
+    test_tolerance = 1e-6
+    n_quad_points = 50
+    
+    # Test 1: Identity transformation (should recover original normal moments)
+    print("1. Testing Identity Transformation (X = Y)")
+    print("-" * 50)
+    
+    def identity_transform(x, direction):
+        return x  # No transformation
+    
+    mu = 2.0
+    sigma2 = 1.5
+    
+    result = compute_variance_gauss_hermite(mu, sigma2, identity_transform, n_quad_points)
+    
+    # For identity, we should recover the original normal distribution moments
+    expected_mean = mu
+    expected_variance = sigma2
+    
+    print(f"  Expected mean: {expected_mean:.6f}, Got: {result['mean']:.6f}")
+    print(f"  Expected var:  {expected_variance:.6f}, Got: {result['variance']:.6f}")
+    
+    mean_error = abs(result['mean'] - expected_mean)
+    var_error = abs(result['variance'] - expected_variance)
+    
+    print(f"  Mean error: {mean_error:.2e}")
+    print(f"  Var error:  {var_error:.2e}")
+    
+    assert mean_error < test_tolerance, f"Identity mean test failed: error = {mean_error}"
+    assert var_error < test_tolerance, f"Identity variance test failed: error = {var_error}"
+    print("  ✓ Identity transformation test PASSED")
+    print()
+    
+    # Test 2: Log transformation (log-normal distribution)
+    print("2. Testing Log Transformation (Y = exp(X))")
+    print("-" * 50)
+    
+    def log_transform(x, direction):
+        if direction == "forward":
+            return xp.log(x)
+        elif direction == "backward":
+            return xp.exp(x)
+        else:
+            raise ValueError(f"Unknown direction: {direction}")
+    
+    mu = 0.5
+    sigma2 = 0.25
+    
+    result = compute_variance_gauss_hermite(mu, sigma2, log_transform, n_quad_points)
+    
+    # Analytical moments for log-normal: if log(Y) ~ N(μ, σ²)
+    expected_mean = np.exp(mu + sigma2/2)
+    expected_variance = (np.exp(sigma2) - 1) * np.exp(2*mu + sigma2)
+    
+    print(f"  Expected mean: {expected_mean:.6f}, Got: {result['mean']:.6f}")
+    print(f"  Expected var:  {expected_variance:.6f}, Got: {result['variance']:.6f}")
+    
+    mean_error = abs(result['mean'] - expected_mean) / expected_mean
+    var_error = abs(result['variance'] - expected_variance) / expected_variance
+    
+    print(f"  Relative mean error: {mean_error:.2e}")
+    print(f"  Relative var error:  {var_error:.2e}")
+    
+    assert mean_error < 1e-4, f"Log-normal mean test failed: rel error = {mean_error}"
+    assert var_error < 1e-4, f"Log-normal variance test failed: rel error = {var_error}"
+    print("  ✓ Log transformation test PASSED")
+    print()
+    
+    # Test 3: Linear transformation (Y = aX + b)
+    print("3. Testing Linear Transformation (Y = aX + b)")
+    print("-" * 50)
+    
+    a, b = 3.0, -1.5
+    
+    def linear_transform(x, direction):
+        if direction == "forward":
+            return (x - b) / a
+        elif direction == "backward":
+            return a * x + b
+        else:
+            raise ValueError(f"Unknown direction: {direction}")
+    
+    mu = 1.0
+    sigma2 = 0.8
+    
+    result = compute_variance_gauss_hermite(mu, sigma2, linear_transform, n_quad_points)
+    
+    # For Y = aX + b where X ~ N(μ, σ²): E[Y] = aμ + b, Var[Y] = a²σ²
+    expected_mean = a * mu + b
+    expected_variance = a**2 * sigma2
+    
+    print(f"  Transform: Y = {a}X + {b}")
+    print(f"  Expected mean: {expected_mean:.6f}, Got: {result['mean']:.6f}")
+    print(f"  Expected var:  {expected_variance:.6f}, Got: {result['variance']:.6f}")
+    
+    mean_error = abs(result['mean'] - expected_mean)
+    var_error = abs(result['variance'] - expected_variance)
+    
+    print(f"  Mean error: {mean_error:.2e}")
+    print(f"  Var error:  {var_error:.2e}")
+    
+    assert mean_error < test_tolerance, f"Linear mean test failed: error = {mean_error}"
+    assert var_error < test_tolerance, f"Linear variance test failed: error = {var_error}"
+    print("  ✓ Linear transformation test PASSED")
+    print()
+    
+    # Test 4: Quadratic transformation (Y = X²)
+    print("4. Testing Quadratic Transformation (Y = X²)")
+    print("-" * 50)
+    
+    def quadratic_transform(x, direction):
+        if direction == "forward":
+            return xp.sqrt(x)  # Only works for x >= 0
+        elif direction == "backward":
+            return x**2
+        else:
+            raise ValueError(f"Unknown direction: {direction}")
+    
+    mu = 0.0  # Centered to avoid issues with sqrt
+    sigma2 = 0.5
+    
+    result = compute_variance_gauss_hermite(mu, sigma2, quadratic_transform, n_quad_points)
+    
+    # For Y = X² where X ~ N(0, σ²): E[Y] = σ², Var[Y] = 2σ⁴
+    expected_mean = sigma2
+    expected_variance = 2 * sigma2**2
+    
+    print(f"  X ~ N({mu}, {sigma2})")
+    print(f"  Expected mean: {expected_mean:.6f}, Got: {result['mean']:.6f}")
+    print(f"  Expected var:  {expected_variance:.6f}, Got: {result['variance']:.6f}")
+    
+    mean_error = abs(result['mean'] - expected_mean) / expected_mean
+    var_error = abs(result['variance'] - expected_variance) / expected_variance
+    
+    print(f"  Relative mean error: {mean_error:.2e}")
+    print(f"  Relative var error:  {var_error:.2e}")
+    
+    assert mean_error < 1e-3, f"Quadratic mean test failed: rel error = {mean_error}"
+    assert var_error < 1e-2, f"Quadratic variance test failed: rel error = {var_error}"
+    print("  ✓ Quadratic transformation test PASSED")
+    print()
+    
+    # Test 5: Probit transformation (Y = Φ(X), where Φ is the standard normal CDF)
+    print("5. Testing Probit Transformation (Y = Φ(X))")
+    print("-" * 50)
+    
+    from scipy.stats import norm
+    
+    def probit_transform(x, direction):
+        if direction == "forward":
+            # Inverse probit: Φ⁻¹(x) = norm.ppf(x)
+            # Clip to avoid numerical issues at boundaries
+            x_clipped = xp.clip(x, 1e-15, 1-1e-15)
+            return norm.ppf(x_clipped)
+        elif direction == "backward":
+            # Probit: Φ(x) = norm.cdf(x)
+            return norm.cdf(x)
+        else:
+            raise ValueError(f"Unknown direction: {direction}")
+    
+    # For probit transformation, we need to be careful about the domain
+    mu_probit = 0.0  # Center at 0 for symmetry
+    sigma2_probit = 0.5  # Moderate variance to avoid extreme values
+    
+    result = compute_variance_gauss_hermite(mu_probit, sigma2_probit, probit_transform, n_quad_points)
+    
+    # For Y = Φ(X) where X ~ N(0, σ²), we can compute this numerically
+    # Since there's no closed form, we'll use a high-precision reference calculation
+    
+    # Reference calculation using many quadrature points
+    ref_result = compute_variance_gauss_hermite(mu_probit, sigma2_probit, probit_transform, 200)
+    
+    print(f"  X ~ N({mu_probit}, {sigma2_probit})")
+    print(f"  Reference mean (200 pts): {ref_result['mean']:.6f}")
+    print(f"  Test mean ({n_quad_points} pts):     {result['mean']:.6f}")
+    print(f"  Reference var (200 pts):  {ref_result['variance']:.6f}")
+    print(f"  Test var ({n_quad_points} pts):      {result['variance']:.6f}")
+    
+    mean_error = abs(result['mean'] - ref_result['mean']) / ref_result['mean']
+    var_error = abs(result['variance'] - ref_result['variance']) / ref_result['variance']
+    
+    print(f"  Relative mean error: {mean_error:.2e}")
+    print(f"  Relative var error:  {var_error:.2e}")
+    
+    # For probit, we expect the mean to be close to 0.5 due to symmetry
+    expected_mean_approx = 0.5
+    mean_deviation = abs(result['mean'] - expected_mean_approx)
+    
+    print(f"  Expected mean ≈ 0.5, deviation: {mean_deviation:.4f}")
+    
+    assert mean_error < 1e-3, f"Probit mean test failed: rel error = {mean_error}"
+    assert var_error < 1e-2, f"Probit variance test failed: rel error = {var_error}"
+    assert mean_deviation < 0.1, f"Probit mean should be close to 0.5: deviation = {mean_deviation}"
+    print("  ✓ Probit transformation test PASSED")
+    print()
+    
+    # Test 6: Gamma prior rescaling
+    print("6. Testing Gamma Prior Rescaling")
+    print("-" * 50)
+    
+    config = DummyConfig(alpha=2.0, beta=1.0)
+    gamma_prior = DummyGammaPriorHyperparameters(config=config)
+    
+    def gamma_rescale(x, direction):
+        return gamma_prior.rescale_hyperparameters_to_internal(x, direction)
+    
+    mu = 0.2
+    sigma2 = 0.3
+    
+    result = compute_variance_gauss_hermite(mu, sigma2, gamma_rescale, n_quad_points)
+    
+    # This is the same as log-normal since rescale uses exp transformation
+    expected_mean = np.exp(mu + sigma2/2)
+    expected_variance = (np.exp(sigma2) - 1) * np.exp(2*mu + sigma2)
+    
+    print(f"  Expected mean: {expected_mean:.6f}, Got: {result['mean']:.6f}")
+    print(f"  Expected var:  {expected_variance:.6f}, Got: {result['variance']:.6f}")
+    
+    mean_error = abs(result['mean'] - expected_mean) / expected_mean
+    var_error = abs(result['variance'] - expected_variance) / expected_variance
+    
+    print(f"  Relative mean error: {mean_error:.2e}")
+    print(f"  Relative var error:  {var_error:.2e}")
+    
+    assert mean_error < 1e-4, f"Gamma rescale mean test failed: rel error = {mean_error}"
+    assert var_error < 1e-4, f"Gamma rescale variance test failed: rel error = {var_error}"
+    print("  ✓ Gamma prior rescaling test PASSED")
+    print()
+    
+    # Test 7: Convergence with increasing quadrature points
+    print("7. Testing Convergence with Quadrature Points")
+    print("-" * 50)
+    
+    mu_conv = 0.1
+    sigma2_conv = 0.4
+    expected_mean_conv = np.exp(mu_conv + sigma2_conv/2)
+    
+    n_points_list = [5, 10, 15, 20, 30, 50, 75]
+    errors = []
+    
+    print("  n_points  |   Mean      |  Rel. Error")
+    print("  ----------|-------------|------------")
+    
+    for n in n_points_list:
+        result = compute_variance_gauss_hermite(mu_conv, sigma2_conv, log_transform, n)
+        rel_error = abs(result['mean'] - expected_mean_conv) / expected_mean_conv
+        errors.append(rel_error)
+        
+        print(f"  {n:8d}  | {result['mean']:10.6f} | {rel_error:.3e}")
+    
+    # Check that errors generally decrease (allowing some numerical noise)
+    improving = sum(errors[i+1] < errors[i] * 1.1 for i in range(len(errors)-1))
+    improvement_rate = improving / (len(errors) - 1)
+    
+    print(f"  Improvement rate: {improvement_rate:.1%}")
+    assert improvement_rate > 0.6, f"Convergence test failed: improvement rate = {improvement_rate}"
+    print("  ✓ Convergence test PASSED")
+    print()
+    
+    # Test 8: Edge cases
+    print("8. Testing Edge Cases")
+    print("-" * 50)
+    
+    # Small variance
+    result_small = compute_variance_gauss_hermite(1.0, 1e-6, log_transform, 20)
+    expected_small = np.exp(1.0)  # When σ² → 0, E[exp(X)] → exp(μ)
+    error_small = abs(result_small['mean'] - expected_small) / expected_small
+    
+    print(f"  Small variance test: rel error = {error_small:.3e}")
+    assert error_small < 1e-3, f"Small variance test failed: rel error = {error_small}"
+    
+    # Large variance (but not too large to avoid overflow)
+    result_large = compute_variance_gauss_hermite(0.0, 2.0, log_transform, 100)
+    expected_large = np.exp(1.0)  # E[exp(X)] = exp(μ + σ²/2) = exp(0 + 2/2) = e
+    error_large = abs(result_large['mean'] - expected_large) / expected_large
+    
+    print(f"  Large variance test: rel error = {error_large:.3e}")
+    assert error_large < 1e-2, f"Large variance test failed: rel error = {error_large}"
+    
+    print("  ✓ Edge cases test PASSED")
+    print()
+    
+    # Summary
+    print("=" * 80)
+    print("ALL TESTS PASSED! ✓")
+    print("=" * 80)
+
+
+if __name__ == "__main__":
+    test_gaussian_quadrature()
\ No newline at end of file
diff --git a/src/dalia/utils/multiprocessing.py b/src/dalia/utils/multiprocessing.py
index 9ceeb6d6..43a27206 100644
--- a/src/dalia/utils/multiprocessing.py
+++ b/src/dalia/utils/multiprocessing.py
@@ -134,6 +134,7 @@ def bcast(
     comm (CommunicatorType), optional:
         The communication group. Default is MPI.COMM_WORLD.
     """
+
     if backend_flags["mpi_avail"]:
         comm.Bcast(data, root=root)
 
@@ -209,3 +210,44 @@ def smartsplit(
         color_new_group = 0
 
     return active_comm, comm_new_group, color_new_group
+
+def check_vector_consistency(
+    theta: ArrayLike,
+    comm,
+):
+    """
+    Check if all processes have the same theta.
+
+    Parameters:
+    -----------
+    theta (ArrayLike):
+        The theta to check.
+    comm (CommunicatorType), optional:
+        The communication group. Default is MPI.COMM_WORLD.
+    """
+
+    synchronize(comm = comm)
+
+    theta_ref = theta.copy()
+    bcast(theta_ref, root=0, comm=comm)
+
+    if backend_flags["cupy_avail"]:
+        norm_diff = cp.linalg.norm(theta - theta_ref)
+    else:
+        norm_diff = np.linalg.norm(theta - theta_ref)
+
+    if norm_diff > 1e-10:
+        # Print indices and values where theta and theta_ref differ
+        if backend_flags["cupy_avail"]:
+            diff_indices = cp.where(theta != theta_ref)[0]
+            for idx in diff_indices:
+                print(f"Process {comm.Get_rank()} difference at index {idx}: theta_ref={theta_ref[idx]}, theta={theta[idx]}")
+        else:
+            diff_indices = np.where(theta != theta_ref)[0]
+            for idx in diff_indices:
+                print(f"Process {comm.Get_rank()} difference at index {idx}: theta_ref={theta_ref[idx]}, theta={theta[idx]}")
+        raise ValueError(
+            f"Process {comm.Get_rank()} has a different theta than the reference process."
+            f" Expected: {theta_ref}, but got:  {theta}. diff = {norm_diff:.4e}"
+            f""
+        )
diff --git a/src/dalia/utils/reparametrizations.py b/src/dalia/utils/reparametrizations.py
new file mode 100644
index 00000000..ade3d130
--- /dev/null
+++ b/src/dalia/utils/reparametrizations.py
@@ -0,0 +1,348 @@
+from scipy.stats import norm
+
+from dalia import NDArray, xp
+
+
+"""
+Computing Quantiles for Transformed Distributions
+
+Theory:
+If X ~ f(x) and Y = φ(X), then to find quantiles of Y:
+1. For a given probability p, find q_p such that P(Y ≤ q_p) = p
+2. This is equivalent to P(φ(X) ≤ q_p) = p
+3. If φ is monotone increasing:     x_original = np.linspace(original_lower, original_upper, 1000)
+    x_internal_array = np.array([transform_func(x, "forward") for x in x_original])
+    pdf_original = np.array([
+        compute_transformed_pdf(mean_internal, std_internal**2, x_int, transform_func) 
+        for x_int in x_internal_array
+    ]) ≤ φ⁻¹(q_p)) = p
+4. So φ⁻¹(q_p) = F_X⁻¹(p), where F_X⁻¹ is the quantile function of X
+5. Therefore: q_p = φ(F_X⁻¹(p))
+
+Key insight: Quantiles transform directly through φ!
+"""
+
+
+def compute_transformed_quantiles(mean_internal, var_internal, percentiles, transform):
+    """
+    Compute quantiles for a transformed distribution
+    
+    Parameters:
+    - original_dist_params: (mean, std) for the internal/transformed distribution
+    - percentiles: array of probability values (0, 1)
+    - transform: TransformationFunction object
+    
+    Returns:
+    - quantiles in original scale
+    """
+    
+    # Step 1: Compute quantiles in internal scale
+    internal_quantiles = norm.ppf(percentiles, loc=mean_internal, scale=var_internal**0.5)
+    
+    # Step 2: Transform back to original scale
+    # If φ: original → internal, then original quantiles = φ⁻¹(internal quantiles)
+    original_quantiles = transform(internal_quantiles, direction='backward')
+    
+    return original_quantiles
+
+def compute_transformed_pdf(mean_internal, var_internal, x_internal, transform):
+    """
+    Compute PDF of transformed distribution using change of variables
+    
+    If Y = φ(X), then f_Y(y) = f_X(φ⁻¹(y)) * |dφ⁻¹/dy|
+    But we want f_X(x) where x is in original scale, so:
+    f_X(x) = f_Y(φ(x)) * |dφ/dx|
+    """ 
+        
+    # PDF in internal scale
+    pdf_internal = norm.pdf(x_internal, loc=mean_internal, scale=var_internal**0.5)
+    
+    # Jacobian: derivative of transformation
+    x_original = transform(x_internal, direction='backward')
+    jacobian = xp.abs(transform(x_original, direction='forward_jacobian'))
+    
+    # PDF in original scale
+    pdf_original = pdf_internal * jacobian
+    
+    return pdf_original
+
+# Automatic bound calculation based on 3 standard deviations in internal scale
+def compute_bounds(mean_internal, var_internal, transform, n_std=3):
+    """
+    Compute plotting bounds based on n standard deviations in internal scale
+    """
+    # Internal scale bounds (±n standard deviations)
+    internal_lower = mean_internal - n_std * var_internal**0.5
+    internal_upper = mean_internal + n_std * var_internal**0.5
+
+    # Transform to original scale
+    original_lower = transform(internal_lower, direction='backward')
+    original_upper = transform(internal_upper, direction='backward')
+
+    return (internal_lower, internal_upper), (original_lower, original_upper)
+
+
+if __name__ == "__main__":
+    """
+    Test reparametrization functions using a dummy gamma prior hyperparameter class.
+    
+    This test demonstrates how the reparametrization functions work with 
+    transformation functions that have forward, backward, and jacobian directions.
+    """
+    import numpy as np
+    import matplotlib.pyplot as plt
+    
+    # Dummy Gamma Prior Hyperparameter Class
+    class DummyGammaPrior:
+        """
+        Dummy implementation of gamma prior rescaling for testing purposes.
+        Implements log transformation: forward = log(x), backward = exp(x)
+        """
+        def rescale_hyperparameters_to_internal(self, theta, direction):
+            """Log transformation between positive (external) and unconstrained (internal) space"""
+            if direction == "forward":
+                return np.log(theta)  # theta -> log(theta)
+            elif direction == "backward": 
+                return np.exp(theta)  # log(theta) -> theta
+            elif direction == "forward_jacobian":
+                return 1.0 / theta   # d(log(theta))/d(theta) = 1/theta
+            elif direction == "backward_jacobian":
+                return theta         # d(exp(theta))/d(theta) = exp(theta) = theta
+            else:
+                raise ValueError(f"Unknown direction: {direction}")
+    
+    print("=" * 80)
+    print("Testing Reparametrization Functions with Dummy Gamma Prior")
+    print("=" * 80)
+    
+    # Create dummy gamma prior instance
+    gamma_prior = DummyGammaPrior()
+    
+    # Define transform function compatible with reparametrization functions
+    def transform_func(x, direction):
+        return gamma_prior.rescale_hyperparameters_to_internal(x, direction)
+    
+    # Test parameters (internal space: log-normal distribution)
+    mean_internal = 0.5   # mean of log(theta)
+    std_internal = 0.8    # std of log(theta)
+    
+    print(f"Internal distribution parameters:")
+    print(f"  Mean (log scale): {mean_internal:.3f}")
+    print(f"  Std (log scale):  {std_internal:.3f}")
+    print()
+    
+    # Test 1: Transform validation
+    print("1. Testing transformation consistency:")
+    test_values = [0.1, 0.5, 1.0, 2.0, 5.0, 10.0]
+    
+    print("   Original -> Internal -> Original (round-trip test)")
+    for theta in test_values:
+        # Forward transformation
+        log_theta = transform_func(theta, "forward")
+        
+        # Backward transformation  
+        theta_recovered = transform_func(log_theta, "backward")
+        
+        # Check error
+        error = abs(theta - theta_recovered)
+        
+        print(f"   {theta:5.1f} -> {log_theta:6.3f} -> {theta_recovered:6.3f}, "
+              f"error = {error:.2e}")
+    print()
+    
+    # Test 2: Jacobian validation using finite differences
+    print("2. Testing Jacobian accuracy (forward direction):")
+    
+    test_theta_vals = [0.5, 1.0, 2.0, 3.0]
+    eps = 1e-8
+    
+    print("   θ     | Analytical | Numerical  | Error")
+    print("   ------|------------|------------|----------")
+    
+    for theta in test_theta_vals:
+        # Analytical jacobian
+        jac_analytical = transform_func(theta, "forward_jacobian")
+        
+        # Numerical jacobian using finite differences
+        f_plus = transform_func(theta + eps, "forward")
+        f_minus = transform_func(theta - eps, "forward")
+        jac_numerical = (f_plus - f_minus) / (2 * eps)
+        
+        error = abs(jac_analytical - jac_numerical)
+        
+        print(f"   {theta:4.1f} | {jac_analytical:10.6f} | {jac_numerical:10.6f} | {error:.2e}")
+    print()
+    
+    # Test 3: Quantile computation
+    print("3. Testing quantile computation:")
+    
+    percentiles = np.array([0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975])
+    original_quantiles = compute_transformed_quantiles(
+        mean_internal, std_internal**0.5, percentiles, transform_func
+    )
+    
+    print("   Percentile | Quantile (original scale)")
+    print("   -----------|-----------------------")
+    for p, q in zip(percentiles, original_quantiles):
+        print(f"   {p:8.2f}   | {q:18.6f}")
+    print()
+    
+    # Test 4: PDF computation and validation
+    print("4. Testing PDF computation:")
+    
+    # Compute bounds for plotting
+    (internal_lower, internal_upper), (original_lower, original_upper) = compute_bounds(
+        mean_internal, std_internal, transform_func, n_std=3
+    )
+    
+    print(f"   Internal bounds: [{internal_lower:.3f}, {internal_upper:.3f}]")
+    print(f"   Original bounds: [{original_lower:.3f}, {original_upper:.3f}]")
+    
+    # Test PDF at specific points
+    test_x_original = np.array([0.5, 1.0, 2.0, 3.0, 5.0])
+    
+    print("   x (orig) | PDF (orig) | log(x)   | PDF (int)")
+    print("   ---------|------------|----------|----------")
+    
+    for x in test_x_original:
+        x_internal = transform_func(x, "forward")
+        pdf_orig = compute_transformed_pdf(mean_internal, std_internal**2, x_internal, transform_func)
+        pdf_internal = norm.pdf(x_internal, loc=mean_internal, scale=std_internal)
+        
+        print(f"   {x:7.2f}  | {pdf_orig:10.6f} | {x_internal:8.3f} | {pdf_internal:8.6f}")
+    print()
+    
+    # Test 5: Analytical validation for log-normal distribution
+    print("5. Analytical validation (log-normal distribution):")
+    
+    # For log-normal distribution, we can compute analytical moments
+    mu = mean_internal
+    sigma = std_internal
+    
+    # Analytical log-normal statistics
+    analytical_mean = np.exp(mu + sigma**2/2)
+    analytical_var = (np.exp(sigma**2) - 1) * np.exp(2*mu + sigma**2)
+    analytical_std = np.sqrt(analytical_var)
+    
+    # Numerical verification using quantiles
+    # Mean ≈ 50th percentile for log-normal (approximately)
+    median_quantile = compute_transformed_quantiles(
+        mean_internal, std_internal**0.5, np.array([0.5]), transform_func
+    )[0]
+    
+    print(f"   Analytical mean: {analytical_mean:.6f}")
+    print(f"   Analytical std:  {analytical_std:.6f}")
+    print(f"   Median quantile: {median_quantile:.6f}")
+    print(f"   Mean/Median ratio: {analytical_mean/median_quantile:.6f}")
+    print("   (Should be > 1 for log-normal due to skewness)")
+    print()
+    
+    # Test 6: PDF integration check (numerical verification)
+    print("6. PDF integration check:")
+    
+    # Create fine grid for integration
+    x_grid = np.linspace(original_lower, original_upper, 1000)
+    x_internal_grid = np.array([transform_func(x, "forward") for x in x_grid])
+    pdf_values = np.array([
+        compute_transformed_pdf(mean_internal, std_internal**2, x_int, transform_func) 
+        for x_int in x_internal_grid
+    ])
+    
+    # Numerical integration using trapezoidal rule
+    dx = x_grid[1] - x_grid[0]
+    integral = np.trapz(pdf_values, dx=dx)
+    
+    print(f"   Numerical integral of PDF: {integral:.6f}")
+    print(f"   Should be close to 1.0, error: {abs(1.0 - integral):.6f}")
+    print()
+    
+    # Test 7: Bounds computation for different n_std values
+    print("7. Testing bounds computation:")
+    
+    for n_std in [1, 2, 3, 4]:
+        (int_lower, int_upper), (orig_lower, orig_upper) = compute_bounds(
+            mean_internal, std_internal, transform_func, n_std=n_std
+        )
+        
+        print(f"   {n_std}σ bounds:")
+        print(f"     Internal: [{int_lower:7.3f}, {int_upper:7.3f}]")
+        print(f"     Original: [{orig_lower:7.3f}, {orig_upper:7.3f}]")
+    print()
+    
+    # Test 8: Edge case handling
+    print("8. Testing edge cases:")
+    
+    # Test very small values
+    small_vals = [1e-6, 1e-4, 1e-2]
+    print("   Small values transformation:")
+    for val in small_vals:
+        try:
+            internal = transform_func(val, "forward")
+            recovered = transform_func(internal, "backward")
+            error = abs(val - recovered)
+            print(f"     {val:.2e} -> {internal:8.3f} -> {recovered:.2e}, error = {error:.2e}")
+        except Exception as e:
+            print(f"     {val:.2e} -> Error: {e}")
+    
+    # Test large values
+    large_vals = [1e2, 1e4, 1e6]
+    print("   Large values transformation:")
+    for val in large_vals:
+        try:
+            internal = transform_func(val, "forward")
+            recovered = transform_func(internal, "backward")
+            rel_error = abs(val - recovered) / val
+            print(f"     {val:.2e} -> {internal:8.3f} -> {recovered:.2e}, rel_error = {rel_error:.2e}")
+        except Exception as e:
+            print(f"     {val:.2e} -> Error: {e}")
+    print()
+    
+    # Test 9: Plotting PDFs in both scales
+    print("9. Plotting PDFs in internal and original scales:")
+        
+    # Plot 1: PDF in internal scale (log-scale, normal distribution)
+    x_internal = np.linspace(internal_lower, internal_upper, 500)
+    pdf_internal = norm.pdf(x_internal, loc=mean_internal, scale=std_internal)
+    # Create visualization
+    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
+
+    ax1.plot(x_internal, pdf_internal, 'b-', linewidth=2, label='Internal PDF')
+    internal_quantiles = norm.ppf(percentiles, loc=mean_internal, scale=std_internal)
+    for i, (p, q) in enumerate(zip(percentiles, internal_quantiles)):
+        color = 'red' if p in [0.025, 0.975] else 'orange'
+        ax1.axvline(q, color=color, linestyle='--', alpha=0.7)
+        if i % 2 == 0:
+            ax1.text(q, max(pdf_internal)*0.8, f'{p*100:.1f}%', rotation=90, ha='right', va='top')
+
+    ax1.set_title(f'Internal Scale: N(mean = {mean_internal}, std = {std_internal:.3f})')
+    ax1.set_xlabel('θ (internal)')
+    ax1.set_ylabel('PDF')
+    ax1.set_xlim(x_internal[0], x_internal[-1])
+    ax1.grid(True, alpha=0.3)
+    ax1.legend()
+    
+    # Plot 2: Original distribution with quantiles
+    x_original = np.linspace(orig_lower, orig_upper, 1000)
+    x_internal = transform_func(x_original, "forward")
+    pdf_original = compute_transformed_pdf(mean_internal, std_internal**2, x_internal, transform_func)
+
+    ax2.plot(x_original, pdf_original, 'g-', linewidth=2, label='Original PDF')
+    for i, (p, q) in enumerate(zip(percentiles, original_quantiles)):
+        color = 'red' if p in [0.025, 0.975] else 'orange'
+        ax2.axvline(q, color=color, linestyle='--', alpha=0.7)
+        if i % 2 == 0:
+            ax2.text(q, max(pdf_original)*0.8, f'{p*100:.1f}%', rotation=90, ha='right', va='top')
+
+    ax2.set_title('Original Scale: Log-Normal Distribution')
+    ax2.set_xlabel('θ_outer (original)')
+    ax2.set_ylabel('PDF')
+    ax2.set_xlim(orig_lower, 15)
+    ax2.grid(True, alpha=0.3)
+    ax2.legend()
+    plt.tight_layout()
+    plt.show()
+        
+    print("=" * 80)
+    print("All reparametrization functions are working correctly with")
+    print("the gamma prior rescaling transformation!")
+    print("=" * 80)
\ No newline at end of file