Bayesian Analysis of Many-Pole Fits of Hadron Propagators in Lattice QCD

TL;DR

This paper applies Bayesian probability to analyze many-pole fits of hadron propagators in lattice QCD, providing a new method for estimating parameters and their uncertainties with advantages over traditional approaches.

Contribution

It introduces a Bayesian approach to fit analysis in lattice QCD, offering a novel method for parameter estimation and model comparison.

Findings

Bayesian method yields more reliable parameter estimates.

The probability distribution of fit parameters is explicitly calculated.

Model probabilities are quantitatively assessed.

Abstract

We use Bayes' probability theorem to analyze many-pole fits of hadron propagators. An alternative method of estimating values and uncertainties of the fit parameters is offered, which has certain advantages over the conventional methods. The probability distribution of the parameters of a fit is calculated. The relative probability of various models is calculated.

Tables

Table 1. Table 1: Values of χ 2 superscript 𝜒 2 \chi^{2} for fitting 1 and 2 pole data of different noise levels with 1 and 2 pole models.

data	# of propagators	model	${\chi }^2/dof$
		1 pole	0.61
	360	2 pole	0.57
1 pole		1 pole	1.2
	3600	2 pole	1.3
		1 pole	0.87
	360	2 pole	0.99
2 pole		1 pole	2.0
	3600	2 pole	1.2

Table 2. Table 2: Total probabilities ratio R 𝑅 R .

data	# of propagators	R
	360	3
1 pole	3600	30
	360	3
2 pole	3600	0.02

Equations

$$P(T|D,I)P(D|I)=P(D|T,I)P(T|I),$$

$$P(\{c,E\}|D)=\frac{P(\{c,E\})}{P(D)}P(D|\{c,E\}|).$$

$$P(E_n|D)=\int \Pi id{c}_i\Pi i\ne nd{E}_{i}P(\{c,E\}|D).$$

$${\overline{E}}_n=\int d{E}_n{E}_{n}P(E_n|D),$$

$${\sigma }_{E_i}^2=\int d{E}_n(E_n-\overline{E_n}{)}^{2}P(E_n|D),$$

$$\frac{P(T_i|D)}{P(T_j|D)}=\frac{P(T_i)P(D|T_i)}{P(T_j)P(D|T_j)}.$$

$$P(D|T)=\int \{dc\}\{dE\}P(\{c,E\})P(D|\{c,E\}).$$

$$P(D|\{c,E\})=e^{-{\chi }^2(\{c,E\})/2},$$

$$P(E_1|D)=\frac{\int d{c}_1{e}^{-{\chi }^2(c_1,E_1)/2}}{\int d{c}_{1}d{E}_1{e}^{-{\chi }^2(c_1,E_1)/2}}.$$

$$P(D|1pole)=\int d{c}_{1}d{E}_1\frac{e^{-{\chi }^2(c_1,E_1)}}{c_1{E}_1}$$

Full text

1 Introduction

A common procedure in Lattice QCD is to calculate a correlation function in a certain channel and then fit it as a sum of several exponentials [1, 2, 3, 4]. The parameters of the fits are estimated by minimizing $\chi^{2}$. To find the errors, ideally the calculations should be repeated many times, but this is impractical. Usually the jacknife or the bootstrap is performed. Instead here we use Bayes’s theorem to derive the parameters’ probability distribution for given data from the probability of the data for given parameters. Usually one determines the number of poles by comparing $\chi^{2}$ for different models. It is very often ambiguous. Here we use the parameters probability distribution to calculate relative probabilities of the possible models with given data. We perform the Bayesian analysis for some model functions imitating Lattice QCD propagators. Then we apply this method to analyze SU(2) hadron propagators.

2 Bayes’ theorem and its applications

$P(A|B)$ is the conditional probability that proposition $A$ is true, given that proposition $B$ is true. Bayes’ theorem reads:

[TABLE]

where $T$ is the theoretical model to be tested, $D$ is the data, and $I$ is the prior information. $P(T|D,I)$ is the posterior probability of the theoretical model. $P(T|I)$ is the prior probability of the theoretical model. $P(D|I)$ is the prior probability of the data; it will be always absorbed into the normalization constant. $P(D|T,I)$ is the direct probability of the data. For shortness $I$ will be implicit in all formulae henceforth.

We are interested in calculating the posterior probability of a theoretical model and its parameters $\{c,E\}$:

[TABLE]

Given $P(\{c,E\})$, and $P(D|\{c,E\})$ we can [6]:

