The problem to compute the multiple integral \(Z_ K=\int^{b}_{a}...\int^{b}_{a}f(x_ 1,...,x_ K)dx_ 1...dx_ K\) is considered for some constants a and b, and large K, the multiplicity of the integral. The crude Monte Carlo method does not work well for large K. For this reason, the investigation is interested in estimating \(\log (Z_ K)\) directly rather than \(Z_ K\) itself. This way is motivated by a previous paper of the author [A Monte Carlo method for the objective Bayesian procedure, Res. Memorandum No.347, Inst. Statistical Math., Tokyo (1988)] which provides a solution to the objective Bayesian procedure. A new numerical integration method is proposed. This method is an appropriate one for very high dimensional functions, while its implementation is based on the Metropolis Monte Carlo algorithm. The computation of \(\log (Z_ K)\) is reduced to a simple integration of a certain statistical function with respect to a scale parameter over the range of unit interval. This new method ensures a substantial improvement in the accuracy comparing to the crude Monte Carlo integration. Results of some numerical experiments are given. A numerical example illustrates how the high dimensional integration on the infinite domain can be reasonably calculated. A FORTRAN program for estimating \(\log (Z_ K)\) is also presented.