This paper proposes applying Bayesian inference to the problem of estimating shallow S wave velocity structures-in particular, those containing a structural singularity (a low-velocity layer sandwiched between high-velocity layers or a high-velocity layer between low-velocity layers)-by using Rayleigh-wave phase velocities obtained in a microtremor survey. The method proposed here uses the empirical Bayesian approach, whereby field measurement data are used to build the prior distribution. An optimal velocity structure model is selected, from among multiple candidate models with different layer thickness parameter values, on the basis of a Bayes factor. The use of the proposed method allows the layer thickness parameter values in a one-dimensional velocity structure model to be determined objectively for explaining an observed Rayleigh-wave phase velocity dispersion curve. It also stabilizes the inversion process by the use of the prior distribution and strengthens constraints on the unknown parameter, or S wave velocity in each layer, by using multiple Rayleigh-wave modes. These characteristics of the proposed method allow flexible and stable inverse analysis of highly variable shallow subsurface structures. The present study conducts numerical experiments to show that our method allows the assumed structural singularities are reproduced. Then field measurement data from three sites are used to show that the method has allowed shallow structural singularities (depths 5-30 m) to be identified appropriately.Research on model selection in microtremor array methods has only started in recent years. Giulio et al. (2012) applied a quantity called the corrected Akaike information criterion (c-AIC;Hurvich & Tsai, 1989;Sugiura, 1978) to selecting a velocity structure model. When identifying parameter values that maximized the likelihood, however, they set limits on their range of parameter search and did not allow velocity reversal layers in some cases. That is not an appropriate way to use the c-AIC, because there is no way to count the number of "free" parameters, which is needed to calculate the c-AIC, under those limitations. Molnar et al. (2010), by contrast, avoided the problem by setting limits on their range of search and using those limits as the prior distribution for Bayesian inference. Dettmer et al. (2012) adopted a Bayesian CHO AND IWATA 527