Bayesian analysis of NPD data

Bayesian analysis of NPD data#

This notebook continues the NPD refinement tutorial. In that tutorial, you built a two-phase model containing La₀.₅Ba₀.₅CoO₃ (LBCO) and a small Si impurity, refined the model against simulated time-of-flight neutron data, and saved the project.

A conventional least-squares refinement gives a best-fit point and uncertainty estimates based on the local shape of the objective function. Bayesian analysis instead uses Markov chain Monte Carlo (MCMC) sampling to explore the joint posterior distribution of the free parameters. This allows us to investigate questions such as:

Where is each posterior distribution centered, and how does that compare with the deterministic best-fit value?
How broad are their credible intervals?
Do the local parameter correlations from the deterministic fit persist across the posterior, and are any relationships nonlinear or asymmetric?
How does parameter uncertainty propagate into the calculated diffraction pattern?

The deterministic fit already provides a correlation matrix derived from the local covariance estimate near the best-fit point. Bayesian sampling extends this view by exploring how parameters vary together throughout the joint posterior. It can therefore reveal curved, asymmetric, bounded, or multimodal relationships that one local correlation coefficient cannot describe.

This notebook follows the same teaching structure as the refinement tutorial:

Introduction: a complete, guided MCMC analysis of the refined LBCO+Si model.
Exercise: repeat the analysis after fixing one member of a strongly correlated peak-profile pair, then compare the results.

We will use the DREAM sampler provided by the bumps (dream) minimizer. The short chains used here keep the tutorial practical. A scientific analysis requires longer chains and careful convergence checks.

🛠️ Import Library#

import easydiffraction as edi

📘 Introduction: MCMC Analysis of LBCO+Si#

We will first work through the complete Bayesian workflow without exercises. This reference analysis samples the same seven scientific parameters that remained free at the end of the refinement tutorial:

the LBCO lattice parameter,
the LBCO and Si scale factors, and
four peak-profile parameters.

The resulting project will be called project_1, following the convention used for the introductory Si fit in the previous tutorial.

📂 Load the Refined Project#

Rather than rebuilding the experiment and structures, load the project saved at the end of the refinement tutorial. This restores the measured data, both structures, the refined values, and the free-parameter settings.

The following cell reads the refined project from the local projects directory. Complete the NPD refinement tutorial first so that projects/fitting-exercise-si-lbco-main is available.

📖 See documentation for more details about loading a saved project.

refinement_project_dir = 'projects/fitting-exercise-si-lbco-main'
project_1 = edi.Project.load(refinement_project_dir)

Save a copy under a new name before changing the analysis. This keeps the deterministic refinement unchanged and gives the sampler its own location for posterior summaries and chain data.

project_1.metadata.title = 'Reference Bayesian Analysis of LBCO with Si Impurity'
project_1.metadata.description = (
    'Reference MCMC analysis of a two-phase LBCO and Si powder diffraction model.'
)
project_1.save_as(dir_path='projects/exercise-bayesian-si-lbco-reference')

Saving project 📦 'main' to 'projects/exercise-bayesian-si-lbco-reference'

├── 📄 project.edi
├── 📁 structures/
│   └── 📄 lbco.edi
│   └── 📄 si.edi
├── 📁 experiments/
│   └── 📄 sim_lbco.edi
├── 📁 analysis/
│   └── 📄 analysis.edi
└── 📁 reports/
    └── 📄 main.html

Confirm that the restored project contains the LBCO and Si structures and the sim_lbco experiment.

project_1.structures.show_names()
project_1.experiments.show_names()

Defined structures 🧩

['lbco', 'si']

Defined experiments 🔬

['sim_lbco']

project_1.display.structure(struct_name='lbco')
project_1.display.structure(struct_name='si')

Structure 🧩 'lbco' (Atom view type: 'covalent')

Loading plot…

drag = rotate
wheel = zoom
right-drag = pan

Structure 🧩 'si' (Atom view type: 'covalent')

Loading plot…

drag = rotate
wheel = zoom
right-drag = pan

project_1.display.pattern(expt_name='sim_lbco')

🎯 Choose the Sampled Parameters#

MCMC varies every free parameter simultaneously. Before sampling, it is worth asking which parameters are needed to answer the scientific question.

First, display all free parameters in the refined project.

project_1.display.parameters.free()

Free parameters for both structures (🧩 data blocks) and experiments (🔬 data blocks)

	datablock	category	entry	parameter	value	uncertainty	min	max	units
1	lbco	cell		length_a	3.89191	0.00028	-inf	inf	Å
2	sim_lbco	linked_structure	lbco	scale	43.65000	0.39000	-inf	inf
3	sim_lbco	linked_structure	si	scale	0.32100	0.02200	-inf	inf
4	sim_lbco	peak		rise_alpha_1	0.01232	0.00017	-inf	inf	μs/Å
5	sim_lbco	peak		decay_beta_1	0.01147	0.00047	-inf	inf	μs/Å
6	sim_lbco	peak		broad_gauss_sigma_1	-78400.00000	2900.00000	-inf	inf	μs/Å
7	sim_lbco	peak		broad_gauss_sigma_2	24500.00000	1700.00000	-inf	inf	μs²/Å²
8	sim_lbco	background	1	intensity	0.28600	0.03000	-inf	inf
9	sim_lbco	background	2	intensity	0.25650	0.00650	-inf	inf
10	sim_lbco	background	3	intensity	0.25550	0.00570	-inf	inf
11	sim_lbco	background	4	intensity	0.23300	0.00530	-inf	inf
12	sim_lbco	background	5	intensity	0.23900	0.00390	-inf	inf
13	sim_lbco	background	6	intensity	0.22040	0.00470	-inf	inf
14	sim_lbco	background	7	intensity	0.23500	0.01600	-inf	inf

