n <- 1000
split <- 0.8
# Auxiliary functions
mu <- function(x1){
10 + 5*x1^2
}
sigma_v <- function(x1){
30*x1
}
# Generating heteroscedastic data.
x <- runif(n, 1, 10)
y <- rnorm(n, mu(x), sigma_v(x))
# Train set
x_train <- x[1:(n*split)]
y_train <- y[1:(n*split)]
# Calibration/Validation set.
x_cal <- x[(n*split+1):n]
y_cal <- y[(n*split+1):n]
# New observations or the test set.
x_new <- runif(n/5, 1, 10)
# Fitting a simple linear regression, which will not capture the heteroscedasticity
model <- lm(y_train ~ x_train)
y_hat_cal <- predict(model, newdata=data.frame(x_train=x_cal))
MSE_cal <- mean((y_hat_cal - y_cal)^2)
y_hat_new <- predict(model, newdata=data.frame(x_train=x_new))
pit <- PIT_global(ycal=y_cal, yhat= y_hat_cal, mse=MSE_cal)
recalibrate(
space_cal=x_cal,
space_new=x_new,
yhat_new=y_hat_new,
pit_values=pit,
mse= MSE_cal,
type="local")
Run the code above in your browser using DataLab