Monte Carlo method is a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. They are often used in physical and mathematical problems and are most useful when it is difficult or impossible to use other mathematical methods. Basically, many statisticians have been increasingly drawn to Monte Carlo method in three distinct problem classes: optimization, numerical integration, and generating draws from a probability distribution. In this paper, we will introduce the Monte Carlo method for calculating regression coefficients in Generalized Linear Model (GLM), especially for Logistic Regression. Our main methods are Metropolis Hastings (MH) Algorithms and Stochastic Approximation in Monte Carlo Computation (SAMC). For comparison, we also get results automatically using MLE method in R software. Then we apply some Monte Carlo algorithms to a real example study and compare the efficiency of each method.

Keywords: Monte Carlo method, GLM, MH-algorithms, SAMC

Markov chain Monte Carlo (MCMC) methods are a class of algorithms for sampling from a probability distribution based on constructing a Markov chain that has the desired distribution as its equilibrium distribution. The Metropolis-Hastings (MH) algorithm is one of the best known of these methods, which was developed by Metropolis et al.¹ and subsequently generalized by Hastings.² Firstly, MH algorithm has been used in physics and was little known to statisticians until Muller³ and Tierney⁴ provided important applications using this algorithm. Then Chib & Greenberg⁵ gave some recent examples including relevant Markov chain theories which made more applications appeared in the recently literature. In the next section, we will provide a brief introduction to the MH-Algorithms and its work principles. We will also discuss the three different MH-Algorithms based on different sampling methods.

As for Stochastic Approximation in Monte Carlo Computation (SAMC), it was first introduced by Liang et al.⁶ As we know, the Metropolis-Hastings (MH) algorithm and the Gibbs sampler Geman⁷ are prone to get trapped into local energy minima in simulations from a system for which the energy landscape is rugged. In terms of physics, the negative of the logarithmic density/mass function is called the energy function of the system. To overcome this problem, many advanced Monte Carlo algorithms have been proposed, such as parallel tempering Geyer,⁸ Hukushima & Nemoto,⁹ simulated tempering Marinari & Parisi,¹⁰ Geyer & Thompson,¹¹ evolutionary Monte Carlo Liang & Wong,¹² dynamic weighting Wong & Liang,¹³ multicanonical sampling Berg & Neuhaus,¹⁴ 1/k-ensemble sampling Hesselbo & Stinchcombe,¹⁵ the Wang-Landau algorithm Wang & Landau,¹⁶ equi-energy sampler Mitsutake et al.¹⁷ stochastic approximation Monte Carlo Liang et al.⁹ among others. In this paper, we adopt the SAMC algorithm to our problem, which has been proofed to be a sophisticated algorithm in both theory and applications. The basic idea of SAMC stems from the Wang-Landau algorithm and extends it to the continuum systems and it is designed to simulate data with a complex model structure.

The rest of the paper is organized as follows. In section 2, we will introduce three different sampling methods based on Metropolis Hastings Algorithms and SAMC algorithm for parameter estimations in the logistic regression. In section 3, we will give a simple simulation on the all of the methods and compare the results with the MLE method provided in R package. In section 4, we apply the proposed method to a Adult Intelligence data. In section 5, we conclude the paper with a brief discussion.

Metropolis hastings algorithms

The working principle of Markov Chain Monte Carlo methods is quite straight forward to describe. Given a target density, we build a Markov kernel  with stationary distribution  and then generate a Markov chain ($X^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}$ ) using this kernel so that the limiting distribution of ($X^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}$ ) is and integrals can be approximated according to the Ergodic Theorem. The difficulty should thus be in constructing a kernel  that is associated with an arbitrary density. But, quite miraculously, there exist methods for deriving such kernels that are universal in that they are theoretically valid for any density $f$ .

The Metropolis Hastings algorithm is an example of those methods. Given the target density, it is associated with a working conditional density $q\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$ that, in practice, is easy to simulate. In addition,  can be almost arbitrary in that the only theoretical requirements are that the ratio $f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo><mo>/</mo>}q\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$ is known up to a constant independent of  and that  has enough dispersion to lead to an exploration of the entire support of $p_{{\theta }_t}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}\propto {\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^m}\frac{\psi \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}}{e^{{\theta }_{ti}}}I\mathrm{<mo\; stretchy="false">(</mo>}x\in E_i\mathrm{<mo\; stretchy="false">)</mo>}$ . For every given , we can then construct a Metropolis Hastings kernel such that  is its stationary distribution. Here is the general step of the Metropolis Hastings algorithm.