The seven line-segment background intensities were refined in the previous tutorial. If we left them free, this analysis would sample 14 parameters instead of 7. The larger space would require more model evaluations, take longer to explore, and generally need a longer chain to mix well. Background parameters can also correlate with phase scales and broad peaks, making the posterior harder to interpret.

We therefore fix the background at its refined values. This is a practical simplification for a short teaching example, not a universal rule. Fixing the background means that its uncertainty is not propagated into the final credible intervals. In a complete scientific analysis, you should sample relevant background parameters or otherwise account for their uncertainty when it can affect the result.

experiment_1 = project_1.experiments['sim_lbco']

for line_segment in experiment_1.background:
    line_segment.intensity.free = False

Display the remaining free parameters to verify the sampling problem.

project_1.display.parameters.free()

Free parameters for both structures (🧩 data blocks) and experiments (🔬 data blocks)

	datablock	category	entry	parameter	value	uncertainty	min	max	units
1	lbco	cell		length_a	3.89191	0.00028	-inf	inf	Å
2	sim_lbco	linked_structure	lbco	scale	43.65000	0.39000	-inf	inf
3	sim_lbco	linked_structure	si	scale	0.32100	0.02200	-inf	inf
4	sim_lbco	peak		rise_alpha_1	0.01232	0.00017	-inf	inf	μs/Å
5	sim_lbco	peak		decay_beta_1	0.01147	0.00047	-inf	inf	μs/Å
6	sim_lbco	peak		broad_gauss_sigma_1	-78400.00000	2900.00000	-inf	inf	μs/Å
7	sim_lbco	peak		broad_gauss_sigma_2	24500.00000	1700.00000	-inf	inf	μs²/Å²

🚀 Run a Local Refinement#

Fixing the background changes the optimization problem slightly. We first run a fast local least-squares refinement to update the best-fit values and their uncertainty estimates. These values will provide the starting point and finite bounds for MCMC.

We use the BUMPS Levenberg-Marquardt minimizer, bumps (lm), so the local and Bayesian stages use the same minimizer library.

📖 See documentation for more details about the available minimizers.

project_1.analysis.minimizer.show_supported()

Minimizer types

		Type	Description
1		bumps	BUMPS library using the default Levenberg-Marquardt method
2		bumps (amoeba)	BUMPS library with Nelder-Mead simplex method
3		bumps (de)	BUMPS library with differential evolution method
4		bumps (dream)	BUMPS library with DREAM Bayesian sampling
5		bumps (lm)	BUMPS library with Levenberg-Marquardt method
6		dfols	DFO-LS library for derivative-free least-squares optimization
7		emcee	emcee affine-invariant ensemble Bayesian sampling
8		lmfit	LMFIT library using the default Levenberg-Marquardt method
9		lmfit (least_squares)	LMFIT library with SciPy's trust region reflective algorithm
10	*	lmfit (leastsq)	LMFIT library with Levenberg-Marquardt least squares method

project_1.analysis.minimizer.type = 'bumps (lm)'

Current minimizer changed to

bumps (lm)

project_1.analysis.fit()
project_1.display.fit.results()

Standard fitting

📋 Using experiment 🔬 'sim_lbco' for 'single' fitting

🚀 Starting fit process with 'bumps (lm)'...

📈 Goodness-of-fit progress:

	iteration	time (s)	χ²	change / status
1	1	0.02	1.24
2	27	0.76	1.24

🏆 Best goodness-of-fit (reduced χ²) is 1.24 at iteration 27

✅ Fitting complete.

⚙️ Settings used:

	Name	Value	Description
1	max_iterations	1000	Maximum solver iterations.

📋 Least-squares fit results:

	Metric	Value
1	🧪 Minimizer	bumps (lm)
2	✅ Overall status	success
3	⏱️ Fitting time (seconds)	0.76
4	📏 Goodness-of-fit (reduced χ²)	1.24
5	📏 R-factor (Rf, %)	6.89
6	📏 R-factor squared (Rf², %)	5.58
7	📏 Weighted R-factor (wR, %)	4.76

📈 Refined parameters:

	datablock	category	entry	parameter	units	start	value	s.u.	change
1	lbco	cell		length_a	Å	3.8919	3.8919	0.0003	0.00 % ↓
2	sim_lbco	linked_structure	lbco	scale		43.6500	43.6497	0.3770	0.00 % ↓
3	sim_lbco	linked_structure	si	scale		0.3210	0.3211	0.0209	0.04 % ↑
4	sim_lbco	peak		rise_alpha_1	μs/Å	0.0123	0.0123	0.0002	0.05 % ↓
5	sim_lbco	peak		decay_beta_1	μs/Å	0.0115	0.0115	0.0005	0.04 % ↑
6	sim_lbco	peak		broad_gauss_sigma_1	μs/Å	-78400.0000	-78430.2812	2835.8370	0.04 % ↑
7	sim_lbco	peak		broad_gauss_sigma_2	μs²/Å²	24500.0000	24485.3821	1722.2569	0.06 % ↓

  • start      = parameter value before refinement
  • value      = refined value from least-squares minimization
  • s.u.       = standard uncertainty (one sigma), from the covariance matrix
  • change     = relative change from start, in %; ↑ = increase, ↓ = decrease

🔗 Understand the Local Correlation Chart#

The correlation chart summarizes how pairs of refined parameters can change together near the least-squares optimum. Each off-diagonal value is a correlation coefficient between -1 and +1:

a value near +1 means that the two parameters tend to increase together;
a value near -1 means that one tends to decrease when the other increases; and
a value near 0 means that there is little linear relationship.

Correlation does not mean that one parameter physically causes the other. It means that the measured pattern cannot distinguish their effects independently within this model.

project_1.display.fit.correlations(max_parameters=5)

The matrix shows only one triangular half because the other half would contain the same values in reverse order. With max_parameters=5, EasyDiffraction filters the chart to keep the strongest relationships readable. In an interactive Jupyter display, hover over a colored cell to see the two full parameter names and the numerical correlation coefficient.

The strongest relationship is typically the negative correlation between broad_gauss_sigma_1 and broad_gauss_sigma_2, close to -0.94. Both parameters contribute to the d-spacing dependence of the Gaussian TOF peak width. An increase in one can be partly compensated by a decrease in the other while producing a similar calculated peak shape.

This chart comes from the local covariance estimate, so it describes only the neighborhood around the best-fit point. MCMC will show whether the relationship remains linear and approximately elliptical across a wider region of parameter space.

🎲 Define the Sampling Region#

DREAM requires finite lower and upper bounds for every sampled parameter. Here, we derive them from the uncertainty estimated by the local fit. set_fit_bounds_from_uncertainty() places each bound four estimated standard uncertainties from the current value by default, while respecting any physical parameter limits.

These bounds act as bounded prior support in this example. They must be inspected rather than accepted blindly: narrow bounds can truncate the posterior, while unnecessarily wide bounds can make sampling less efficient.

📖 See documentation for more details about parameter values, uncertainties, and fit bounds.

for param in project_1.free_parameters:
    param.set_fit_bounds_from_uncertainty()

project_1.display.parameters.free()

Free parameters for both structures (🧩 data blocks) and experiments (🔬 data blocks)

	datablock	category	entry	parameter	value	uncertainty	min	max	units
1	lbco	cell		length_a	3.89191	0.00028	3.89079	3.89303	Å
2	sim_lbco	linked_structure	lbco	scale	43.64972	0.37698	42.14178	45.15766
3	sim_lbco	linked_structure	si	scale	0.32113	0.02085	0.23772	0.40454
4	sim_lbco	peak		rise_alpha_1	0.01231	0.00017	0.01165	0.01298	μs/Å
5	sim_lbco	peak		decay_beta_1	0.01148	0.00046	0.00965	0.01330	μs/Å
6	sim_lbco	peak		broad_gauss_sigma_1	-78430.28123	2835.83698	-89773.62915	-67086.93331	μs/Å
7	sim_lbco	peak		broad_gauss_sigma_2	24485.38213	1722.25685	17596.35472	31374.40955	μs²/Å²

🎲 Run DREAM Sampling#

In Bayesian analysis, the posterior combines the likelihood of the measured data with prior assumptions. MCMC constructs chains of correlated samples whose long-run distribution approximates that posterior.

DREAM uses multiple chains and differential-evolution proposals to explore correlated parameter spaces. This makes it useful for the strong peak-profile correlation observed above.

📖 See documentation for more details about Bayesian minimizers and posterior displays.

project_1.analysis.minimizer.type = 'bumps (dream)'
project_1.analysis.minimizer.sampling_steps = 300  # lower than the default 3000
project_1.analysis.minimizer.burn_in_steps = 60  # lower than the default 600
project_1.analysis.minimizer.random_seed = 42

⚠️ Switching minimizer type removes these settings:                                                                               
   • max_iterations                                                                                                               

⚠️ Switching minimizer type adds these settings with defaults:                                                                    
   • burn_in_steps=600                                                                                                            
   • initialization_method='latin_hypercube'                                                                                      
   • parallel_workers=0                                                                                                           
   • population_size=4                                                                                                            
   • random_seed=None                                                                                                             
   • sampling_steps=3000                                                                                                          
   • thinning_interval=1                                                                                                          

Current minimizer changed to

bumps (dream)

Burn-in samples allow the chains to move away from their initial positions before the retained posterior is summarized. The fixed seed makes the tutorial output reproducible.

project_1.analysis.fit()

Standard fitting

📋 Using experiment 🔬 'sim_lbco' for 'single' fitting

🚀 Starting fit process with 'bumps (dream)'...

📈 Bayesian sampling progress:

	step	progress	time (s)	log posterior	phase
1			0.60	-435.88	pre-processing
2	21/360	5.8%	5.67	-275.53	burn-in
3	40/360	11.1%	10.57	-268.83	burn-in
4	60/360	16.7%	15.72	-267.41	burn-in
5	61/360	16.9%	16.03	-267.29	sampling
6	76/360	21.1%	19.93	-267.10	sampling
7	91/360	25.3%	24.20	-266.65	sampling
8	106/360	29.4%	28.30	-266.10	sampling
9	121/360	33.6%	32.37	-266.35	sampling
10	136/360	37.8%	36.35	-266.12	sampling
11	151/360	41.9%	40.40	-267.18	sampling
12	166/360	46.1%	44.43	-266.73	sampling
13	181/360	50.3%	48.78	-266.95	sampling
14	196/360	54.4%	52.76	-266.56	sampling
15	211/360	58.6%	57.02	-266.54	sampling
16	226/360	62.8%	61.04	-266.85	sampling
17	241/360	66.9%	65.12	-266.57	sampling
18	256/360	71.1%	68.77	-266.49	sampling
19	271/360	75.3%	73.22	-266.52	sampling
20	286/360	79.4%	77.25	-266.56	sampling
21	301/360	83.6%	81.37	-266.17	sampling
22	316/360	87.8%	85.18	-266.81	sampling
23	331/360	91.9%	89.30	-266.49	sampling
24	346/360	96.1%	93.24	-266.30	sampling
25	360/360	100.0%	97.17	-267.24	sampling
26			124.07		post-processing

