•Objective method for crustal modeling based on Bayesian inference is developed and applied to Japan arc.•Biases in conventional geochemical methods are investigated and a new bias-free rock composition model is presented.•Geo-neutrino flux from the Japan crust is calculated with fully probabilistic uncertainty estimation. Geoneutrino observations, first achieved by KamLAND in 2005 and followed by Borexino in 2010, have accumulated statistics and improved sensitivity for more than ten years. The uncertainty of the geoneutrino flux at the surface is now reduced to a level small enough to set useful constraints on U and Th abundances in the Bulk Silicate Earth (BSE). However, in order to make inferences on earth’s compositional model, the contributions from the local crust need to be understood within a similar uncertainty. Here we develop a new method to construct a stochastic crustal composition model utilizing Bayesian inference. While the methodology has general applicability, it incorporates all the local uniqueness in its probabilistic framework. Unlike common approaches for this type of problem, our method does not depend on crustal segmentation into upper, (middle) and lower, whose classification and boundaries are not always well defined. We also develop a new modeling method to infer rock composition distributions that conserve mass balance and therefore do not bias the results. Combined with a new vast collection of geochemical data for rock samples in the Japan arc, we apply this method to geoneutrino observation at Kamioka, Japan. Currently a difficulty remains in the handling of correlations in the flux integration; we conservatively assume maximum correlation, which leads to large flux estimation errors of 60–70%. Despite the large errors, this is the first local crustal model for geoneutrino flux prediction with probabilistic error estimation in a reproducible way.