## 1 Introduction  
In near-infrared spectroscopy (NIRS), we estimate optical properties of biological tissue by solving inverse diffuse problems [1; 2]. Such inverse problems are commonly solved by means of iterative methods. In the case of a homogeneous medium, absorption and reduced scattering coefficients of the medium can be obtained with iterative methods such as the Levenberg-Marquardt algorithm [3; 4] from time-resolved measurements (for example, see the review article [5]). Neuroimaging for the human brain by NIRS through the neurovascular coupling has been developed and is called functional NIRS (fNIRS) [6; 7; 8; 9]. Since the region of interest on the head can be small, it is possible to assume a simple geometry of the half space. Quantitative measurements of inter-regional differences in neuronal activity requires accurate estimates of optical properties. In Ref. [10], the Nelder-Mead simplex method was used to retrieve optical parameters in layered tissue. Heterogeneity of optical properties can be obtained by diffuse optical tomography [1; 11; 12]. Iterative numerical schemes are used to minimize the cost function when solving these inverse problems. A gradient-based approach [13] was used to detect breast cancer [14]. The brain activity of a newborn infant was investigated [15] with diffuse optical tomography by TOAST (temporal optical absorption and scattering tomography) [16; 17], in which iterative algorithms such as the nonlinear conjugate gradients, damped Gauss-Newton method, and Levenberg-Marquardt method are implemented. Diffuse optical tomography was performed on human lower legs and a forearm with the algebraic reconstruction algorithm in the framework of the modified generalized pulse spectrum technique [18]. See Refs. [19; 20] for numerical techniques for diffuse optical tomography. For these iterative numerical schemes to work, it is important to choose a good initial guess for the initial value of the iteration. Solving inverse problems by the Bayesian approach has been sought as an alternative way. In Ref. [21], a novel use of the Bayesian approach was considered to take modeling error into account. The Bayesian inversion with the Metropolis-Hastings Markov chain Monte Carlo was used in Refs. [22; 23]. The Bayesian approach was used to determine optical parameters of the human head [24]. Although the Markov chain Monte Carlo (MCMC) approach is in principle able to escape from local minima, it is computationally time consuming. Hence, despite the above-mentioned efforts, the use of Monte Carlo for inverse problems in NIRS has been extremely limited. In this paper, we shed light on Markov chain Monte Carlo once again by combining it with an iterative scheme, and investigate the use of it in NIRS. In particular, we test a hybrid numerical scheme of Markov chain Monte Carlo and an iterative method. In the proposed hybrid scheme, the Markov chain Monte Carlo algorithm is first used to provide an initial guess using jumps in the landscape of the cost function, and then an iterative method is used after the initial large fluctuation. Thus, the proposed method realizes a fast reconstruction while, in the beginning, obtained values at Monte Carlo steps jump around in the landscape of the cost function. The MCMC simulation is necessary only for the first 100 steps. Then the hybrid method starts to use an iterative scheme and reconstructs optical properties by searching the global minimum. The computation time of the iterative scheme is negligible compared with that of the Monte Carlo simulation. Since the Monte Carlo simulation can be stopped way before the burn-in time for the hybrid method, the proposed scheme is at least ten times faster than simulated annealing, which is implemented by the naive use of the Metropolis-Hastings Markov chain Monte Carlo. The rest of the paper is organized as follows. We develop diffusion theory in Section 2. A polyurethane-based phantom and numerical phantom are also described in Section 2. Section 3 is devoted to experimental and numerical results.