✅ Bayesian sampling complete.

⚠️ Convergence diagnostics indicate the posterior may be poorly mixed.

📋 Understand the Bayesian Fit Summary#

The result table now reports posterior medians and 95% credible intervals in addition to fit-quality metrics. A 95% credible interval is the interval containing 95% of the retained marginal posterior samples under this model and bounded sampling setup.

The table also reports convergence diagnostics:

r-hat compares variation within and between chains. Values close to 1 are desirable; EasyDiffraction uses r-hat <= 1.01 as its displayed convergence criterion.
ess bulk is the effective sample size after accounting for autocorrelation. A larger value means that the chain contains more independent information; the display recommends at least 400.

project_1.display.fit.results()

⚙️ Settings used:

	Name	Value	Description
1	sampling_steps	300	Total sampler iterations per chain.
2	burn_in_steps	60	Sampler iterations discarded as warm-up.
3	thinning_interval	1	Sampler thinning interval.
4	population_size	4	Number of chains or walkers.
5	parallel_workers	0	Worker count; 0 uses all available CPUs.
6	initialization_method	latin_hypercube	Sampler initialization method.
7	random_seed	42	Random seed; None uses a system-derived seed.

📋 Bayesian fit results:

	Metric	Value
1	🧪 Sampler	bumps (dream)
2	❌ Overall status	failed
3	💬 Engine message	DREAM sampling completed
4	⏱️ Fitting time (seconds)	124.07
5	📏 Goodness-of-fit (reduced χ²)	1.24
6	📏 R-factor (Rf, %)	6.89
7	📏 R-factor squared (Rf², %)	5.58
8	📏 Weighted R-factor (wR, %)	4.76
9	📉 Best log-posterior	-263.03
10	📊 Convergence status	failed
11	📊 Max r-hat	1.180
12	📊 Min ess bulk	296.603
13	📊 Draws per chain	300
14	📊 Chains	28

📈 Committed parameters:

	datablock	category	entry	parameter	units	start	value	s.u.	change
1	lbco	cell		length_a	Å	3.8919	3.8919	0.0003	0.00 % ↓
2	sim_lbco	linked_structure	lbco	scale		43.6497	43.6497	0.3386	0.00 % ↓
3	sim_lbco	linked_structure	si	scale		0.3211	0.3211	0.0187	0.00 % ↓
4	sim_lbco	peak		rise_alpha_1	μs/Å	0.0123	0.0123	0.0002	0.00 % ↓
5	sim_lbco	peak		decay_beta_1	μs/Å	0.0115	0.0115	0.0004	0.00 % ↓
6	sim_lbco	peak		broad_gauss_sigma_1	μs/Å	-78430.2812	-78430.2812	2557.9102	0.00 % ↓
7	sim_lbco	peak		broad_gauss_sigma_2	μs²/Å²	24485.3821	24485.3821	1585.8085	0.00 % ↓

  • start      = parameter value before sampling
  • value      = estimate written back to the project (best posterior sample)
  • s.u.       = standard uncertainty (one sigma), posterior standard deviation
  • change     = relative change from start, in %; ↑ = increase, ↓ = decrease

📊 Posterior distribution:

	datablock	category	entry	parameter	units	median	95% CI	r-hat	ess bulk
1	lbco	cell		length_a	Å	3.8919	[3.8914, 3.8924]	1.100	296.6
2	sim_lbco	linked_structure	lbco	scale		43.6910	[43.0364, 44.3511]	1.113	355.7
3	sim_lbco	linked_structure	si	scale		0.3221	[0.2826, 0.3578]	1.128	347.1
4	sim_lbco	peak		rise_alpha_1	μs/Å	0.0123	[0.0120, 0.0127]	1.180	312.4
5	sim_lbco	peak		decay_beta_1	μs/Å	0.0115	[0.0107, 0.0124]	1.138	322.1
6	sim_lbco	peak		broad_gauss_sigma_1	μs/Å	-78701.7593	[-84008.4869, -73831.7126]	1.105	312.5
7	sim_lbco	peak		broad_gauss_sigma_2	μs²/Å²	24677.0434	[21778.4041, 27996.7606]	1.110	318.5

  • median       = 50th percentile of the marginal posterior
  • 95% CI       = 95% credible interval (2.5%-97.5%, asymmetric)
  • r-hat        = Gelman-Rubin diagnostic (good convergence: r-hat <= 1.01)
  • ess bulk     = bulk effective sample size (typically >= 400)

⚠️ r-hat > 1.01: Consider longer sampling, better initialization, or reparameterization.
⚠️ ess bulk < 400: Consider longer sampling or reparameterization.

A short teaching run may fail these convergence criteria even when it finishes normally. That means the numerical posterior summaries are provisional. Increase the number of sampling steps and inspect the diagnostics again before drawing scientific conclusions.

🔗 Understand Posterior Correlations#

Calling the same method after MCMC now builds the matrix from posterior samples rather than from the local least-squares covariance estimate. It answers: across the sampled posterior, which parameter pairs vary together, and how strongly?

project_1.display.fit.correlations(max_parameters=5)

Compare this chart with the local chart above. Similar coefficients suggest that the local approximation captured the main linear relationship. A substantial difference can indicate curvature, asymmetry, bounds, or another feature that a local covariance matrix cannot represent.