Algorithm (Metropolis Hastings)

1. Initialize $x^{\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}}\sim q\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$ .
2. Given current value of $x^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}$ , generate a candidate $y\sim q\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}{x}^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}$ .
3. Calculate acceptance probability $\alpha \mathrm{<mo\; stretchy="false">(</mo>}{x}^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}y\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}min\mathrm{<mo>\{</mo>}\frac{f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}q\mathrm{<mo\; stretchy="false">(</mo>}{x}^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>|</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}}{f\mathrm{<mo\; stretchy="false">(</mo>}{x}^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}q\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}{x}^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,1<mo>\}</mo>}$ .
4. Take $x^{\mathrm{<mo\; stretchy="false">(</mo>}t+\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}y$ with probability $\alpha \mathrm{<mo\; stretchy="false">(</mo>}{x}^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}y\mathrm{<mo\; stretchy="false">)</mo>}$ and $x^{\mathrm{<mo\; stretchy="false">(</mo>}t+\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{x}^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}$ otherwise.
5. Iterate between step (ii) and step (iv) until converge.

In this paper, we consider the logistic regression. We want to use Metropolis Hastings algorithm to fit the parameter of the model. i.e ${\beta }_0$ and ${\beta }_1$ in the following function.

$Y_i\sim Bin\mathrm{<mo\; stretchy="false">(</mo>1,}{\pi }_i\mathrm{<mo\; stretchy="false">)</mo>,}log\frac{{\pi }_i}{1-{\pi }_i}\mathrm{<mo>=</mo>}{\beta }_0+x_i{\beta }_1\mathrm{,}i\mathrm{<mo>=</mo>1,..,}n$ (1)

And the likelihood function is

$f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}{\beta }_0\mathrm{,}{\beta }_1\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{\displaystyle \prod _{i\mathrm{<mo>=</mo>1}}^n}{\left(\frac{e^{{\beta }_0+x_i{\beta }_1}}{1+e^{{\beta }_0+x_i{\beta }_1}}\right)}^{y_i}{\left(\frac{1}{1+e^{{\beta }_0+x_i{\beta }_1}}\right)}^{1-y_i}$ (2)

$2.0cm\mathrm{<mo>=</mo>}exp\mathrm{<mo>\{</mo>}n\overline{y}{\beta }_0+{\beta }_1{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n}{x}_i{y}_i-log\mathrm{<mo\; stretchy="false">(</mo>1}+e^{{\beta }_0+x_i{\beta }_1}\mathrm{<mo\; stretchy="false">)</mo><mo>\}</mo>}$ (3)

Consider the prior distribution of ${\beta }_0$ and ${\beta }_1$ as the independent normal distribution

${\beta }_j\sim N\mathrm{<mo\; stretchy="false">(</mo>}{\mu }_j\mathrm{,}{\sigma }_j^2\mathrm{<mo\; stretchy="false">)</mo>,}j\mathrm{<mo>=</mo>0,1.}$ (4)

Then the posterior distribution is

$f\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0\mathrm{,}{\beta }_1\mathrm{<mo>|</mo>}y\mathrm{<mo\; stretchy="false">)</mo>}\propto f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}{\beta }_0\mathrm{,}{\beta }_1\mathrm{<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0\mathrm{,}{\beta }_1\mathrm{<mo\; stretchy="false">)</mo>}$

$2.0cm\propto exp\mathrm{<mo>\{</mo>}n\overline{y}{\beta }_0+{\beta }_1{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n}\mathrm{<mo\; stretchy="false">[</mo>}{x}_i{y}_i-log\mathrm{<mo\; stretchy="false">(</mo>1}+e^{{\beta }_0+x_i{\beta }_1}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}$

$$2.0cm-\frac{{\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0-{\mu }_0\mathrm{<mo\; stretchy="false">)</mo>}}^2}{2{\sigma }_0^2}-\frac{{\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_1-{\mu }_1\mathrm{<mo\; stretchy="false">)</mo>}}^2}{2{\sigma }_1^2}\mathrm{<mo>\}</mo>}$$

