Gaussian processes (GPs) provide a Bayesian framework for regression that allows exact predictive analysis using matrix operations. This paper introduces GP priors over functions, enabling exact Bayesian inference for fixed hyperparameters. Two methods—optimization and averaging via Hybrid Monte Carlo (HMC)—are tested on challenging problems, yielding excellent results. The paper discusses the use of GPs for regression, focusing on parameterizing covariance functions and estimating hyperparameters from data. It shows how hyperparameters can be estimated using maximum likelihood or Bayesian approaches, leading to "Automatic Relevance Determination" (ARD). The predictive distribution for a test case is derived from the joint Gaussian distribution of training and test outputs, conditioned on observed targets. The predictive mean and variance are given by equations involving the covariance matrix and training data. The paper also presents a covariance function with multiple terms, including a linear regression term and noise term, and discusses its validity for various input dimensions. The log likelihood of the training data is given by a formula involving the determinant of the covariance matrix and the training targets. The paper compares GPs to previous work, including ARMA models, spline smoothing, and regularization networks. It describes two methods for training GPs: maximum likelihood and HMC. The HMC method uses a prior distribution over hyperparameters and samples from the posterior distribution using a Markov chain Monte Carlo approach. The paper presents experimental results on a modified robot arm problem and five real-world datasets, showing that GPs perform well compared to other methods. The results indicate that GPs can achieve performance close to the theoretical minimum error. The paper concludes by discussing future directions, including the use of GPs for classification and non-stationary covariance functions. The authors hope to make their GP prediction code publicly available.Gaussian processes (GPs) provide a Bayesian framework for regression that allows exact predictive analysis using matrix operations. This paper introduces GP priors over functions, enabling exact Bayesian inference for fixed hyperparameters. Two methods—optimization and averaging via Hybrid Monte Carlo (HMC)—are tested on challenging problems, yielding excellent results. The paper discusses the use of GPs for regression, focusing on parameterizing covariance functions and estimating hyperparameters from data. It shows how hyperparameters can be estimated using maximum likelihood or Bayesian approaches, leading to "Automatic Relevance Determination" (ARD). The predictive distribution for a test case is derived from the joint Gaussian distribution of training and test outputs, conditioned on observed targets. The predictive mean and variance are given by equations involving the covariance matrix and training data. The paper also presents a covariance function with multiple terms, including a linear regression term and noise term, and discusses its validity for various input dimensions. The log likelihood of the training data is given by a formula involving the determinant of the covariance matrix and the training targets. The paper compares GPs to previous work, including ARMA models, spline smoothing, and regularization networks. It describes two methods for training GPs: maximum likelihood and HMC. The HMC method uses a prior distribution over hyperparameters and samples from the posterior distribution using a Markov chain Monte Carlo approach. The paper presents experimental results on a modified robot arm problem and five real-world datasets, showing that GPs perform well compared to other methods. The results indicate that GPs can achieve performance close to the theoretical minimum error. The paper concludes by discussing future directions, including the use of GPs for classification and non-stationary covariance functions. The authors hope to make their GP prediction code publicly available.