Probabilistic models were developed to provide predictive distributions of daily maximum surface level ozone concentrations, which were calculated as the daily maximum of 8-hour running average of hourly data. Several forecast models were compared at two stations (Chilliwack and Surrey) in the Lower Fraser Valley of British Columbia, Canada, with NCEP's Global Forecasting System (GFS) Reforecast data used as predictors. The models were of two types, conditional density models and Bayesian models. The conditional density models used neural networks to estimate the parameters of the given distribution (the Gaussian distribution or the lower-bounded Johnson distribution) as functions of the predictor variables. The Bayesian models include both the Gaussian process model (with nonlinear or linear kernel) and the Bayesian neural network model. The Bayesian models (especially the Gaussian process model with nonlinear kernel) gave better forecasts for extreme events, namely poor air quality events defined as having ozone concentration >= 65 ppb.