And we need to get the Markov sequence of the parameters. We propose to consider three sampling methods based on the general Metropolis Hastings algorithm, which is independent sampling, dependent sampling and individual sampling

Independent sampling Assuming that the parameters ${\beta }_0$  and ${\beta }_1$ are independent, we can generate ${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}$  from the proposal distribution$N\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}diag\mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_0^2\mathrm{,}{\sigma }_1^2\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$ . Then the independent sampling can be stated as following.

Algorithm (Independent sampling)

1. Initialize${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}$ .
2. Generate  from the proposal distribution$N\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}diag\mathrm{<mo\; stretchy="false">(</mo>}{\sigma }_0^2\mathrm{,}{\sigma }_1^2\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$ .
3. Calculate acceptance probability
4.               $\alpha \mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}min\mathrm{<mo>\{</mo>}\frac{f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}{\beta }_{0^{\prime }}\mathrm{,}{\beta }_{1^{\prime }}\mathrm{<mo\; stretchy="false">)</mo>}q\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_{0^{\prime }}\mathrm{,}{\beta }_{1^{\prime }}\mathrm{<mo\; stretchy="false">)</mo>}}{f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}q\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,1<mo>\}</mo>}$
5. We accept ${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\beta }^{\prime }$  with probability $\alpha \mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}$ or ${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}$  otherwise.
6. Iterate between step (ii) and step (iv) until converge.

Dependent sampling If the parameters are not independent, the convergence of stochastic sequence may be less efficient. So we consider the multinomial prior for ${\beta }_0$ . In this case, we use the Fisher information matrix $H\mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{<mo\; stretchy="false">)</mo>}$  to build the proposal distribution.

               ${\beta }^{\prime }\sim q\mathrm{<mo>=</mo>}N\mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{,}{c}_{\beta }^2{\mathrm{<mo\; stretchy="false">[</mo>}H\mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}}^{-1}\mathrm{<mo\; stretchy="false">)</mo>}$     (5)

Where $c_{\beta }$  is a regulation parameter, which can be adjust to reach the target acceptance rate. Based on the likelihood function, we get the information matrix

  $H\mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{X}^{T}diag\mathrm{<mo\; stretchy="false">(</mo>}{h}_i\mathrm{<mo\; stretchy="false">)</mo>}X+{\Sigma }_{\beta }^{-1}$    (6)

Where${\Sigma }_{\beta }$  is the prior covariance matrix of $\beta$ , $h_i\mathrm{<mo>=</mo>}exp\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0+{\beta }_1{x}_i\mathrm{<mo\; stretchy="false">)</mo><mo>/</mo><mo\; stretchy="false">(</mo>1}+exp\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0+{\beta }_1{x}_i{\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}}^2$ , $X\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{1}_n\mathrm{,}x\mathrm{<mo\; stretchy="false">)</mo>}$  is a $2\times n$  matrix.

Then we give the algorithm for this method:

Algorithm (Dependent Sampling)

1. Initialize${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}$ .
2. Calculate Fisher Information matrix
3.                $diag\mathrm{<mo\; stretchy="false">(</mo>}{h}_i\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}diag\mathrm{<mo\; stretchy="false">(</mo>}\frac{exp\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}+{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}{x}_i\mathrm{<mo\; stretchy="false">)</mo>}}{{\mathrm{<mo\; stretchy="false">(</mo>1}+exp\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}+{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}{x}_i\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}}^2}\mathrm{<mo\; stretchy="false">)</mo>}$
4.               $H\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}{X}^{T}diag\mathrm{<mo\; stretchy="false">(</mo>}{h}_i\mathrm{<mo\; stretchy="false">)</mo>}X+{\Sigma }_{{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}}^{-1}{S}_{{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}}\mathrm{<mo>=</mo>}{c}_{{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}}^2{\mathrm{<mo\; stretchy="false">[</mo>}H\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">]</mo>}}^{-1}$

 

5. Generate  from the proposal distribution $N\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{S}_{{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}}\mathrm{<mo\; stretchy="false">)</mo>}$ .
6. Calculate acceptance probability
7.                 $\alpha \mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}min\mathrm{<mo>\{</mo>}\frac{f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}{\beta }_{0^{\prime }}\mathrm{,}{\beta }_{1^{\prime }}\mathrm{<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>}{\beta }_{0^{\prime }}\mathrm{,}{\beta }_{1^{\prime }}\mathrm{<mo\; stretchy="false">)</mo>}q\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo>|</mo>}{\beta }^{\prime }\mathrm{,}{S}_{{\beta }^{\prime }}\mathrm{<mo\; stretchy="false">)</mo>}}{f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}q\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\prime }\mathrm{<mo>|</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{S}_{{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}}\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,1<mo>\}</mo>}$
8. We accept  with probability $\alpha \mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}$  or ${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}$  otherwise.
9. Iterate between step (ii) and step (v) until converge.

