## 2 Materials and Methods
### Numerical Phantom  
To examine when the iterative scheme fails by being trapped by a local minimum and how the Markov chain Monte Carlo is capable of escaping from such local minima, a toy model is devised which is simple enough to explicitly understand the structure of the cost function but is complicated enough that the cost function has one local minimum and one global minimum.  
We consider diffuse optical tomography in the two-dimensional space. In our toy model we suppose that the diffusion coefficient \(D\) is constant everywhere in the medium but there is absorption inhomogeneity and the absorption coefficient \(\delta\mu_{a}\) is unknown. Moreover, we assume that \(\delta\mu_{a}(\mathbf{\rho})\) is given by  
\[\delta\mu_{a}(\mathbf{\rho})=\eta f_{a}(x)\delta(y-y_{0}), \tag{35}\]  
where \(\eta,y_{0}\) are given positive constants. Here, \(f_{a}(x)\) is given by  
\[f_{a}(x)=\left[a^{3}+3\left(1+\frac{\tanh x^{2}}{10}\right)a^{2}\right]\left( 1-\tanh x^{2}\right), \tag{36}\]  
where \(a>1.1\) is a constant. Thus, \(a\) is the parameter to be reconstructed in our toy inverse problem. As is shown in Figure 2, the function \(f_{a}(x)\) has one peak at \(x=0\) and the maximum value is \(f_{a}(0)=a^{2}(a+3)\), \(f_{a}\) monotonically decays for \(|x|>0\), and \(f_{a}\to 0\) as \(|x|\to\infty\).  
Now we describe how the forward data is computed. The unit of length and unit of time are taken to be mm and ps, respectively. On the \(x\)-axis, we place two sources (\(M_{s}=2\)) at \(\mathbf{\rho}_{s}^{1}=(-20,0)\), \(\mathbf{\rho}_{s}^{2}=(20,0)\) and three detectors (\(M_{d}=3\)) at \(\mathbf{\rho}_{d}^{1}=(-40,0)\), \(\mathbf{\rho}_{d}^{2}=(0,0)\), \(\mathbf{\rho}_{d}^{3}=(40,0)\). We set \(\mu_{s}^{\prime}=1\), \(\mu_{a}=0.02\), \(n=1.37\). Suppose that there is absorption inhomogeneity at depth \(5\). For \(\delta\mu_{a}\), we put \(\eta=300/c\), \(y_{0}=5\), and  
\[a=1.5. \tag{37}\]  
To distinguish, hereafter the true value of \(a\) is denoted by \(\bar{a}\). We assume 3% Gaussian noise and give the measured data as  
\[y_{ijk}=u(\mathbf{\rho}_{d}^{j},t^{k};\mathbf{\rho}_{s}^{i};\bar{a})(1+e_{ijk}), \tag{38}\]  
where \(e_{ijk}\sim\mathcal{N}(0,0.03^{2})\). In this numerical experiment, \(M_{t}=500\) and \(t^{k+1}-t^{k}=t^{1}=5\) (\(k=1,\ \ldots,\ M_{t}-1\)).  
Figure 1: Measured data from TRS-80. The instrument response function and measured reflected light are shown in the upper and lower panels, respectively. In the lower panel, the dashed lines show \(t_{1}=2.00\,\mathrm{ns}\) and \(t_{M_{t}}=8.00\,\mathrm{ns}\) (\(M_{t}=601\)).  
By substituting the form (35) of \(\delta b\) in (21), we can express the energy density as  
\[\begin{split} u(\rho_{d^{\prime}}^{j}t^{k};\rho_{s^{\prime}}^{i}a)& =u(\rho_{d^{\prime}}^{j}t^{k};\rho_{s}^{i})\\ &=u_{0}(\rho_{d^{\prime}}^{j}t^{k};\rho_{s}^{i})\exp\left[-\frac{ \eta e^{-\mu_{0}ct^{k}}}{u_{0}(\rho_{d^{\prime}}^{j}t^{k};\rho_{s}^{i})}\int_{0 }^{t^{k}}g(0,t^{k};y_{0},s)\right.\\ &\times g(y_{0},s;0,0)\left.\left(\int_{-\infty}^{\infty}f_{a}(x^ {\prime})e^{-\frac{(\sigma_{a}^{j}-x^{\prime})^{2}}{4 Dc^{k}(k-\rho_{s})}}e^{- \frac{(x^{\prime}-x^{\prime}_{s})^{2}}{4 Dcs}}\,dx^{\prime}\right)\ ds\right],\end{split} \tag{39}\]  
where  
\[u_{0}=\frac{e^{-\mu_{0}ct^{k}}}{2\pi Dt^{k}}e^{-\frac{(x^{\prime}_{s}-x^{ \prime}_{s})^{2}}{4 Dct^{k}}}\left[1-\frac{\sqrt{\pi Dtct^{k}}}{\ell}e^{\left( \frac{\sqrt{Dct^{k}}}{\ell}\right)^{2}}\text{erfc}\left(\frac{\sqrt{Dct^{k}}} {\ell}\right)\right]. \tag{40}\]