🗺️ Understand the Posterior Pair Plot#

The pair plot shows more information than a matrix of single coefficients:

Diagonal panels show each parameter’s one-dimensional marginal posterior density. A narrow peak indicates greater precision; a wide, skewed, truncated, or multimodal shape indicates more uncertainty or a more complicated posterior.
Lower-triangle panels show joint posterior samples for pairs of parameters. The contours summarize where the sampled density is concentrated.
A compact, nearly round cloud suggests weak correlation. An elongated upward cloud suggests positive correlation, and an elongated downward cloud suggests negative correlation. Curved or split contours cannot be summarized well by one correlation coefficient.
The upper triangle is intentionally blank because it would duplicate the lower triangle.

project_1.display.posterior.pairs(max_parameters=5)

In Jupyter, the default plotting engine resolves to interactive Plotly. Hover over a diagonal density curve to see the full parameter name, parameter value, and probability density. Hover over a visible sample point in a lower-triangle panel to see the exact values of both parameters for that posterior draw. The shaded contour itself has no hover tooltip; it is a smoothed two-dimensional density guide. Use the Plotly toolbar to zoom, pan, and reset the view, and click legend items to hide or show samples, contours, or marginal densities.

The broad_gauss_sigma_1–broad_gauss_sigma_2 panel should form a long downward band. This reveals the range of compensating parameter combinations, whereas the correlation chart reduces the relationship to one number.

📈 Understand Marginal Posterior Distributions#

A separate distribution plot gives a more detailed view of each diagonal panel. It includes a histogram, a smoothed marginal density, the median, the best posterior sample, and the 95% credible interval. Hover over the traces and interval markers to read their values.

project_1.display.posterior.distribution()

A density pressed against a fit bound warns that the allowed region may be too narrow or that the parameter is poorly identified. A best sample far from the median can occur for a skewed or irregular posterior and is another reason not to summarize MCMC with only one point estimate.

📊 Understand the Posterior-Predictive Plot#

Posterior prediction propagates an evenly spaced subset of retained parameter combinations through the diffraction calculation. The best-posterior-sample curve shows one calculated pattern, while the 95% band shows parameter uncertainty propagated into the pattern.

project_1.display.posterior.predictive(expt_name='sim_lbco')

project_1.display.posterior.predictive(
    expt_name='sim_lbco',
    x_min=80000,
    x_max=81500,
)

Compare the width of the band with the experimental uncertainty and remaining residuals. A narrow band does not prove that the model is correct: systematic deviations can indicate model inadequacy, underestimated experimental uncertainty, or effects not included in the refinement. In the interactive view, hover over the measured and best-sample curves to inspect individual x and intensity values. The shaded 95% band is read from its upper and lower boundaries and does not itself show a hover tooltip.

Save the completed reference Bayesian project. Its MCMC chain and posterior data are stored with the project and remain available for comparison with the exercise below.

project_1.save()

Saving project 📦 'main' to 'projects/exercise-bayesian-si-lbco-reference'

├── 📄 project.edi
├── 📁 structures/
│   └── 📄 lbco.edi
│   └── 📄 si.edi
├── 📁 experiments/
│   └── 📄 sim_lbco.edi
├── 📁 analysis/
│   ├── 📄 analysis.edi
│   └── 📄 mcmc.h5
└── 📁 reports/
    └── 📄 main.html

💪 Exercise: Fix One Correlated Parameter#

The reference analysis showed that broad_gauss_sigma_1 and broad_gauss_sigma_2 are strongly negatively correlated. In this exercise, you will fix broad_gauss_sigma_2 at its refined value and repeat the Bayesian workflow.

Removing one member of the pair reduces the sampling dimension and removes that pair from the correlation and pair plots. This can improve sampling efficiency, but it changes the scientific question: the new posterior is conditional on the chosen fixed value of broad_gauss_sigma_2. Its uncertainty is no longer propagated. The exercise demonstrates the computational and visual effect of fixing a parameter; it does not establish that fixing it is always the correct scientific choice.

📂 Exercise 1: Create a Fresh Project#

Load the original deterministic refinement again as project_2 and save it under a new name. Starting from the deterministic project ensures that the second MCMC run does not reuse the first posterior state or its uncertainty estimates.

Hint:

Use edi.Project.load() with the local refinement_project_dir defined in the introduction.

Solution:

Saving project 📦 'main' to 'projects/exercise-bayesian-si-lbco-main'

├── 📄 project.edi
├── 📁 structures/
│   └── 📄 lbco.edi
│   └── 📄 si.edi
├── 📁 experiments/
│   └── 📄 sim_lbco.edi
├── 📁 analysis/
│   └── 📄 analysis.edi
└── 📁 reports/
    └── 📄 main.html

🎯 Exercise 2: Reduce the Free-Parameter Set#

Exercise 2.1: Fix the Background#

Fix every background intensity for the same computational reason as in the reference analysis.

Hint:

Iterate over project_2.experiments['sim_lbco'].background and set each point’s intensity.free attribute to False.

Solution:

Exercise 2.2: Fix `broad_gauss_sigma_2`#

Remove broad_gauss_sigma_2 from the refined and sampled parameter set, then display the remaining free parameters. How many remain?

Hint:

Set the parameter’s free attribute to False, then call project_2.display.parameters.free().

Solution:

Free parameters for both structures (🧩 data blocks) and experiments (🔬 data blocks)

	datablock	category	entry	parameter	value	uncertainty	min	max	units
