Let u i (x), with i ranging from 1 to m, be the reduced potential energy functions 13 of m states. We aim to compute the free energy differences among these states by sampling configurations from their Boltzmann distributions. For the i-th state, we use {x ik , k = 1, . . . , n i } to represent the n i configurations sampled from its Boltzmann distribution and we assume that these configurations are uncorrelated. We use F = (F 1 , . . . , F m ) to denote the free energies of the m states. In BayesMBAR, we introduce y ik as the index of the state from which configuration x ik is sampled. Therefore, y ik is equal to i. Although the indices of states for sampled configurations are determined during sampling, they are treated as random variables in BayesMBAR. Specifically, y ik is viewed as a sample from a categorical distribution with parameters π = (n 1 /n, . . . , n m /n), meaning that the probability of sampling a configuration from the i-th state is p(y = i) = π i = n i /n, where n = m i=1 n i is the total number of configurations. The concatenation of state indices and configurations, denoted as (y, x), is viewed as samples from the conditional distribution p(y, x|F ), defined as