I. Calculate the posterior probability density $P(E_{n}|D)$ [$P(c_{n}|D)$] for the parameters $E_{n}$ [$c_{n}$] :

[TABLE]

II. Calculate the average values $\overline{E}_{i}$ and the standard deviations $\sigma_{E_{i}}$ (similarly for $c_{i}$)

[TABLE]

provided $P(E_{n}|D)$ is normalized.

III. Compare several models $T_{i}$, (for example one pole, two pole and three pole models). We cannot find the absolute probability of a theory, since we do not have the “complete set” of theories. But we can calculate the relative probabilities of two theories:

[TABLE]

$P(D|T)$ can be obtained from $P(D|\{c,E\})$ by integrating over all parameters of the theory:

[TABLE]

For the prior probability $P(\{c,E\})$ of the parameters of a model in section 4 we make the “least informative” assumption [6] $P(c,E)$ $dc$ $dE$ $\sim$ $dc/c$ $dE/E$. This form is scale invariant. Priors $P(\{c,E\})$ should be normalized.

The direct probability $P(D|\{c,E\})$ of the data $D$ can be calculated relatively easily if the data is Gaussian distributed. We generate “fake” data to be used in the analysis. We use an n-pole model, add noise $e(t)$ to generate a sample of $N$ “propagators” $g_{\alpha}(t)$= $\sum_{i=1}^{n}$ $c_{i}e^{-E_{i}t}+e(t)$ ($\alpha=1,..,N$), calculate the average $G(t)$ and estimate the covariance matrix $C(t,t^{\prime})$ from the data. Here t is the discrete, “lattice” time. We vary the number of “propagators” to control the noise level in the data. Here we use $N=360$ and $N=3600$, which corresponds to a decrease in the noise level by a factor of 3.

The probability distribution of $G(t)$ is [7]

[TABLE]

where $\chi^{2}$ is calculated using the full covariance matrix [7]. The individual $g_{\alpha}(t)$ need not be Gaussian distributed, as long as we average over enough of them so that $G(t)$ are. Gaussian distribution of the “fake” data is ensured.

3 Estimating the Parameters.

3.1 1 pole data.

We generate data $D$ for the one pole model with $c_{1}^{in}=0.15$ and $E_{1}^{in}=0.485$. We use the one pole model to fit the data. Here we assume the prior probability of the data $P(\{c,E\})$ to be constant. Then the posterior probability of the parameters $P(\{c,E\}|D)$ is up to a constant equal to the direct probability of the data given by equation (5). The posterior probability density for $E_{1}$

[TABLE]

It begs for the Monte Carlo integration with the Metropolis algorithm.

We generate a set of points $(c_{1},E_{1})$. Every point is characterized by $\chi^{2}(c_{1},E_{1})$. We sample the vicinity of the minimum of $\chi^{2}(c_{1},E_{1})$ maximum of $exp{(-\chi^{2}(c_{1},E_{1}))}$.

\begin{picture}(0.0,130.0) \put(71.0,30.0){\leavevmode} \put(71.0,5.0){\leavevmode} \end{picture}

We make a scatter plot $c_{1}$ vs.$E_{1}$ (Fig. 2) to visualize this distribution. The density of the points is proportional to the weight $exp{(-\chi^{2}(c_{1},E_{1}))}$.

\begin{picture}(0.0,160.0) \put(71.0,40.0){\leavevmode} \put(71.0,0.0){\leavevmode} \end{picture}

Taking integrals (6) is equivalent to making a histogram with steps big enough to make the distribution smooth (Fig.3).

\begin{picture}(0.0,130.0) \put(71.0,30.0){\leavevmode} \put(71.0,0.0){\leavevmode} \end{picture}

Both equation (2) and jacknife give the same values for the parameters and errors of $E=0.4851(1)$ and $c=0.1499(2)$. The present approach conserves computational time compared to the jacknife. If one uses simulated annealing to fit the data, one covers the same regions in the ($c_{i}$,$E_{i}$) space as needed to calculate probability distributions (2). With the probability distributions one immediately obtains the parameters and errors, whereas with the jacknife the fitting has to be repeated $N$ times.

3.2 2 Pole data

We repeat the analysis performed in section 3.1 for two-pole data when the poles are well separated: $c_{1}^{in}$=0.1, $c_{2}^{in}$=0.1 $E_{1}^{in}$=0.5 $E_{2}^{in}$=0.6. We use the two-pole model to fit the data. The probability distributions for the energies and coefficients are obtained and the parameters are estimated just as in section 3.1.

