Documentation:Monte Carlo Equilibration

Monte Carlo equilibration

Theory

We have a timeseries of N measurements obtained from a Monte Carlo simulation, i.e. $y_0,y_1,\cdots,y_{N-1}$.

Suppose $\bar{y}_i = \beta_0 + \beta_1 x_i$ (s.t. $i = 0, 1, \cdots, N-1$) is the least-squares best fitted line, we attempt to minimize $S = \sum_i (y_i - \bar{y}_i)^2$ w.r.t. $\beta_0$ and $\beta_1$.

$\frac{\partial S}{\partial \beta_0 } = 0$ , $\frac{\partial S}{\partial \beta_1 } = 0$ :

$\Rightarrow \left( \begin{array}{cc} N & \sum_i x_i \\ \sum_i x_i & \sum_i x_i^2 \end{array} \right) \left( \begin{array}{c} \beta_0 \\ \beta_1 \end{array} \right) = \left( \begin{array}{c} \sum_i y_i \\ \sum_i x_i y_i \end{array} \right)$

$\Rightarrow \left( \begin{array}{c} \beta_0 \\ \beta_1 \end{array} \right) = \frac{1}{N\sum_i x_i^2 - (\sum_i x_i)^2} \left( \begin{array}{cc} \sum_i x_i^2 & -\sum_i x_i \\ -\sum_i x_i & N \end{array} \right) \left( \begin{array}{c} \sum_i y_i \\ \sum_i x_i y_i \end{array} \right)$

$\Rightarrow \beta_1 = \frac{N \sum_i x_i y_i - \sum_i x_i \sum_i y_i}{N\sum_i x_i^2 - (\sum_i x_i)^2}$

Slope of best-fitted line

$\Rightarrow \beta_1 = \frac{\sum_i (x_i - \bar{x}_i)( y_i - \bar{y}_i) }{\sum_i (x_i - \bar{x}_i)^2} \,\,\,\,\, \left( = \frac{s_{xy}}{s_{xx}} \right)$

Error in slope of best-fitted line

$\mathrm{Var}({\beta_1}) = \mathrm{Var} \left( \frac{\sum_i (x_i - \bar{x}_i)( y_i - \bar{y}_i) }{\sum_i (x_i - \bar{x}_i)^2} \right) = \frac{\sum_i (x_i - \bar{x}_i)^2 \mathrm{Var}(y_i) }{(\sum_i (x_i - \bar{x}_i)^2)^2}$

Denoting $\mathrm{Var}(y_i) = \sigma_{y}^2$, we have:

$\Rightarrow \mathrm{Var}({\beta_1}) = \sigma_{\beta_1}^2 = \frac{\sigma_y^2}{\sum_i (x_i - \bar{x}_i)^2}$

$\Rightarrow \mathrm{Var}({\beta_1}) = \sigma_{\beta_1}^2 = \frac{ 12 \, \sigma_y^2 }{ N(N^2 - 1) }$

Hypothesis testing

$\begin{array}{l} H_0 : \beta_1 = 0 \\ H_1 : \beta_1 \ne 0 \end{array}$

Using the standard z-test, we reject $H_0$ at confidence interval $\alpha$ if

$\left| \frac{\beta_1 - 0 }{ \sigma_{\beta_1} } \right| > z_{\frac{1-\alpha}{2}}$

Implementation in Python

Synposis

pyalps.checkSteadyState(outfile, observables, confidenceInterval=0.63, simplified=False, includeLog=False)

Arguments

argument	default	type	remark
outfile	-	Python str	ALPS hdf5 output file(name)
observables	-	(list of ) Python str	(list of) measurement observable(s)
confidenceInterval	0.63	Python float	confidence interval in which steady state has not been reached
simplified	False	Python bool	shall we combine the checks of all observables as 1 final boolean answer?
includeLog	False	Python bool	shall we print the detailed log?