In this paper, we provide a novel regression algorithm based on a Gaussian random field (GRF) indexed by a Riemannian manifold (M, g). We utilize both the labeled and unlabeled data sets to exploit the geometric structure of M. We use the recovered heat (H) kernel as the covariance function for the GRF (HGRF). We propose a Monte Carlo integral theorem on Riemannian manifolds and derive the corresponding convergence rate and approximation error. Based on this theorem, we correctly normalize the recovered eigenvector to make it compatible with Riemannian measure. More importantly, we prove that the HGRF is intrinsic to the original data manifold by comparing the pullback geometry and the original geoemtry of M. Essentially it is a semi-supervised learning method, which means the unlabeled data can be utilized to help identify the geometry structure of the unknown manifold M.