Goal: Write down a probability model for the data-generating process for \(\mathbf{y}\).
Direct statistical model, \[\mathbf{y} \sim \mathcal{D}(\theta).\]
Model for the residuals of a numerical model, \[\mathbf{r} = \mathbf{y} - F(\mathbf{x}) \sim \mathcal{D}(\theta).\]
Model Fitting as Maximum Likelihood Estimation
We can interpret fitting a model (reducing error according to some loss or error metric) as maximizing the probability of observing our data from this data generating process.
Accounting for Model Residuals
The Energy Balance Model Revisited
Last class, we introduced the Energy Balance Model (EBM):
If we ignore all of the constants, this is proportional to the negative mean-squared error.
Maximum Likelihood Estimate
# p are the model parameters, σ the standard deviation of the normal errors, y is the data, m the model functionfunctionlog_likelihood(p, σ, y, m) y_pred =m(p) ll =sum(logpdf.(Normal.(y_pred, σ), y))endebm_wrap(params) =ebm(forcing_non_aerosol_85[idx], forcing_aerosol_85[idx], p = params)# maximize log-likelihood within some range# important to make everything a Float instead of an Int lower = [1.0, 50.0, 0.0, 0.0]upper = [4.0, 150.0, 2.0, 10.0]p0 = [2.0, 100.0, 1.0, 1.0]result = Optim.optimize(params ->-log_likelihood(params[1:end-1], params[end], temp_obs, ebm_wrap), lower, upper, p0)θ = result.minimizer
# set number of sampled simulationsn_samples =1000residuals =rand(Normal(0, θ[end]), (n_samples, length(temp_obs)))model_out =ebm_wrap(θ[1:end-1])# this uses broadcasting to "sweep" the model simulation across the sampled residual matrixmodel_sim = residuals .+ model_out'# need to transpose the model output vector due to how Julia treats vector dimensions
Figure 3: Comparison of best fit with uncertain realization for the EBM with normal residuals.
Checking for Calibration
For an \(\alpha\)% projection interval \(\mathcal{I}_\alpha\), we would expect ~\(\alpha\)% of the data to be contained in this interval.
This rate is quantified by the surprise index: \(1 - \frac{1}{n} \sum_{i=1}^n \mathbb{I}_{\mathcal{I}_\alpha}(y_i).\)
Surprise Index
Code
surprises =0# initialize surprise counter# go through the data and check which points are outside of the 90% intervalfor i =1:length(temp_obs)## The || operator is an OR, so returns true if either of the terms are trueif (temp_obs[i] < q_90[1, i]) || (q_90[2, i] < temp_obs[i]) surprises +=1endendsurprises /length(temp_obs)
# p are the model parameters, σ the standard deviation of the AR(1) errors, ρ is the autocorrelation coefficient, y is the data, m the model functionfunctionar_log_likelihood(p, σ, ρ, y, m) y_pred =m(p) ll =0# initialize log-likelihood counter residuals = y_pred .- y ll +=logpdf(Normal(0, σ/sqrt(1-ρ^2)), residuals[1])for i =1:length(residuals)-1 residuals_whitened = residuals[i+1] - ρ * residuals[i] ll +=logpdf(Normal(0, σ), residuals_whitened)endreturn llendebm_wrap(params) =ebm(forcing_non_aerosol_85[idx], forcing_aerosol_85[idx], p = params)# maximize log-likelihood within some range# important to make everything a Float instead of an Int lower = [1.0, 50.0, 0.0, 0.0, -1.0]upper = [4.0, 150.0, 2.0, 10.0, 1.0]p0 = [2.0, 100.0, 1.0, 1.0, 0.0]result = Optim.optimize(params ->-ar_log_likelihood(params[1:end-2], params[end-1], params[end], temp_obs, ebm_wrap), lower, upper, p0)θ_ar1 = result.minimizer
Figure 9: Comparison of best fit with uncertain realization for the EBM with normal residuals.
AR(1) Surprise Index
Code
surprises =0# initialize surprise counter# go through the data and check which points are outside of the 90% intervalfor i =1:length(temp_obs)## The || operator is an OR, so returns true if either of the terms are trueif (temp_obs[i] < q90_ar1[1, i]) || (q90_ar1[2, i] < temp_obs[i]) surprises +=1endendsurprises /length(temp_obs)