Individual sampling

Algorithm (Individual Sampling)

1. Initialize${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>0<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}$ .
2. Generate  from the proposal distribution$N\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{\sigma }_0^2\mathrm{<mo\; stretchy="false">)</mo>}$ .
3. Set ${\beta }^{\prime }\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{\beta }_{0^{\prime }}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}$  and calculate acceptance rate
4.                ${\alpha }_0\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}min\mathrm{<mo>\{</mo>}\frac{f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}{\beta }_{0^{\prime }}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>}{\beta }_{0^{\prime }}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}{f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,1<mo>\}</mo>}$
5. We accept ${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\beta }^{\prime }$  with probability ${\alpha }_0\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}$  or ${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}$  otherwise.
6. Generate from the proposal distribution$N\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}{\sigma }_1^2\mathrm{<mo\; stretchy="false">)</mo>}$ .
7. Set${\beta }^{\prime }\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_{1^{\prime }}\mathrm{<mo\; stretchy="false">)</mo>}$  and calculate acceptance rate
8.                ${\alpha }_1\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}min\mathrm{<mo>\{</mo>}\frac{f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_{1^{\prime }}\mathrm{<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_{1^{\prime }}\mathrm{<mo\; stretchy="false">)</mo>}}{f\mathrm{<mo\; stretchy="false">(</mo>}y\mathrm{<mo>|</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}\pi \mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,1<mo>\}</mo>}$
9. We accept ${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t+\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\beta }^{\prime }$ with probability ${\alpha }_1\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}$  or ${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t+\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}$  otherwise.
10. Iterate between step (ii) and step (vii) until converge.

Stochastic approximation Monte Carlo computation (SAMC)

Consider the sampling distribution that$p\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}c{p}_0\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$ , where c is a constant, x is generated from the sample space. We let $E_1\mathrm{,...,}{E}_m$  donate  disjoint regions that from a partition of $\chi$ , which can be partitioned according to any function of x such as the energy function $U\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}$  as follows: $E_1\mathrm{<mo>=</mo><mo>\{</mo>}x\mathrm{<mo>:</mo>}U\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}\le u_1\mathrm{<mo>\}</mo>,}{E}_2\mathrm{<mo>=</mo><mo>\{</mo>}x\mathrm{<mo>:</mo>}{u}_1\mathrm{<mo><</mo>}U\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}\le u_2\mathrm{<mo>\}</mo>,}\dots \mathrm{,}{E}_m\mathrm{<mo>=</mo><mo>\{</mo>}x\mathrm{<mo>:</mo>}U\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}\ge u_m\mathrm{<mo>\}</mo>}$ . Then we let  donate the estimate of  and rewritten the invariant distribution  in the generalized Wang-Landau (GWL) algorithm as

                $p_{{\theta }_t}\mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}\propto {\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^m}\frac{\psi \mathrm{<mo\; stretchy="false">(</mo>}x\mathrm{<mo\; stretchy="false">)</mo>}}{e^{{\theta }_{ti}}}I\mathrm{<mo\; stretchy="false">(</mo>}x\in E_i\mathrm{<mo\; stretchy="false">)</mo>}$     (7)

For theoretical simplicity, we assume that  for all t, where  is a compact set. Let  be an m-vector with  and , which defines the desired sampling frequency for each of the sub regions. Henceforth,  is called the desired sampling distribution. Let  be a positive, no decreasing sequence satisfying