The only complication is that now we have to deal with the 4-d space of parameters $(c_{1},E_{1},c_{2},E_{2})$. We perform a Monte Carlo integration as described above. The 4-d probability distribution is visualized by projecting it onto two planes $(c_{1},E_{1})$ and $(c_{2},E_{2})$. Each blob in Fig. 4 is a projection of the 4-d distribution on a 2-d plane. Each blob represents the probability distribution for one pair of $(c_{i},E_{i})$ after the second pair has been integrated out.

\begin{picture}(0.0,150.0) \put(71.0,40.0){\leavevmode} \put(71.0,0.0){\leavevmode} \end{picture}

4 Model selection

Unless one has some prior knowledge, determining the number of poles present in the data based on comparing $\chi^{2}$ for models with different numbers of poles is very often ambiguous.

Table 1 contains the values of $\chi^{2}$ obtained when 1 and 2 pole data are fit with 1 and 2 pole models. We use the one pole data from section 3.1 and generate 2 pole data with the poles close to each other: $c_{1}^{in}=0.05$ $c_{2}^{in}=0.1$ $E_{1}^{in}=0.49$ $E_{2}^{in}=0.5$. With one exception the corresponding values of $\chi^{2}$ are close and do not allow for a definite answer.

We need to calculate the total probabilities ratio (3) to determine the number of poles in the data. We assume $P(2pole)=P(1pole)$, i.e. a priori these two models are equally probable. To calculate $P(D|1pole)$ we substitute (5) into equation (4) (and similarly for $P(D|2pole)$)

[TABLE]

The integration is tricky since we are dealing with a function that varies rapidly in the multidimensional space. We use scatter plots to determine the areas of integration. Table 2 contains integration results for the total probability ratio $R=\frac{P(1pole|D)}{P(2pole|D)}$. With a sufficiently low noise level the total probabilities ratio picks the correct model.

If we estimate the parameters of the 2 pole data with 3600 propagators using equation (2), we get the input parameters back within the error bars. The jacknife here gives unreasonably large errors because $\chi^{2}$ has several minima in the ($c_{1},c_{2},E_{1},E_{2}$) space. It can be seen from the graph of $\chi^{2}$ vs. $n$ (Fig.7), and they can be also clearly identified on the scatter plot (Fig.8). Two minima ( 1 and 3) have the same $\chi^{2}$. When we perform the jacknife, the blind fitting finds either minimum 1 or 3, which results in the unreasonably large error estimates.

\begin{picture}(0.0,165.0) \put(71.0,50.0){\leavevmode} \put(71.0,20.0){\leavevmode} \end{picture}

\begin{picture}(0.0,135.0) \put(71.0,40.0){\leavevmode} \put(71.0,0.0){\leavevmode} \end{picture}

5 Analysis of SU(2) data

Here we analyze hadron propagators in the pseudoscalar channel. The detailed description of these data is given in [4, 5]. The coupling constant $\beta=2.5$, lattice spacing $a=0.09\pm 0.012$ fm, the lattice size is $12^{3}\times 24$. The propagators analyzed here were calculated with $\kappa=0.146$ for 360 configurations. The point source is at $t=5$. The fit is performed in the time range 6-20.

We repeat the analysis for 60 and 360 configurations to study the data with different noise levels. For 60 configurations the $\chi^{2}/dof$ is 1.0 for the 3 pole fit and 1.17 for the 4 pole fit, and the total probabilities ratio is $\frac{P(3pole|D)}{P(4pole|D)}\sim 1$. For 360 configurations the $\chi^{2}/dof$ is 0.52 for the 3 pole model fit and 0.69 for the 4 pole model fit, and the total probabilities ratio for 360 configurations is $\frac{P(3pole|D)}{P(4pole|D)}\sim 10$. Here again, for low noise data we are able to choose between two models based on a qualitative estimate given by the total probabilities ratio.

6 Conclusion

A new method has been introduced that can be used to analyze many-pole fits of hadron propagators in Lattice QCD. It has been used to estimate the many-pole model parameters and their uncertainties. It works in the presence of multiple minima when the jacknife (at least in its simpleminded form) fails. It cuts the computational time. The new method has been used to calculate relative probabilities for different models, which can be crucial in making the optimal choice of a model. The scatter plots, which have been introduced as an auxiliary tool for the multidimensional integration, can be used as an independent tool for the many pole fit analysis.

I wish to thank Greg Kilcup for his suggestion to use Bayesian analysis.