1	lbco	cell		length_a	3.89191	0.00028	-inf	inf	Å
2	sim_lbco	linked_structure	lbco	scale	43.65000	0.39000	-inf	inf
3	sim_lbco	linked_structure	si	scale	0.32100	0.02200	-inf	inf
4	sim_lbco	peak		rise_alpha_1	0.01232	0.00017	-inf	inf	μs/Å
5	sim_lbco	peak		decay_beta_1	0.01147	0.00047	-inf	inf	μs/Å
6	sim_lbco	peak		broad_gauss_sigma_1	-78400.00000	2900.00000	-inf	inf	μs/Å

Six parameters remain. The background intensities and broad_gauss_sigma_2 will stay at their refined values throughout the following local fit and MCMC run.

🚀 Exercise 3: Repeat the Local Refinement#

Select bumps (lm), run the fit, display the result table, and inspect the local correlation chart. Is the original strongly correlated pair still present?

Hint:

Repeat the local-refinement sequence from the introduction with project_2.

Solution:

Current minimizer changed to

bumps (lm)

Standard fitting

📋 Using experiment 🔬 'sim_lbco' for 'single' fitting

🚀 Starting fit process with 'bumps (lm)'...

📈 Goodness-of-fit progress:

	iteration	time (s)	χ²	change / status
1	1	0.02	1.24
2	24	0.89	1.24

🏆 Best goodness-of-fit (reduced χ²) is 1.24 at iteration 24

✅ Fitting complete.

⚙️ Settings used:

	Name	Value	Description
1	max_iterations	1000	Maximum solver iterations.

📋 Least-squares fit results:

	Metric	Value
1	🧪 Minimizer	bumps (lm)
2	✅ Overall status	success
3	⏱️ Fitting time (seconds)	0.89
4	📏 Goodness-of-fit (reduced χ²)	1.24
5	📏 R-factor (Rf, %)	6.89
6	📏 R-factor squared (Rf², %)	5.58
7	📏 Weighted R-factor (wR, %)	4.76

📈 Refined parameters:

	datablock	category	entry	parameter	units	start	value	s.u.	change
1	lbco	cell		length_a	Å	3.8919	3.8919	0.0003	0.00 % ↓
2	sim_lbco	linked_structure	lbco	scale		43.6500	43.6498	0.3763	0.00 % ↓
3	sim_lbco	linked_structure	si	scale		0.3210	0.3211	0.0208	0.05 % ↑
4	sim_lbco	peak		rise_alpha_1	μs/Å	0.0123	0.0123	0.0002	0.04 % ↓
5	sim_lbco	peak		decay_beta_1	μs/Å	0.0115	0.0115	0.0004	0.05 % ↑
6	sim_lbco	peak		broad_gauss_sigma_1	μs/Å	-78400.0000	-78452.4995	1074.7614	0.07 % ↑

  • start      = parameter value before refinement
  • value      = refined value from least-squares minimization
  • s.u.       = standard uncertainty (one sigma), from the covariance matrix
  • change     = relative change from start, in %; ↑ = increase, ↓ = decrease

The broad_gauss_sigma_1–broad_gauss_sigma_2 pair is absent because broad_gauss_sigma_2 was not varied. This does not show that the physical ambiguity disappeared; it shows that the fixed parameter is no longer part of the estimated covariance or posterior.

🎲 Exercise 4: Set New Sampling Bounds#

Derive finite bounds from the new local-fit uncertainties and verify them. Why is it better to recalculate the bounds than to copy them from project_1?

Hint:

Iterate over project_2.free_parameters and call set_fit_bounds_from_uncertainty().

Solution:

Free parameters for both structures (🧩 data blocks) and experiments (🔬 data blocks)

	datablock	category	entry	parameter	value	uncertainty	min	max	units
1	lbco	cell		length_a	3.89191	0.00028	3.89079	3.89303	Å
2	sim_lbco	linked_structure	lbco	scale	43.64981	0.37634	42.14446	45.15517
3	sim_lbco	linked_structure	si	scale	0.32115	0.02075	0.23815	0.40415
4	sim_lbco	peak		rise_alpha_1	0.01231	0.00016	0.01167	0.01296	μs/Å
5	sim_lbco	peak		decay_beta_1	0.01148	0.00045	0.00969	0.01327	μs/Å
6	sim_lbco	peak		broad_gauss_sigma_1	-78452.49947	1074.76142	-82751.54514	-74153.45380	μs/Å

Fixing one correlated parameter changes the local covariance matrix and therefore the uncertainty estimates of the remaining parameters. The new bounds should be based on this new local problem.

🎲 Exercise 5: Repeat DREAM Sampling#

Configure DREAM with the same short-chain settings used for project_1, then sample the six-parameter posterior.

Hint:

Use 300 sampling steps, 60 burn-in steps, and random seed 42 so the two runs use comparable settings.

Solution:

⚠️ Switching minimizer type removes these settings:                                                                               
   • max_iterations                                                                                                               

⚠️ Switching minimizer type adds these settings with defaults:                                                                    
   • burn_in_steps=600                                                                                                            
   • initialization_method='latin_hypercube'                                                                                      
   • parallel_workers=0                                                                                                           
   • population_size=4                                                                                                            
   • random_seed=None                                                                                                             
   • sampling_steps=3000                                                                                                          
   • thinning_interval=1                                                                                                          

Current minimizer changed to

bumps (dream)

Standard fitting

📋 Using experiment 🔬 'sim_lbco' for 'single' fitting

🚀 Starting fit process with 'bumps (dream)'...

📈 Bayesian sampling progress:

	step	progress	time (s)	log posterior	phase