$\mathrm{<mo\; stretchy="false">(</mo>}a\mathrm{<mo\; stretchy="false">)</mo>}{\displaystyle \sum _{t\mathrm{<mo>=</mo>1}}^{\infty }}{\gamma }_t\mathrm{<mo>=</mo>}\infty and \mathrm{<mo\; stretchy="false">(</mo>}b\mathrm{<mo\; stretchy="false">)</mo>}{\displaystyle \sum _{t\mathrm{<mo>=</mo>1}}^{\infty }}{\gamma }_t\mathrm{<mo><</mo>}\infty$ (8)

For some. For example, in this article we set

${\gamma }_t\mathrm{<mo>=</mo>}\frac{t_0}{max\mathrm{<mo\; stretchy="false">(</mo>}{t}_0\mathrm{,}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{.}t\mathrm{<mo>=</mo>1,2...,}$ (9)

For some specified value of .

In the logistic regression model, the invariant distribution $p_{{\theta }_t}\mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{<mo\; stretchy="false">)</mo>}\propto exp\mathrm{<mo>\{</mo>}n\overline{y}{\beta }_0+{\beta }_1{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>1}}^n}{x}_i{y}_i-log\mathrm{<mo\; stretchy="false">(</mo>1}+e^{{\beta }_0+x_i{\beta }_1}\mathrm{<mo\; stretchy="false">)</mo><mo>\}</mo>}$ and the proposal distribution $q\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}N\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t-\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{,}{1}_2\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$ . With the foregoing notation, one iteration of SAMC can be described as follows:

Algorithm (SAMC)

Generate a sample ${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo><mo\; stretchy="false">(</mo>}{\beta }_0^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}{\beta }_1^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}$ by a single MH update with the target distribution given by

$p_{{\theta }_t}\mathrm{<mo\; stretchy="false">(</mo>}\beta \mathrm{<mo\; stretchy="false">)</mo>}\propto exp\mathrm{<mo>\{</mo>}n\overline{y}{\beta }_0+{\beta }_1{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n}{x}_i{y}_i-log\mathrm{<mo\; stretchy="false">(</mo>1}+e^{{\beta }_0+x_i{\beta }_1}\mathrm{<mo\; stretchy="false">)</mo><mo>\}</mo>}$

Generateaccording to the proposal distribution $q\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\prime }\mathrm{<mo>|</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo><mo>=</mo>}N\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{,}{1}_2\mathrm{<mo\; stretchy="false">)</mo><mo\; stretchy="false">)</mo>}$ . If$J\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}\notin S$ , then update  to, where  denote the collection of indices of the sub regions from which a sample has been proposed and  denote the index of the sub region of sample.

Calculate the ratio

$r\mathrm{<mo>=</mo>}{e}^{{\theta }_{tJ\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}-{\theta }_{tJ\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}}}\frac{\psi \mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}q\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo>|</mo>}{\beta }^{\prime }\mathrm{<mo\; stretchy="false">)</mo>}}{\psi \mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}q\mathrm{<mo\; stretchy="false">(</mo>}{\beta }^{\prime }\mathrm{<mo>|</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}\mathrm{<mo\; stretchy="false">)</mo>}}$

Accept the proposal with probability. If accepted ${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t+\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\beta }^{\prime }$ , set, otherwise, set ${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t+\mathrm{1<mo\; stretchy="false">)</mo>}}\mathrm{<mo>=</mo>}{\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t\mathrm{<mo\; stretchy="false">)</mo>}}$ .

For all, update to as follows:

${\theta }_{t+\mathrm{1,}i}\mathrm{<mo>=</mo>}{\theta }_{t\mathrm{,}i}+{\gamma }_{t+1}\mathrm{<mo\; stretchy="false">(</mo>}{e}_{t+\mathrm{1,}i}-{\pi }_i\mathrm{<mo\; stretchy="false">)</mo>}$

