Update Bayesian_for_mPWP_PlioMIP.py

Bayesian_for_mPWP_PlioMIP.py

@@ -104,20 +104,39 @@ pliomip2y = [0.9211941,2.1174774, 1.3736095]
 ## MCMC model
 
 with Model() as model: # model specifications in PyMC3 are wrapped in a with-statement
-    # Define priors
+    # Define prior means, sd or precision (see conjugate for precision)
+    
+    mn_intercept = 0.0
+    mn_slope = 1.0
+    sd_reg = 1.0
+    
+    precis_intercept = np.sqrt(0.5)
+    precis_slope = np.sqrt(0.5)
+        
+    # Define priors: Non-conjugate approach
+    
     sigma = HalfCauchy('Sigma', beta=5, testval=1.)
-    intercept = Normal('Intercept', 0, sd=1)
-    x_coeff = Normal('x', 1, sd=1)
+    intercept = Normal('Intercept', mn_intercept, sd=sd_reg)
+    x_coeff = Normal('x', mn_slope, sd=sd_reg)
+
+    # Priors for a Normal-Inverse Gamma conjugate approach. 
+    # In conjugate approach, priors on intercept and slope depends on a scaled sigma
+    # and need to be defined that way to match mathematical equations.
+    # This part should be mainly used for comparison / check with the conjugate approach code
+    
+    # sigma = InverseGamma('Sigma',alpha=0.5,beta=0.5)
+    # intercept = Normal('Intercept', mn_intercept, sd=sigma/precis_intercept)
+    # x_coeff = Normal('x', mn_slope, sd=sigma/precis_slope)
 
     # Define likelihood
+    
     likelihood = Normal('y', mu=intercept + x_coeff * x, 
                         sd=sigma, observed=y)
     
-    # 4 chains in parallel
+    # Inference! 4 jobs in parallel (convergence check)
     # By default, the sampling method is NUTS
-    # The following line will not work with PyMC3 older than 3.8. If you use an
-    # older version, replace "cores=4" by "njobs=4"
-    trace = sample(progressbar=False, draws=100000, cores=4)
+    
+    trace = sample(progressbar=True, draws=30000, cores=2, tune=2000)
 
 # Extract the data of the trace
 values_x = trace['x']