Constructing an occupancy representation of the environment is a fundamental problem for robot autonomy. Many accurate and efficient methods exist that address this problem but most assume that the occupancy states of different elements in the map representation are statistically independent. The focus of this paper is to provide a model that captures correlation of the occupancy of map elements. Correlation is important not only for improved accuracy but also for quantifying uncertainty in the map and for planning autonomous mapping trajectories based on the correlation among known and unknown areas. Recent work proposes Gaussian Process (GP) regression to capture covariance information and enable resolution-free occupancy estimation. The drawback of techniques based on GP regression (or classification) is that the computation complexity scales cubically with the length of the measurement history. Our main contribution is a new approach for occupancy mapping that models the binary nature of occupancy measurements precisely, via a Bernoulli distribution, and provides an efficient approximation of GP classification with complexity that does not scale with time. We prove that the error between the estimates provided by our method and those provided by GP classification is negligible. The proposed method is evaluated using both simulated data and real data collected using a Velodyne Puck 3-D range sensor.