where $e_{t+\mathrm{1,}i}\mathrm{<mo>=</mo>1}$ if ${\beta }^{\mathrm{<mo\; stretchy="false">(</mo>}t+\mathrm{1<mo\; stretchy="false">)</mo>}}\in E_i$ , and 0 otherwise.

 In this part, we generate simple logistics regression data and use the Monte Carlo method to fit the model.

                $log\frac{{\pi }_i}{1-{\pi }_i}\mathrm{<mo>=</mo>}{\beta }_0+x_i{\beta }_1\mathrm{,}i\mathrm{<mo>=</mo>1,..,}n$              (10)

 Based on the model and given the value of  and, we generate 1000 samples of.

                $x_i\sim N\mathrm{<mo\; stretchy="false">(</mo>1,1<mo\; stretchy="false">)</mo>}{\pi }_i\mathrm{<mo>=</mo>}\frac{exp\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0+x_i{\beta }_1\mathrm{<mo\; stretchy="false">)</mo>}}{1+exp\mathrm{<mo\; stretchy="false">(</mo>}{\beta }_0+x_i{\beta }_1\mathrm{<mo\; stretchy="false">)</mo>}}i\mathrm{<mo>=</mo>1,...}n$          (11)

 Then  is generated from binomial  distribution. We do iterations for 1000 times and calculate the mean and variance of the$\mathrm{<mo\; stretchy="false">(</mo>}{\widehat{\beta }}_0\mathrm{,}{\widehat{\beta }}_1\mathrm{<mo\; stretchy="false">)</mo>}$  in different methods (Table 1).

Based on the Table 1, we find that both MH algorithms and SAMC have good performance on calculating the parameters in logistics regression while SAMC has smallest variance which indicates a better convergence rate and robustness. Since the independence of  also the  is close to 0, all the three MH algorithms perform similarly. As for SAMC algorithm, we have to partition the sample space based on the energy function. The energy regions can set up the initial values spanned in the full model space and therefore can have faster convergence rate and lower variance than other MH algorithm. Finally, we use the MLE method, i.e the ’glm’ function in R to estimate the parameters for comparison. It shows that SAMC algorithm even achieve better performance than the MLE method in terms of the variance of the estimators. Therefore, both MH and SAMC algorithms are eligible to obtain the consistent estimators in the logistic regression but even for this simple problem, SAMC still can enjoy better performance than others.

 In this part, we apply Monte Carlo methods into a real example data set. Our dataset is obtianed from the Wechsler Adult Intelligence Scale (WAIS), which is a test designed to measure intelligence in adults and older adolescents and it is currently in its fourth edition (WAIS-IV). The original WAIS was published in February 1955 by David Wechsler, as a revision of the Wechsler-Bellevue Intelligence Scale that had been released in 1939. The propose of this test is to find the relationship between intelligence and Senile Dementia among older people. In our score scale, people (with and without Senile Dementia) were asked to do some tasks individually. When all tasks completed, they will get a score, ranging from 0-20. Therefore, the response  is a binary response that reflects whether people have Senile Dementia. The WAIS score is the explanation variables. Based on data, we need to build the logistic model and estimate the parameters.

Firstly, we apply the MH Independent Sampling Algorithm. We set the length of chain equals to 10000 and get the plot of parameters of, . As showed in the lower two figures of Figure 1, the two lines does not coincide, which indicate that independent sampling algorithm does not converge. The independent sampling algorithm assumes that parameters are independent; however, it cannot be easily satisfied in the real application. We calculate the covariance of  and, which is . Therefore, the independent sampling method cannot be applied here. Then we use another two MH algorithms, i.e dependent sampling and individual sampling as well as SAMC algorithm. We set the length of chain equals to 5000 and plot the cumulative mean curves for each method. As showed in Figure 2, all the methods converge when iteration grows while the SAMC converges faster than other, i.e; the two curves for SAMC coincide at around iteration 2000 while others need more iteration. To access the accuracy of estimations of  and, we use the results from MLE method as benchmark as compare them with the ones from the Monte Carlo methods. As showed in Table 2, the results from SAMC are more close to the MLE and also it has the smallest variance among all other methods. According to the result,  is significantly less than 0, which means people who get less score have higher risk of Senile Dementia.

Figure 1 The path and mean plots of , for MH Independent Sampling algorithm. The upper two plots are the trace plot of estimators of and at each iteration. The lower two plots are the cumulative mean of and , respectively. The solid line calculates the cumulative mean starting at iteration 1, while the dash lines starts at step 1500.

Figure 2 Comparison the cumulative mean plot of and for multiple Monte Carlo algorithms. The upper two panels denote the MH Dependent Sampling method; the middle two panels denote the MH individual sampling method; the lower two panels denote the SAMC method. The solid line calculates the cumulative mean starting at iteration 1, while the dash lines starts at step 1500.

