{'fraction_non_alphanumeric': 0.0396520253213717, 'fraction_numerical': 0.013670433803254495, 'mean_word_length': 4.047040376323011, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hypothesis testing for Markov chain Monte Carlo

Testing between hypotheses, when independent sampling is possible, is a well developed subject. In this paper, we propose hypothesis tests that are applicable when the samples are obtained using Markov chain Monte Carlo. These tests are useful when one is interested in deciding whether the expected value of a certain quantity is above or below a given threshold. We show non-asymptotic error bounds and bounds on the expected number of samples for three types of tests, a fixed sample size test, a sequential test with indifference region, and a sequential test without indifference region. Our tests can lead to significant savings in sample size. We illustrate our results on an example of Bayesian parameter inference involving an ODE model of a biochemical pathway.


Introduction
The goal of Markov chain Monte Carlo estimation is to calculate an expected value with respect to a probability distribution from which sampling directly is impossible or impractical.
It is often only of interest whether the expected value is above or below a certain threshold (such as whether E π f > r, where π is a probability distribution and f is the function of interest). This problem can be posed as a decision between two hypotheses, and one must make a decision with a bounded error probability. This has applications in the verification of stochastic systems (here one aims to verify whether a model satisfies a property with at least a given probability), arising in software testing, robotics and systems biology (Legay et al, 2010).
In the case of independent samples, optimal fixed length and sequential hypothesis tests are available for deciding between two hypotheses (see Wald (1945), Lai (1973), Lai (1988), Lehmann and Romano (2005)). However, such tests are not available in the case of samples obtained by MCMC simulation.
Our main contribution in this paper is the introduction of one fixed length test and two sequential tests (one with indifference region, and one without indifference region) that allow us to decide whether the expected value of a quantity exceeds a certain threshold. We prove non-asymptotic bounds for the probability of error (choosing the incorrect hypothesis) and the expected running times of these tests. The advantage of our approach is that the sample size needed to make the decision between the two hypothesis can be much shorter than the one needed to precisely estimate the expected value.
The paper is organized as follows. In Section 2, we briefly review the necessary preliminaries from the theory of Markov chains. In Section 3, we state the three hypothesis tests, and their theoretical properties. Finally, in Section 4, we evaluate these tests on an ODE model of a biochemical pathway, whose parameter posterior is explored using an MCMC chain.

Markov chain preliminaries
In this section, we review some basic definitions about Markov chains.

Spectral gap of general state Markov chains
Firstly, we state the definition of the spectral gap of reversible Markov chains following Roberts and Rosenthal (2004). The spectral gap is a measure of how fast the chain is mixing, and will be needed to state our error bounds for the hypothesis tests in this paper. The main reason we are using the spectral gap is that it is the basis of existing non-asymptotic bounds for MCMC empirical averages (see León and Perron (2004)) that are known to be sharp in some cases.
( 2.1) Define L 2 (π) as the Hilbert space of complex valued measurable functions that are square integrable with respect to π, endowed with the inner product f, g π = f g * dπ. P can be then viewed as a linear operator on L 2 (π), denoted by P , defined as and reversibility is equivalent to the self-adjointness of P . The operator P acts on measures to the left, i.e. for every measurable subset A of Ω, For a Markov chain with stationary distribution π, we define the spectrum of the chain as S 2 := {λ ∈ \ 0 :(λI − P ) −1 does not exist as a bounded linear operator on L 2 (π)}.
For reversible chains, S 2 lies on the real line. The following is the main definition of this section.
Definition 2.1. The spectral gap for reversible chains is defined as if eigenvalue 1 has multiplicity 1, and 0 otherwise.
Similary, we define the absolute spectral gap as if eigenvalue 1 has multiplicity 1, and 0 otherwise.
It follows from the definition that γ ≥ γ * . They are often equal in practice, and it is simple to modify the MCMC steps such that γ = γ * (by considering the so-called "lazy" version of the chain). For reversible chains, it is known that γ * > 0 implies geometric ergodicity, and the CLT for any function f ∈ L 2 (π). Moreover, for any initial distribution ν, and any k, so the absolute spectral gap is related to the speed of convergence to equilibrium (in total variational distance).
In the case of non-reversible chains, Paulin (2013) defines the pseudo-spectral gap, and shows that it has similar properties as the spectral gap has for reversible chains.

Estimating the spectral gap
In practice, the spectral gap is often not known, and has to be estimated from the output of the chain. In the Appendix, we propose an estimator based on the following simple fact.
From the output of the chain we can estimate ρ η (f ), which in turn leads to estimates for γ * . For further details on the exact procedure, as well as some numerical results, we refer the reader to the Appendix. One important point to note is that the proposed procedure is quite fast, and does not affect the overall running time significantly.

Hypothesis tests and error bounds
In this section we present bounds on the error of the hypothesis tests on MCMC estimates. First we review a Hoeffding-type inequality for reversible Markov chains, and then introduce several hypothesis tests based on it.
An important assumption of the concentration inequalities is that the Markov chain is stationary. To ensure this, we are going to discard the first t 0 terms of the Markov chain (t 0 is commonly called the burn-in time (Gilks et al, 1996)), and only take into account the terms f (X t0+1 ), . . . , f (X n ). Similarly to Gyori and Paulin (2014), we will set t 0 ≥ 30/γ. In what follows, we will not show t 0 explicitly, and assume that sufficient burn-in steps were performed such that f (X 1 ),. . .,f (X n ) is approximately stationary.

Concentration of the MCMC estimate
In this section, we review a Hoeffding-type inequality for reversible Markov chains (this will be used to bound the errors of our hypothesis tests). This result was proven for Markov chains on finite state spaces in León and Perron (2004), and extended to general state spaces in Miasojedow (2014).
Theorem 3.1 (Hoeffding inequality for reversible chains). Let X 1 , . . . , X n be a stationary, reversible Markov chain with spectral gap γ, and unique stationary distribution π. Let f ∈ L 2 (π) such that 0 ≤ f (x) ≤ 1 for every x ∈ Ω. Let S n := n i=1 f (X i ), then for any t ≥ 0, This form of the result follows from equation (3) of León and Perron (2004) by rescaling. A similar inequality can be shown to hold for non-reversible chains with the spectral gap being replaced by the pseudo spectral gap, see Paulin (2013) for more details.
In the case of non-stationary chains, one can show that essentially the same result holds (see Propositions 3.3 and 3.4 of Paulin (2013)).

Hypothesis tests with indifference region
Suppose that X 1 , X 2 , . . . is a reversible Markov chain taking values in a Polish state space Ω, with unique stationary distribution π, and f : Ω → [0, 1] is a bounded function (with a simple scaling argument, our results extend to the case when f : Ω → [α, β]).
Here (r − δ, r + δ) is an indifference region in which choosing either hypothesis is acceptable.
We first discuss two tests to choose between these hypotheses. The first one is a fixed sample size test, while the second one is a sequential test.

Fixed length hypothesis test
Suppose that we have a sample of length n consisting of the values f (X 1 ), . . ., f (X n ). The fixed length hypothesis test is shown in Algorithm 1.
The next proposition bounds the error probability of the test in Algorithm 1, where an error constitutes accepting hypothesis H 1 when in fact hypothesis H 0 holds, and vice-versa.

Algorithm 1 Fixed length hypothesis test
Input: Threshold r, number of samples N Output: Choice of H0 or H1.
Applying the Hoeffding inequality (Theorem 3.1), we get The same holds under the opposite hypothesis.
This implies, in particular, that if we want the error to be smaller than ǫ, then samples are sufficient.

Sequential hypothesis test
A sequential test is shown in Algorithm 2. The main idea behind this test is to sequentially monitor the empirical sum, and stop if the samples collected so far are sufficient to decide between the hypotheses. The following proposition bounds the error probability of the test in Algorithm 2 with a particular choice of M .

Proposition 3.3 (Error bound for sequential hypothesis test).
For the sequential hypothesis test, the probability of an error is bounded by

Algorithm 2 Sequential hypothesis test with indifference region
Input: Threshold r, stopping condition M Output: Choice of H0 or H1.
Set S0 := 0 and n := 1 loop Set n := n + 1 and continue end if end loop It is easy to see that the probability of choosing H 0 is bounded by the sum ∞ n=1 P(S n > nr + M ).
( 3.6) Since n ≤ nr + M for n ≤ M/(1 − r), and S n ≤ n, the first ⌊M/(1 − r)⌋ terms of (3.6) will be 0. By Theorem 3.1, for any n ≥ 1, we have The same holds under the opposite hypothesis.
From (3.5), it is easy to see that the probability of an error with the sequential test will be smaller than ǫ if M is chosen as . (3.7)

Sequential hypothesis test without indifference region
Suppose that, unlike previously, we do not want to use an indifference region. This case arises when our objective is to do a test between the following two hypotheses given r ∈ (0, 1): In this case, we cannot use a fixed length hypothesis test. We propose the following modified version of the sequential hypothesis test of the previous section. Let ǫ > 0 be the specified error probability of the test, and let

Algorithm 3 Sequential hypothesis test without indifference region
Input: Threshold r, error bound ǫ Output: Choice of H0 or H1.
Set S0 := 0 and n := 1 loop Set n := n + 1 and continue end if end loop The following proposition bounds the error probability of the test in Algorithm 3.
Proposition 3.4 (Error bound for sequential hypothesis test without indifference region). For the test, the error probability is bounded by ǫ.
Proof of Proposition 3.4. Suppose H 1 holds, implying that E π f < r. Using Hoeffding's inequality for reversible Markov chains (Theorem 3.1), the probability of choosing H 0 is bounded by The same holds under the opposite hypothesis.

Expected stopping times
In this section we analyze the (expected) number of samples taken in each test, given ǫ, r and δ.
The fixed length hypothesis test will always take log(1/ǫ)/ (γδ 2 ) steps to decide between the hypotheses in (3.2) with error at most ǫ. We now show the expected stopping time of the sequential test with indifference region in Proposition 3.5.
Proposition 3. 5. For the sequential test with indifference region, with M chosen as in (3.7), the expected stopping time satisfies under both hypotheses.
Proof. Under hypothesis H 1 , the probability of the event that T ≥ M/δ + t is less than equal to the event that S n > nr − M at step n : Using the Hoeffding inequatiy we see that for any t > 0 To arrive at the expected stopping time, we use the fact that for any real valued random variable, where the second step comes from an upper bound for the exponential integral.
We now substitute M as in (3.7) to get The same holds under the opposite hypothesis.
By comparing (3.8) with (3.4), we can see that the sequential test typically requires less than half of the samples of the fixed length test with the same errors.
The expected time of the test grows essentially linearly in M . Under hypothesis H 0 , let us denote E π f by E 0 f . Then one can also show that E 0 T ≤ M r−E0f +O( M/δ), which can be much smaller than M/δ if E 0 f is much smaller than r − δ. In such situations, the test will generally use less samples than a test with fixed length would require, resulting in shorter running time.
Simple arguments show that the test will stop in a finite amount of time almost surely, even if none of the two hypotheses is satisfied, that is, if E π f ∈ (r − δ, r + δ). In practice, however, one may need to stop the run after a certain number of steps. In particular, by (3.9), the probability that the chain runs for more than 2M/δ steps is less than
Proof. We first show that and by applying the Hoeffding inequality, this can be further bounded by for any n ∈ AE, with x + denoting the positive part of x ∈ Ê.
By the definition of N , it is easy to see that for n ≥ N , we have g(n, ǫ) ≤ n|E π f − r|/2. Using the fact that E(T ) = ∞ n=1 P(T ≥ n), we have where we have used the definition of N and the fact that for Our bound on the expected time of the test without indifference region is essentially of the same order of magnitude as for the sequential test with indifference region, except that δ is replaced by |E π f − r|.

Case study
Here we present a case study to evaluate various aspects of our hypothesis tests empirically. The case study is from the domain of systems biology. Dynamical system models based on ordinary differential equations (ODEs) are often used to describe the concentration of molecular species such as proteins inside the cell (Klipp et al, 2005). However, the rate constants associated with biological processes are often unknown and not directly measurable. This appears as a set of unknown parameters in the ODE model whose value can only be inferred given limited and noisy experimental data. In a Bayesian inference setting, sampling directly from the posterior distribution of parameters will not be possible, and a Markov chain Monte Carlo method can be used to collect samples from the posterior (Lawrence et al, 2009).
We assume a dynamical system model of the forṁ Here x ∈ Ê dx is a vector of state variables, y ∈ Ê dy is a vector of observables, θ ∈ Ê d θ is a vector of parameters, and ω is a random variable representing measurement noise.
The goal is to construct the posterior distribution of θ given its prior p 0 (θ), and a set of observations Y . Here Y consists of observations of the form Y i,j , where 1 ≤ i ≤ d y denotes the ith observable of the system at time point t j , 1 ≤ j ≤ T . We denote the likelihood of the set of observations with respect to θ as p(Y |θ). The posterior, denoted π(θ|Y ), is then expressed in the standard way as In this case study we use a uniform prior over a bounded set of parameters, and Gaussian likelihood (corresponding to uncorrelated multivariate Gaussian ω) defined as Here y i (t j )| θ denotes the ith observable of the model at time t j , when simulated with parameters θ, and σ i,j is the standard deviation associated with data point Y i,j .
The analysis of such a system involves estimating the expected value of a function f of θ with respect to π(θ|Y ), denoted E π f . Closed form solutions to this problem will, in general, not be available, and sampling independently from the posterior π(θ|Y ) will not be possible. We therefore resort to using a Metropolis-Hastings chain to collect a sequence of parameters θ 0 , θ 1 , . . ., θ n , and use the approximation Each state of the chain is obtained by a proposal and an acceptance step. We use a symmetric Gaussian proposal of the form q(θ i → θ ′ ) = N (θ i , Σ MH ), a d θdimensional multivariate Gaussian with mean identical to the current parameter vector, and covariance matrix Σ MH . Here Σ MH is diagonal with entries σ 2 MH,1 , . . ., σ 2 MH,d θ , representing variances along each dimension independently. The proposed parameter θ ′ is accepted with probability α(θ i → θ ′ ), determined by the posterior as Note that the proposal does not appear in the acceptance ratio due to symmetry.
We apply our method to analyze the dynamics of the JAK-STAT biochemical pathway (for more details, see Swameye et al (2003)). Our goal will be to decide about an important property of the system, namely, whether the concentration of nuclear STAT protein reaches the threshold of 1. The variables in the model and the 4 model parameters (whose values are not known) cannot be directly measured. However, experimental data for two indirect quantities (total phosphorylated STAT, and total STAT in cytoplasm) has been published in Swameye et al (2003). We let when simulating with θ, 1 if nuclear STAT reaches 1 when simulating with θ The equations governing the system dynamics and the parameters used to run the MCMC chain are given in the Appendix.
We ran m = 1000 independent instances of the MCMC chain for a total of 2 · 10 6 steps each (with an additional t 0 = 5 · 10 4 burn-in steps). The method described in the Appendix was used to estimate the value of the spectral gap, independently for each chain. We used the output of the m independent chains as a basis for constructing empirical results in Figure 1(a-f).
To get a reliable estimate of the true expected value, E ≈ E π f , we took the overall average of the estimates from all m chains, and treated the obtained value E = 0.8875 as the reference for E π f .
We first examined the empirical error rate of the fixed sample size hypothesis test. For a fixed sample size n, we define the empirical error rate E n as the ratio of chains choosing H 0 if H 1 holds (or the ratio choosing H 1 if H 0 holds). If neither H 0 nor H 1 holds (when r − δ < E π f < r + δ), then E n := 0. We set r = E − δ and calculated E n for a range of sample sizes up to n = 10 6 . For the same set of sample sizes, we calculated the average error rate bound derived from equation (3.4) as ǫ n = exp(−nγδ) (the average is used as the estimate of γ is different across independent runs). Figure 1(a) shows E n and ǫ n as a function of n for different values of δ. Importantly, Figure 1(a) demonstrates that the theoretical bounds on the error rate are reliable in practice, since the empirical error rate is consistently below this upper bound (E n ≤ ǫ n for all examined n, δ).
We next look at results for sequential hypothesis testing with indifference region. We refer to the number of samples collected in the Markov chain before a decision is made as the stopping time (see also Section 3.4). Figures 1(b-c) show the mean empirical stopping times for a range of r values for different values of δ (b), and different values of ǫ (c). For values of r close to E π f , some chains did not stop within 2 · 10 6 samples, and the corresponding mean values are therefore not determined. In both plots, the average of the fixed sample sizes needed for each chain is also shown. In Figure 1 (d), the empirical cumulative distribution of stopping times is shown for the hypothesis test for a set of r values in (0, 1). Here the value of δ = 0.05 and ǫ = 0.01 is fixed. As a reference, we also show the empirical distribution of sample sizes needed to perform the fixed length test, that differ across chains due to the different estimates of the spectral gap. The plot shows that for values of r distant from the true average, sequential sampling consistently terminates with small variability at low sample sizes. When r is close to the true average, the stopping times show higher variability.
Finally, we used the sequential hypothesis test without indifference region. Figure 1(e) shows empirical stopping times for different values of ǫ and a range of r values. Figure 1 (f) shows the empirical cumulative distribution of stopping times for the hypothesis test for a set of r values. Here the value of ǫ = 0.01 is fixed.
We evaluated the empirical error rate in the sequential hypothesis test with indifference region, and found that out of all 1000 runs, no error was made under any examined choice of r, ǫ, δ. This, shows that the specified error bound (3.5) was indeed met. It also suggests that the bound (3.5) might not be sharp and M could be chosen even smaller than described by (3.7), resulting in earlier stopping. Similarly, for the test without indifference region, no decision error was made for any of the runs under examined choices of r and ǫ.

Conclusion
In this paper we proposed hypothesis tests on MCMC estimates. These tests are useful in cases where one is not interested in the exact value of the estimate, but only whether it is below or above a certain threshold. We have stated three different hypothesis tests, a fixed length test, a sequential test with indifference region, and a sequential test without indifference region. Our main theoretical contribution is rigorous error bounds for these test, and on their expected running times. We have illustrated their usage on a case study from the domain of systems biology, using an ODE model of the JAK-STAT biochemical pathway. In our simulations, the sequential tests have performed well, and their running time can be much shorter than the fixed length test, especially when the true expected value is far away from the threshold.
There are some theoretical and practical questions that remain open for further exploration. We have chosen the details of Algorithms 2 and 3, in particular, the function g(n, ǫ), in this specific way in order to be able to show error bounds using concentration inequalities. However, we do not claim that these choices are optimal, and our error estimates could possibly be sharpened using different techniques (in the independent case, there is a well developed literature on how to do sequential tests with or without indifference regions in an optimal way, see Lai (1973), Lai (1988)).

Estimating the spectral gap
In this section, we propose a procedure to estimate the spectral gap of Markov chains. The procedure is motivated by the following lemma. In the statement of the lemma, we are going to use the scalar product where π is the stationary distribution of a finite state Markov chain.
Lemma 5.1. Suppose that X 1 , X 2 , . . . is a finite state Markov chain with spectral gap γ, and transition kernel P . Suppose that f : Ω → Ê satisfies that E π (f ) = 0, and that f is not orthogonal to the eigenspace corresponding to the eigenvector of P of the second largest absolute value (with respect to the scalar product ·, · π ). Let ρ η (f ) : Remark 5. 2. A similar result can be shown for chains with general state spaces, but for notational simplicity we only consider finite state spaces. In practice, the condition of non-orthogonality is almost always satisfied.
Proof. We can write the eigenvalues of the operator P as 1 = λ 1 ≥ λ 2 ≥ . . . ≥ λ n , with n = |Ω|, and corresponding right eigenvectors f 1 , . . . , f n , which are orthonormal with respect to the scalar product ·, · π . It is clear that f 1 = 1, and using the assumption E π (f ) = 0, f can be decomposed as for some constants {c i } n i=1 . Using this decomposition, we have and using the condition of non-orthogonality, we have In practice, we cannot choose η to be infinity, in particular, there are some issues related to the variance of the estimator of ρ η (f ) that need to be taken in to account. Given burn-in time t 0 , and f (X t0+1 ), . . . , f (X t0+n ), we define the empirical mean as and use the estimator The typical range of dependence of the elements of the Markov chain is of the order of 1/γ * , so we expect the standard deviation ofρ η (f )/Var π (f ) to be of the order of 1 √ nγ * . Now if we want to use the estimator the standard deviation of ρ η (f )/Var π (f ) needs to be much smaller than (1−γ * ) η . Solving the equation for which the standard deviation of ρ η (f )/Var π (f ) is much smaller than (1 − γ * ) η . A slight inconvenience of this choice is that it depends on the unknown parameter γ * . A way to overcome this issue is by computing the value of η iteratively. Based on Lemma 5.1 and the above observation, we propose the following procedure.
The motivation for choosing f 1 , . . . , f m in this way is that we use all the available information about the Markov chain, and thus the estimator is expected to be more accurate than if we would only use a subset of this information. The motivation for Step 5 is to ensure that we have sufficient initial data for the estimate.
An illustration of this iterative procedure based on our case study (see Section 4) is shown in Figure 2. Here we use the components of the Markov chain's state (the value of the model parameters k 1 to k 4 ) for the estimation.  A histogram of estimated spectral gap values is shown in Figure 3 for the 1000 independent chains used in our case study.

Simulation details
Here we provide additional details on the JAK-STAT pathway model case study.
The species in the model are listed in Table 1.
The ODE equations governing the model dynamics are as follows. Here y 1 represents total phosphorylated STAT and y 2 represents total STAT in cytoplasm. Observations are available at 18 discrete time points up to 60 minutes. The data points, as well as standard deviations for each point are available in Swameye et al (2003). The parameter vector of the model is θ = (k 1 , k 2 , k 3 , k 4 ). We assume a uniform prior over a hypercube defined by a bounded interval for each parameter. The parameter ranges and the covariance matrix diagonal entries (σ MH ) used to define the MCMC proposal distribution are provided in Table 2. [0, 5] 0.02 Table 2 Parameter ranges and entries in the proposal covariance matrix.