1			0.21	-325.24	pre-processing
2	21/360	5.8%	4.29	-269.94	burn-in
3	40/360	11.1%	8.09	-268.12	burn-in
4	60/360	16.7%	12.40	-266.68	burn-in
5	61/360	16.9%	12.59	-266.39	sampling
6	76/360	21.1%	15.64	-266.48	sampling
7	91/360	25.3%	18.70	-266.67	sampling
8	106/360	29.4%	21.78	-266.12	sampling
9	121/360	33.6%	24.84	-266.15	sampling
10	136/360	37.8%	28.25	-266.81	sampling
11	151/360	41.9%	31.31	-266.30	sampling
12	166/360	46.1%	34.43	-265.99	sampling
13	181/360	50.3%	37.51	-266.23	sampling
14	196/360	54.4%	40.61	-265.93	sampling
15	211/360	58.6%	43.74	-265.68	sampling
16	226/360	62.8%	46.87	-266.31	sampling
17	241/360	66.9%	49.98	-265.88	sampling
18	256/360	71.1%	53.06	-266.13	sampling
19	271/360	75.3%	56.48	-265.67	sampling
20	286/360	79.4%	59.59	-265.67	sampling
21	301/360	83.6%	62.73	-265.96	sampling
22	316/360	87.8%	65.80	-266.38	sampling
23	331/360	91.9%	68.89	-265.75	sampling
24	346/360	96.1%	72.01	-266.01	sampling
25	360/360	100.0%	75.11	-265.37	sampling
26			106.85		post-processing

✅ Bayesian sampling complete.

⚠️ Convergence diagnostics indicate the posterior may be poorly mixed.

📊 Exercise 6: Compare the Posterior Results#

Exercise 6.1: Check Convergence and Credible Intervals#

Display the Bayesian result table. Compare r-hat, effective sample size, and the credible intervals with the reference run. Did fixing one parameter automatically guarantee convergence?

Hint:

Use project_2.display.fit.results(), then compare the diagnostics with those displayed for project_1 in the introduction.

Solution:

⚙️ Settings used:

	Name	Value	Description
1	sampling_steps	300	Total sampler iterations per chain.
2	burn_in_steps	60	Sampler iterations discarded as warm-up.
3	thinning_interval	1	Sampler thinning interval.
4	population_size	4	Number of chains or walkers.
5	parallel_workers	0	Worker count; 0 uses all available CPUs.
6	initialization_method	latin_hypercube	Sampler initialization method.
7	random_seed	42	Random seed; None uses a system-derived seed.

📋 Bayesian fit results:

	Metric	Value
1	🧪 Sampler	bumps (dream)
2	❌ Overall status	failed
3	💬 Engine message	DREAM sampling completed
4	⏱️ Fitting time (seconds)	106.85
5	📏 Goodness-of-fit (reduced χ²)	1.24
6	📏 R-factor (Rf, %)	6.89
7	📏 R-factor squared (Rf², %)	5.58
8	📏 Weighted R-factor (wR, %)	4.76
9	📉 Best log-posterior	-263.03
10	📊 Convergence status	failed
11	📊 Max r-hat	1.113
12	📊 Min ess bulk	333.492
13	📊 Draws per chain	300
14	📊 Chains	24

📈 Committed parameters:

	datablock	category	entry	parameter	units	start	value	s.u.	change
1	lbco	cell		length_a	Å	3.8919	3.8919	0.0003	0.00 % ↓
2	sim_lbco	linked_structure	lbco	scale		43.6498	43.6498	0.3491	0.00 % ↓
3	sim_lbco	linked_structure	si	scale		0.3211	0.3211	0.0181	0.00 % ↓
4	sim_lbco	peak		rise_alpha_1	μs/Å	0.0123	0.0123	0.0002	0.00 % ↓
5	sim_lbco	peak		decay_beta_1	μs/Å	0.0115	0.0115	0.0004	0.00 % ↓
6	sim_lbco	peak		broad_gauss_sigma_1	μs/Å	-78452.4995	-78452.4995	954.4100	0.00 % ↓

  • start      = parameter value before sampling
  • value      = estimate written back to the project (best posterior sample)
  • s.u.       = standard uncertainty (one sigma), posterior standard deviation
  • change     = relative change from start, in %; ↑ = increase, ↓ = decrease

📊 Posterior distribution:

	datablock	category	entry	parameter	units	median	95% CI	r-hat	ess bulk
1	lbco	cell		length_a	Å	3.8919	[3.8914, 3.8924]	1.090	359.2
2	sim_lbco	linked_structure	lbco	scale		43.5626	[42.8986, 44.2689]	1.091	369.4
3	sim_lbco	linked_structure	si	scale		0.3203	[0.2838, 0.3554]	1.113	379.3
4	sim_lbco	peak		rise_alpha_1	μs/Å	0.0123	[0.0120, 0.0126]	1.102	333.5
5	sim_lbco	peak		decay_beta_1	μs/Å	0.0115	[0.0108, 0.0124]	1.107	343.1
6	sim_lbco	peak		broad_gauss_sigma_1	μs/Å	-78402.8714	[-80143.2173, -76385.3326]	1.110	352.1

  • median       = 50th percentile of the marginal posterior
  • 95% CI       = 95% credible interval (2.5%-97.5%, asymmetric)
  • r-hat        = Gelman-Rubin diagnostic (good convergence: r-hat <= 1.01)
  • ess bulk     = bulk effective sample size (typically >= 400)

⚠️ r-hat > 1.01: Consider longer sampling, better initialization, or reparameterization.
⚠️ ess bulk < 400: Consider longer sampling or reparameterization.

Reducing the dimension can make sampling easier, but a 300-step chain is still short. Convergence must be judged from the diagnostics, not assumed from the number of free parameters.