In this paper, we have proposed three MH algorithms and also SAMC to estimate the parameters of logistic regression. According to the simulation study, when the parameters are uncorrelated, each method has good performance. When the parameters are highly correlated, Independent Sampling method is less efficient and even cannot converge. The results showed that SAMC algorithms can be used in both scenarios and always outperform others. Compared with other MH algorithm, SAMC need less iterations to converge and have smaller variance of estimators even compared with MLE method.

One interesting further work is to apply the Monte Carlo simulations to high-dimensional big data problem. The challenge comes from two aspects. First, the data contains complex structure and therefore, the acceptance rate for most MCMC algorithm can be slow and therefore can be less efficient. Although SAMC has been proved to be a good method for the complex data and can overcome some local trap problem, it still needs to improve to the high-dimensional case where the likelihood function may not be intractable. Secondly, the big data have catastrophically large volumes and most single chain Monte Carlo simulations cannot apply. At each iteration, it should go through all the data which need large memory space for the computer machine and long time to generate the data at each iteration. For the big volume of data, we should consider parallel computing strategies, i.e, generate multiple Monte Carlo chains in parallel and then integrate chains to get the single estimator (Liang and Wu, 2013). Also, the high speed graphics processing unit (GPU) can be helpful for accelerating the speed of many MCMC algorithms.

In conclusion, we introduced several Monte Carlo simulation methods to estimate parameters of generalized linear model, i.e logistic regression. This problem can also be achieved by the MLE method which might be easier to calculate. However, our Monte Carlo algorithm can obtain a chain of estimations instead of a single value and therefore can provide more information about the process of reaching the true values of parameters and can be useful for further analysis.

The authors would like to thank the referee and the Editor for careful reading and comments which greatly improved the paper.

Author declares that there is no conflict of interest.

1. Metropolis N, Rosenbluth AW, Rosenbluth MN, et al. Equation of state calculations by fast computing machines. The journal of chemical physics. 1953;21:1087–1092.
2. Hastings WK. Monte Carlo sampling methods using Markov chains and their applications. Biometrika. 1970;57(1):97–109.
3. Müller P. A generic approach to posterior integration and Gibbs sampling. USA: Purdue University, Department of Statistics. 1991.
4. Tierney L. Markov chains for exploring posterior distributions. The Annals of Statistics. 1994;22(4):1701–1728.
5. Chib S, Greenberg E. Markov chain Monte Carlo simulation methods in econometrics. Econometric theory. 1996;12:409–431.
6. Liang Liu, Carroll, Carroll RJ. Stochastic approximation in Monte Carlo computation. Journal of the American Statistical Association. 2007;102:305–320.
7. Geman S, Geman D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on pattern analysis and machine intelligence. 1994;721–741.
8. Geyer CJ. Markov chain Monte Carlo maximum likelihood. 1991.
9. Hukushima K, Nemoto K. Exchange Monte Carlo method and application to spin glass simulations. Journal of the Physical Society of Japan. 1996;65:1604–1608.
10. Marinari E, Parisi G. Simulated tempering: a new Monte Carlo scheme. Europhysics Letters. 1992;19:451.
11. Geyer CJ, Thompson EA. Annealing Markov chain Monte Carlo with applications to ancestral inference. Journal of the American Statistical Association. 1995;90:909–920.
12. Liang F, Wong WH. Evolutionary Monte Carlo for protein folding simulations. The Journal of Chemical Physics. 2001;115:3374–3380.
13. Wong WH, Liang F. Dynamic weighting in Monte Carlo and optimization. Proceedings of the National Academy of Sciences. 1997;94:14220–14224.
14. Berg BA, Neuhaus T. Multicanonical ensemble: A new approach to simulate first–order phase transitions. Physical Review Letters. 1992;68:1–9.
15. Hesselbo B, Stinchcombe RB. Monte Carlo simulation and global optimization without parameters. Physical review letters. 1995;74:2151.
16. Wang F, Landau D. Efficient, multiple–range random walk algorithm to calculate the density of states. Physical review letters. 2001;86:2050.
17. Mitsutake A, Sugita Y, Okamoto Y. Replica–exchange multicanonical and multicanonical replica–exchange Monte Carlo simulations of peptides. I. Formulation and benchmark test. The Journal of chemical physics. 2003;118:6664–6675.