Exercise 6.2: Inspect Correlations and Pair Relationships#

Display the posterior correlation matrix and pair plot. Confirm that broad_gauss_sigma_2 is absent, then inspect whether broad_gauss_sigma_1 is correlated with any remaining parameter.

Hint:

Use the same display.fit.correlations(max_parameters=5) and display.posterior.pairs(max_parameters=5) calls as in the introduction. Hover over cells and sample points to identify the parameter pairs and values.

Solution:

The original pair is gone because only sampled parameters appear in posterior charts. Look for any remaining elongated contours rather than concluding that all correlations have disappeared. Fixing one parameter can expose or strengthen relationships among the parameters that remain free.

Exercise 6.3: Compare `broad_gauss_sigma_1`#

Plot the marginal posterior for broad_gauss_sigma_1 from both runs. Compare the medians, shapes, and 95% credible intervals. Why might the second interval be narrower?

Hint:

Call display.posterior.distribution() for each project and pass the corresponding broad_gauss_sigma_1 parameter using param.

Solution:

The conditional posterior can be narrower because broad_gauss_sigma_1 no longer shares uncertainty with broad_gauss_sigma_2. That apparent precision comes from assuming the fixed value is exact. It may underestimate the true uncertainty if broad_gauss_sigma_2 is not independently known.

Exercise 6.4: Compare Posterior Predictions#

Plot the same zoomed posterior-predictive region for both analyses. Does fixing broad_gauss_sigma_2 noticeably change the best curve or uncertainty band?

Hint:

Call display.posterior.predictive() for both projects with the same experiment name and the same x_min and x_max values.

Solution:

Two parameterizations can produce similarly good calculated patterns while assigning different uncertainties to individual parameters. This is why parameter correlations, marginal posteriors, and posterior predictions should be interpreted together.

💾 Exercise 7: Save the Project#

Save the second posterior and its MCMC chain.

Hint:

The project directory was set in Exercise 1, so use project_2.save() to update the existing saved project.

Solution:

Saving project 📦 'main' to 'projects/exercise-bayesian-si-lbco-main'

├── 📄 project.edi
├── 📁 structures/
│   └── 📄 lbco.edi
│   └── 📄 si.edi
├── 📁 experiments/
│   └── 📄 sim_lbco.edi
├── 📁 analysis/
│   ├── 📄 analysis.edi
│   └── 📄 mcmc.h5
└── 📁 reports/
    └── 📄 main.html

Final Remarks#

In this part of the notebook, you learned how to:

prepare a refined EasyDiffraction project for MCMC;
reduce runtime by fixing nuisance parameters while recognizing the uncertainty tradeoff;
use a local fit to obtain starting values and finite bounds;
sample a correlated posterior with DREAM;
distinguish a local covariance correlation from a posterior-sample correlation;
read correlation matrices, pair plots, marginal distributions, and posterior-predictive plots; and
understand how fixing one member of a correlated pair changes the statistical question and the reported uncertainty.

For a scientific analysis, run longer chains, verify convergence, examine sensitivity to bounds and fixed values, and reconsider the diffraction model when systematic residuals remain.

🎁 Bonus#

Congratulations — you’ve now completed Part 2 of the diffraction data analysis exercises for the DMSC Summer School!

If you’d like to keep exploring, the EasyDiffraction library offers many additional tutorials and examples on the official documentation site: 👉 https://docs.easydiffraction.org/lib/latest/tutorials

Besides the Python package, EasyDiffraction also comes with a graphical user interface (GUI) for deterministic diffraction refinement workflows. Bayesian analysis is not yet available in the GUI, so MCMC workflows currently require the Python library.

If you prefer a point-and-click interface over coding, the GUI provides a user-friendly way to perform deterministic refinements. You can download it as a standalone application here: 👉 https://easydiffraction.org

We’d love to hear your feedback on EasyDiffraction — both the library and the GUI! 💬

Bayesian analysis of NPD data

Contents

Bayesian analysis of NPD data#

🛠️ Import Library#

📘 Introduction: MCMC Analysis of LBCO+Si#

📂 Load the Refined Project#

🎯 Choose the Sampled Parameters#

🚀 Run a Local Refinement#

🔗 Understand the Local Correlation Chart#

🎲 Define the Sampling Region#

🎲 Run DREAM Sampling#

📋 Understand the Bayesian Fit Summary#

🔗 Understand Posterior Correlations#

🗺️ Understand the Posterior Pair Plot#

📈 Understand Marginal Posterior Distributions#

📊 Understand the Posterior-Predictive Plot#

💪 Exercise: Fix One Correlated Parameter#

📂 Exercise 1: Create a Fresh Project#

🎯 Exercise 2: Reduce the Free-Parameter Set#

Exercise 2.1: Fix the Background#

Exercise 2.2: Fix broad_gauss_sigma_2#

🚀 Exercise 3: Repeat the Local Refinement#

🎲 Exercise 4: Set New Sampling Bounds#

🎲 Exercise 5: Repeat DREAM Sampling#

📊 Exercise 6: Compare the Posterior Results#

Exercise 6.1: Check Convergence and Credible Intervals#

Exercise 6.2: Inspect Correlations and Pair Relationships#

Exercise 6.3: Compare broad_gauss_sigma_1#

Exercise 6.4: Compare Posterior Predictions#

💾 Exercise 7: Save the Project#

Final Remarks#

🎁 Bonus#

Exercise 2.2: Fix `broad_gauss_sigma_2`#

Exercise 6.3: Compare `broad_gauss_sigma_1`#