From 733312a869fe28f87e2187d2f9cfe8a10e1ded72 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sun, 9 Jun 2024 05:49:42 +0200 Subject: [PATCH 01/46] function empirical_frequency --- src/utils.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/utils.jl b/src/utils.jl index 2bfb59ff..6d58743b 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,4 +1,5 @@ using Flux +using Statistics """ get_loss_fun(likelihood::Symbol) @@ -39,3 +40,27 @@ corresponding to the number of neurons on the last layer of the NN. function outdim(model::Chain)::Number return [size(p) for p in Flux.params(model)][end][1] end + +""" + count_fractions(Y_val,array) + +Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the empirical frequency +phat_j = {y_t|F_t(y_t)<= p_j, t= 1,....,T}/T, where T is the number of calibration points, p_j is the confidence level and F_t is the +cumulative density function of the predicted distribution targeting y_t. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 + Arguments: + Y_val: a vector of values + array: an array of sampled distributions stacked column-wise. +""" +function empirical_frequency(Y_cal,array) + + + quantiles= collect(0:0.05:1) + quantiles_matrix = hcat([quantile(samples, quantiles) for samples in array]...) + n_rows = size(bounds_quantiles_matrix,1) + counts = [] + + for i in 1:n_rows + push!(counts, sum(Y_cal .<= quantiles_matrix[i, :]) / length(Y_cal)) + end + return counts +end \ No newline at end of file From 25ea642fc3788fa2f58d09bebd85ff297965f077 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sun, 9 Jun 2024 05:59:02 +0200 Subject: [PATCH 02/46] fixed the docstring. --- src/utils.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index 6d58743b..d5e6caa3 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -42,8 +42,9 @@ function outdim(model::Chain)::Number end """ - count_fractions(Y_val,array) + empirical_frequency(Y_val,array) +FOR REGRESSION MODELS. Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the empirical frequency phat_j = {y_t|F_t(y_t)<= p_j, t= 1,....,T}/T, where T is the number of calibration points, p_j is the confidence level and F_t is the cumulative density function of the predicted distribution targeting y_t. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 From c290ed84e3acd403bf00b202d3bdc9977027fc73 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Fri, 14 Jun 2024 23:18:55 +0200 Subject: [PATCH 03/46] added sharpness and binary classification. i have yet to test them properly due to the issue with mljflux. --- src/utils.jl | 81 ++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 75 insertions(+), 6 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index d5e6caa3..684d4ad4 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -47,16 +47,16 @@ end FOR REGRESSION MODELS. Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the empirical frequency phat_j = {y_t|F_t(y_t)<= p_j, t= 1,....,T}/T, where T is the number of calibration points, p_j is the confidence level and F_t is the -cumulative density function of the predicted distribution targeting y_t. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 +cumulative distribution function of the predicted distribution targeting y_t. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 Arguments: - Y_val: a vector of values - array: an array of sampled distributions stacked column-wise. + Y_val: a vector of values y_t + array: an array of sampled distributions F(x_t) stacked column-wise. """ -function empirical_frequency(Y_cal,array) +function empirical_frequency(Y_cal,sampled_distributions) quantiles= collect(0:0.05:1) - quantiles_matrix = hcat([quantile(samples, quantiles) for samples in array]...) + quantiles_matrix = hcat([quantile(samples, quantiles) for samples in sampled_distributions]...) n_rows = size(bounds_quantiles_matrix,1) counts = [] @@ -64,4 +64,73 @@ function empirical_frequency(Y_cal,array) push!(counts, sum(Y_cal .<= quantiles_matrix[i, :]) / length(Y_cal)) end return counts -end \ No newline at end of file +end + +""" + sharpness(Y_val,array) + +FOR REGRESSION MODELS. +Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the +sharpness of the predicted distributions, i.e., the average of the variances var(F_t) predicted by the forecaster for each x_t +The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 + Arguments: + sampled_distributions: an array of sampled distributions F(x_t) stacked column-wise. +""" +function sharpness(sampled_distributions) + sharpness= mean(var.(sampled_distributions)) + return sharpness + +end + + + + +""" + empirical_frequency-classification(Y_val,array) + +FOR BINARY CLASSIFICATION MODELS. +Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} let p_t= H(x_t)∈[0,1] be the forecasted probability. +We group the p_t into intervals I-j for j= 1,2,...,m that form a partition of [0,1]. The function computes +the observed average p_j= T^-1_j ∑_{t:p_t ∈ I_j} y_j in each interval I_j. +The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 + Arguments: + Y_val: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. + sampled_distributions: an array of sampled distributions stacked column-wise where in the first row + there is the probability for the target class y_1=1 and in the second row y_0=0. +""" +function empirical_frequency_binary_classification(Y_cal,sampled_distributions) + + #unique_elements = unique(Y_cal) + # Create the mapping + #mapping = Dict(unique_elements[1] => 0, unique_elements[2] => 1) + + # Convert categorical data to numeric data + #numeric_array = [mapping[c] for c in categorical_array] + + # Create bins + num_bins=10 + bins = range(0, stop=1, length=num_bins+1) + + # Initialize arrays to hold predicted and empirical averages + pred_avg = zeros(num_bins) + emp_avg = zeros(num_bins) + total_pj_per_intervalj = zeros(num_bins) + + class_indices = (Y_cal .== 1) + + class_probs = sampled_distributions[1, :] + + for j in 0:0.1:0.9 + + push!(total_pj_per_intervalj,sum( j Date: Sat, 15 Jun 2024 07:57:58 +0200 Subject: [PATCH 04/46] added trapz to the list of dependencies. --- Project.toml | 1 + docs/Project.toml | 1 + 2 files changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index 24355ed3..0c3b1a38 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,7 @@ ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1" Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" diff --git a/docs/Project.toml b/docs/Project.toml index 49942c85..03062895 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -9,5 +9,6 @@ RDatasets = "ce6b1742-4840-55fa-b093-852dadbb1d8b" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TaijaPlotting = "bd7198b4-c7d6-400c-9bab-9a24614b0240" +Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" cuDNN = "02a925ec-e4fe-4b08-9a7e-0d78e3d38ccd" From 6a222103f783cd5f781ad78c382d8b429cde23a8 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sat, 15 Jun 2024 08:12:35 +0200 Subject: [PATCH 05/46] added Distributions to theproject --- Project.toml | 1 + docs/Project.toml | 1 + src/baselaplace/predicting.jl | 9 ++++++--- 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 0c3b1a38..535622d3 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.2.1" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" ComputationalResources = "ed09eef8-17a6-5b46-8889-db040fac31e3" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MLJFlux = "094fc8d1-fd35-5302-93ea-dabda2abf845" diff --git a/docs/Project.toml b/docs/Project.toml index 03062895..fedf9ac8 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,5 @@ [deps] +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" LaplaceRedux = "c52c1a26-f7c5-402b-80be-ba1e638ad478" diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index 49773a23..e69f299e 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -39,7 +39,9 @@ function glm_predictive_distribution(la::AbstractLaplace, X::AbstractArray) fμ = reshape(fμ, Flux.outputsize(la.model, size(X))) fvar = functional_variance(la, 𝐉) fvar = reshape(fvar, size(fμ)...) - return fμ, fvar + fstd = sqrt.(fvar) + normal_distr= [Distributions.Normal(fμ[i, j], fstd[i, j]) for i in 1:size(fμ, 1), j in 1:size(fμ, 2)] + return normal_distr end """ @@ -75,11 +77,12 @@ predict(la, hcat(x...)) function predict( la::AbstractLaplace, X::AbstractArray; link_approx=:probit, predict_proba::Bool=true ) - fμ, fvar = glm_predictive_distribution(la, X) + normal_distr = glm_predictive_distribution(la, X) + fμ, fvar = mean.(normal_distr), var.(normal_distr) # Regression: if la.likelihood == :regression - return fμ, fvar + return normal_distr end # Classification: From df3d60d50ff05052f112048edf92afd81b41d305 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sat, 15 Jun 2024 21:50:08 +0200 Subject: [PATCH 06/46] working version --- src/baselaplace/predicting.jl | 1 + src/utils.jl | 50 ++++++++++++++++++----------------- 2 files changed, 27 insertions(+), 24 deletions(-) diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index e69f299e..fcde981d 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -1,3 +1,4 @@ +using Distributions """ functional_variance(la::AbstractLaplace, 𝐉::AbstractArray) diff --git a/src/utils.jl b/src/utils.jl index 684d4ad4..f35f9263 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -94,41 +94,43 @@ We group the p_t into intervals I-j for j= 1,2,...,m that form a partition of [0 the observed average p_j= T^-1_j ∑_{t:p_t ∈ I_j} y_j in each interval I_j. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 Arguments: - Y_val: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. + y_binary: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. sampled_distributions: an array of sampled distributions stacked column-wise where in the first row there is the probability for the target class y_1=1 and in the second row y_0=0. """ -function empirical_frequency_binary_classification(Y_cal,sampled_distributions) +function empirical_frequency_binary_classification(y_binary,sampled_distributions) - #unique_elements = unique(Y_cal) - # Create the mapping - #mapping = Dict(unique_elements[1] => 0, unique_elements[2] => 1) - - # Convert categorical data to numeric data - #numeric_array = [mapping[c] for c in categorical_array] - - # Create bins - num_bins=10 - bins = range(0, stop=1, length=num_bins+1) - - # Initialize arrays to hold predicted and empirical averages - pred_avg = zeros(num_bins) - emp_avg = zeros(num_bins) - total_pj_per_intervalj = zeros(num_bins) + pred_avg= collect(range(0,step=0.1,stop=0.9)) + emp_avg = [] + total_pj_per_intervalj = [] + class_probs = sampled_distributions[1, :] - class_indices = (Y_cal .== 1) + for j in 1:10 + j_float = j / 10.0 -0.1 + push!(total_pj_per_intervalj,sum( j_float. j_float < x Date: Sat, 15 Jun 2024 22:36:31 +0200 Subject: [PATCH 07/46] ops forgot to add sharpness for the classification case --- src/utils.jl | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index f35f9263..883db8fb 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -67,7 +67,7 @@ function empirical_frequency(Y_cal,sampled_distributions) end """ - sharpness(Y_val,array) + sharpness(array) FOR REGRESSION MODELS. Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the @@ -95,6 +95,7 @@ the observed average p_j= T^-1_j ∑_{t:p_t ∈ I_j} y_j in each interval I_j. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 Arguments: y_binary: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. + sampled_distributions: an array of sampled distributions stacked column-wise where in the first row there is the probability for the target class y_1=1 and in the second row y_0=0. """ @@ -126,13 +127,31 @@ function empirical_frequency_binary_classification(y_binary,sampled_distribution end + return (total_pj_per_intervalj,emp_avg,pred_avg) +end - return (total_pj_per_intervalj,emp_avg,pred_avg) +""" + sharpness-classification(array) +FOR BINARY CLASSIFICATION MODELS. +We can also assess sharpness by looking at the distribution of model predictions.When forecasts are sharp, +most predictions are close to 0 or 1; unsharp forecasters make predictions closer to 0.5. +The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 -end + Arguments: + y_binary: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. + sampled_distributions: an array of sampled distributions stacked column-wise where in the first row + there is the probability for the target class y_1=1 and in the second row y_0=0. +""" +function sharpness_classification(y_binary,sampled_distributions) + class_one = sampled_distributions[1,findall(y_binary .== 1)] + class_zero = sampled_distributions[2,findall(y_binary .== 0)] + + return mean(class_one),mean(class_zero) + +end \ No newline at end of file From 07b318f482f0af8ca38eabd0408ba377eae74125 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Fri, 21 Jun 2024 16:38:06 +0200 Subject: [PATCH 08/46] working release.. changed changelog, glm_predictive_distribution, predict --- CHANGELOG.md | 16 ++++++++++++++++ Project.toml | 2 +- src/baselaplace/predicting.jl | 8 +++++--- 3 files changed, 22 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 916da018..6ddcbae2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,22 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), *Note*: We try to adhere to these practices as of version [v0.2.1]. + +## Version [0.3.0] - 2024-06-21 + +### Changed + +- Changed `glm_predictive_distribution` so that return a Normal distribution rather than the tuple (mean,variance). [#90] +- Changed `predict` so that return directly a Normal distribution in the case of regression. [#90] + +### Added + +- Added functions to compute the average empirical frequency for both classification and regression problems in utils.jl. [#90] + + + + + ## Version [0.2.1] - 2024-05-29 ### Changed diff --git a/Project.toml b/Project.toml index 535622d3..270fea5f 100644 --- a/Project.toml +++ b/Project.toml @@ -17,7 +17,6 @@ ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" -Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1" Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" @@ -26,6 +25,7 @@ Aqua = "0.8" ChainRulesCore = "1.23.0" Compat = "4.7.0" ComputationalResources = "0.3.2" +Distributions = "0.25.109" Flux = "0.12, 0.13, 0.14" LinearAlgebra = "1.6, 1.7, 1.8, 1.9, 1.10" MLJFlux = "0.2.10, 0.3, 0.4" diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index fcde981d..e58a5bc0 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -1,4 +1,5 @@ -using Distributions +using Distributions: Distributions +using Statistics: mean, var """ functional_variance(la::AbstractLaplace, 𝐉::AbstractArray) @@ -41,7 +42,8 @@ function glm_predictive_distribution(la::AbstractLaplace, X::AbstractArray) fvar = functional_variance(la, 𝐉) fvar = reshape(fvar, size(fμ)...) fstd = sqrt.(fvar) - normal_distr= [Distributions.Normal(fμ[i, j], fstd[i, j]) for i in 1:size(fμ, 1), j in 1:size(fμ, 2)] + normal_distr = [ + Distributions.Normal(fμ[i], fstd[i]) for i in 1:size(fμ, 1)] return normal_distr end @@ -79,7 +81,7 @@ function predict( la::AbstractLaplace, X::AbstractArray; link_approx=:probit, predict_proba::Bool=true ) normal_distr = glm_predictive_distribution(la, X) - fμ, fvar = mean.(normal_distr), var.(normal_distr) + fμ, fvar = mean.(normal_distr), var.(normal_distr) # Regression: if la.likelihood == :regression From eafa7bd52968e4c593815795697ef42f7feb266b Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sun, 9 Jun 2024 05:49:42 +0200 Subject: [PATCH 09/46] function empirical_frequency --- src/utils.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/src/utils.jl b/src/utils.jl index 2bfb59ff..6d58743b 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,4 +1,5 @@ using Flux +using Statistics """ get_loss_fun(likelihood::Symbol) @@ -39,3 +40,27 @@ corresponding to the number of neurons on the last layer of the NN. function outdim(model::Chain)::Number return [size(p) for p in Flux.params(model)][end][1] end + +""" + count_fractions(Y_val,array) + +Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the empirical frequency +phat_j = {y_t|F_t(y_t)<= p_j, t= 1,....,T}/T, where T is the number of calibration points, p_j is the confidence level and F_t is the +cumulative density function of the predicted distribution targeting y_t. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 + Arguments: + Y_val: a vector of values + array: an array of sampled distributions stacked column-wise. +""" +function empirical_frequency(Y_cal,array) + + + quantiles= collect(0:0.05:1) + quantiles_matrix = hcat([quantile(samples, quantiles) for samples in array]...) + n_rows = size(bounds_quantiles_matrix,1) + counts = [] + + for i in 1:n_rows + push!(counts, sum(Y_cal .<= quantiles_matrix[i, :]) / length(Y_cal)) + end + return counts +end \ No newline at end of file From f66e08ed6ce52b74e44bfbd7db7c27eb2e05637d Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sun, 9 Jun 2024 05:59:02 +0200 Subject: [PATCH 10/46] fixed the docstring. --- src/utils.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index 6d58743b..d5e6caa3 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -42,8 +42,9 @@ function outdim(model::Chain)::Number end """ - count_fractions(Y_val,array) + empirical_frequency(Y_val,array) +FOR REGRESSION MODELS. Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the empirical frequency phat_j = {y_t|F_t(y_t)<= p_j, t= 1,....,T}/T, where T is the number of calibration points, p_j is the confidence level and F_t is the cumulative density function of the predicted distribution targeting y_t. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 From 53552814226af47d68b14e7c11af1911befe3755 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Fri, 14 Jun 2024 23:18:55 +0200 Subject: [PATCH 11/46] added sharpness and binary classification. i have yet to test them properly due to the issue with mljflux. --- src/utils.jl | 81 ++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 75 insertions(+), 6 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index d5e6caa3..684d4ad4 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -47,16 +47,16 @@ end FOR REGRESSION MODELS. Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the empirical frequency phat_j = {y_t|F_t(y_t)<= p_j, t= 1,....,T}/T, where T is the number of calibration points, p_j is the confidence level and F_t is the -cumulative density function of the predicted distribution targeting y_t. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 +cumulative distribution function of the predicted distribution targeting y_t. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 Arguments: - Y_val: a vector of values - array: an array of sampled distributions stacked column-wise. + Y_val: a vector of values y_t + array: an array of sampled distributions F(x_t) stacked column-wise. """ -function empirical_frequency(Y_cal,array) +function empirical_frequency(Y_cal,sampled_distributions) quantiles= collect(0:0.05:1) - quantiles_matrix = hcat([quantile(samples, quantiles) for samples in array]...) + quantiles_matrix = hcat([quantile(samples, quantiles) for samples in sampled_distributions]...) n_rows = size(bounds_quantiles_matrix,1) counts = [] @@ -64,4 +64,73 @@ function empirical_frequency(Y_cal,array) push!(counts, sum(Y_cal .<= quantiles_matrix[i, :]) / length(Y_cal)) end return counts -end \ No newline at end of file +end + +""" + sharpness(Y_val,array) + +FOR REGRESSION MODELS. +Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the +sharpness of the predicted distributions, i.e., the average of the variances var(F_t) predicted by the forecaster for each x_t +The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 + Arguments: + sampled_distributions: an array of sampled distributions F(x_t) stacked column-wise. +""" +function sharpness(sampled_distributions) + sharpness= mean(var.(sampled_distributions)) + return sharpness + +end + + + + +""" + empirical_frequency-classification(Y_val,array) + +FOR BINARY CLASSIFICATION MODELS. +Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} let p_t= H(x_t)∈[0,1] be the forecasted probability. +We group the p_t into intervals I-j for j= 1,2,...,m that form a partition of [0,1]. The function computes +the observed average p_j= T^-1_j ∑_{t:p_t ∈ I_j} y_j in each interval I_j. +The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 + Arguments: + Y_val: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. + sampled_distributions: an array of sampled distributions stacked column-wise where in the first row + there is the probability for the target class y_1=1 and in the second row y_0=0. +""" +function empirical_frequency_binary_classification(Y_cal,sampled_distributions) + + #unique_elements = unique(Y_cal) + # Create the mapping + #mapping = Dict(unique_elements[1] => 0, unique_elements[2] => 1) + + # Convert categorical data to numeric data + #numeric_array = [mapping[c] for c in categorical_array] + + # Create bins + num_bins=10 + bins = range(0, stop=1, length=num_bins+1) + + # Initialize arrays to hold predicted and empirical averages + pred_avg = zeros(num_bins) + emp_avg = zeros(num_bins) + total_pj_per_intervalj = zeros(num_bins) + + class_indices = (Y_cal .== 1) + + class_probs = sampled_distributions[1, :] + + for j in 0:0.1:0.9 + + push!(total_pj_per_intervalj,sum( j Date: Sat, 15 Jun 2024 07:57:58 +0200 Subject: [PATCH 12/46] added trapz to the list of dependencies. --- Project.toml | 1 + docs/Project.toml | 1 + 2 files changed, 2 insertions(+) diff --git a/Project.toml b/Project.toml index 24355ed3..0c3b1a38 100644 --- a/Project.toml +++ b/Project.toml @@ -16,6 +16,7 @@ ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" +Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1" Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" diff --git a/docs/Project.toml b/docs/Project.toml index 49942c85..03062895 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -9,5 +9,6 @@ RDatasets = "ce6b1742-4840-55fa-b093-852dadbb1d8b" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" TaijaPlotting = "bd7198b4-c7d6-400c-9bab-9a24614b0240" +Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" cuDNN = "02a925ec-e4fe-4b08-9a7e-0d78e3d38ccd" From 26643ee25f5d93d06f8ef751135a1c01dd688dbd Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sat, 15 Jun 2024 08:12:35 +0200 Subject: [PATCH 13/46] added Distributions to theproject --- Project.toml | 1 + docs/Project.toml | 1 + src/baselaplace/predicting.jl | 9 ++++++--- 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 0c3b1a38..535622d3 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "0.2.1" ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" Compat = "34da2185-b29b-5c13-b0c7-acf172513d20" ComputationalResources = "ed09eef8-17a6-5b46-8889-db040fac31e3" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MLJFlux = "094fc8d1-fd35-5302-93ea-dabda2abf845" diff --git a/docs/Project.toml b/docs/Project.toml index 03062895..fedf9ac8 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,4 +1,5 @@ [deps] +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" LaplaceRedux = "c52c1a26-f7c5-402b-80be-ba1e638ad478" diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index 49773a23..e69f299e 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -39,7 +39,9 @@ function glm_predictive_distribution(la::AbstractLaplace, X::AbstractArray) fμ = reshape(fμ, Flux.outputsize(la.model, size(X))) fvar = functional_variance(la, 𝐉) fvar = reshape(fvar, size(fμ)...) - return fμ, fvar + fstd = sqrt.(fvar) + normal_distr= [Distributions.Normal(fμ[i, j], fstd[i, j]) for i in 1:size(fμ, 1), j in 1:size(fμ, 2)] + return normal_distr end """ @@ -75,11 +77,12 @@ predict(la, hcat(x...)) function predict( la::AbstractLaplace, X::AbstractArray; link_approx=:probit, predict_proba::Bool=true ) - fμ, fvar = glm_predictive_distribution(la, X) + normal_distr = glm_predictive_distribution(la, X) + fμ, fvar = mean.(normal_distr), var.(normal_distr) # Regression: if la.likelihood == :regression - return fμ, fvar + return normal_distr end # Classification: From b79ca39273a9c11ab41656d4f208cbd38e172136 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sat, 15 Jun 2024 21:50:08 +0200 Subject: [PATCH 14/46] working version --- src/baselaplace/predicting.jl | 1 + src/utils.jl | 50 ++++++++++++++++++----------------- 2 files changed, 27 insertions(+), 24 deletions(-) diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index e69f299e..fcde981d 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -1,3 +1,4 @@ +using Distributions """ functional_variance(la::AbstractLaplace, 𝐉::AbstractArray) diff --git a/src/utils.jl b/src/utils.jl index 684d4ad4..f35f9263 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -94,41 +94,43 @@ We group the p_t into intervals I-j for j= 1,2,...,m that form a partition of [0 the observed average p_j= T^-1_j ∑_{t:p_t ∈ I_j} y_j in each interval I_j. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 Arguments: - Y_val: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. + y_binary: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. sampled_distributions: an array of sampled distributions stacked column-wise where in the first row there is the probability for the target class y_1=1 and in the second row y_0=0. """ -function empirical_frequency_binary_classification(Y_cal,sampled_distributions) +function empirical_frequency_binary_classification(y_binary,sampled_distributions) - #unique_elements = unique(Y_cal) - # Create the mapping - #mapping = Dict(unique_elements[1] => 0, unique_elements[2] => 1) - - # Convert categorical data to numeric data - #numeric_array = [mapping[c] for c in categorical_array] - - # Create bins - num_bins=10 - bins = range(0, stop=1, length=num_bins+1) - - # Initialize arrays to hold predicted and empirical averages - pred_avg = zeros(num_bins) - emp_avg = zeros(num_bins) - total_pj_per_intervalj = zeros(num_bins) + pred_avg= collect(range(0,step=0.1,stop=0.9)) + emp_avg = [] + total_pj_per_intervalj = [] + class_probs = sampled_distributions[1, :] - class_indices = (Y_cal .== 1) + for j in 1:10 + j_float = j / 10.0 -0.1 + push!(total_pj_per_intervalj,sum( j_float. j_float < x Date: Sat, 15 Jun 2024 22:36:31 +0200 Subject: [PATCH 15/46] ops forgot to add sharpness for the classification case --- src/utils.jl | 25 ++++++++++++++++++++++--- 1 file changed, 22 insertions(+), 3 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index f35f9263..883db8fb 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -67,7 +67,7 @@ function empirical_frequency(Y_cal,sampled_distributions) end """ - sharpness(Y_val,array) + sharpness(array) FOR REGRESSION MODELS. Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the @@ -95,6 +95,7 @@ the observed average p_j= T^-1_j ∑_{t:p_t ∈ I_j} y_j in each interval I_j. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 Arguments: y_binary: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. + sampled_distributions: an array of sampled distributions stacked column-wise where in the first row there is the probability for the target class y_1=1 and in the second row y_0=0. """ @@ -126,13 +127,31 @@ function empirical_frequency_binary_classification(y_binary,sampled_distribution end + return (total_pj_per_intervalj,emp_avg,pred_avg) +end - return (total_pj_per_intervalj,emp_avg,pred_avg) +""" + sharpness-classification(array) +FOR BINARY CLASSIFICATION MODELS. +We can also assess sharpness by looking at the distribution of model predictions.When forecasts are sharp, +most predictions are close to 0 or 1; unsharp forecasters make predictions closer to 0.5. +The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 -end + Arguments: + y_binary: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. + sampled_distributions: an array of sampled distributions stacked column-wise where in the first row + there is the probability for the target class y_1=1 and in the second row y_0=0. +""" +function sharpness_classification(y_binary,sampled_distributions) + class_one = sampled_distributions[1,findall(y_binary .== 1)] + class_zero = sampled_distributions[2,findall(y_binary .== 0)] + + return mean(class_one),mean(class_zero) + +end \ No newline at end of file From 5f772cf7e0fe9299d3e0453f3a23fbcbaa01bc7c Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Fri, 21 Jun 2024 16:38:06 +0200 Subject: [PATCH 16/46] working release.. changed changelog, glm_predictive_distribution, predict --- CHANGELOG.md | 16 ++++++++++++++++ Project.toml | 2 +- src/baselaplace/predicting.jl | 8 +++++--- 3 files changed, 22 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 916da018..6ddcbae2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,22 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), *Note*: We try to adhere to these practices as of version [v0.2.1]. + +## Version [0.3.0] - 2024-06-21 + +### Changed + +- Changed `glm_predictive_distribution` so that return a Normal distribution rather than the tuple (mean,variance). [#90] +- Changed `predict` so that return directly a Normal distribution in the case of regression. [#90] + +### Added + +- Added functions to compute the average empirical frequency for both classification and regression problems in utils.jl. [#90] + + + + + ## Version [0.2.1] - 2024-05-29 ### Changed diff --git a/Project.toml b/Project.toml index 535622d3..270fea5f 100644 --- a/Project.toml +++ b/Project.toml @@ -17,7 +17,6 @@ ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" -Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1" Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" @@ -26,6 +25,7 @@ Aqua = "0.8" ChainRulesCore = "1.23.0" Compat = "4.7.0" ComputationalResources = "0.3.2" +Distributions = "0.25.109" Flux = "0.12, 0.13, 0.14" LinearAlgebra = "1.6, 1.7, 1.8, 1.9, 1.10" MLJFlux = "0.2.10, 0.3, 0.4" diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index fcde981d..e58a5bc0 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -1,4 +1,5 @@ -using Distributions +using Distributions: Distributions +using Statistics: mean, var """ functional_variance(la::AbstractLaplace, 𝐉::AbstractArray) @@ -41,7 +42,8 @@ function glm_predictive_distribution(la::AbstractLaplace, X::AbstractArray) fvar = functional_variance(la, 𝐉) fvar = reshape(fvar, size(fμ)...) fstd = sqrt.(fvar) - normal_distr= [Distributions.Normal(fμ[i, j], fstd[i, j]) for i in 1:size(fμ, 1), j in 1:size(fμ, 2)] + normal_distr = [ + Distributions.Normal(fμ[i], fstd[i]) for i in 1:size(fμ, 1)] return normal_distr end @@ -79,7 +81,7 @@ function predict( la::AbstractLaplace, X::AbstractArray; link_approx=:probit, predict_proba::Bool=true ) normal_distr = glm_predictive_distribution(la, X) - fμ, fvar = mean.(normal_distr), var.(normal_distr) + fμ, fvar = mean.(normal_distr), var.(normal_distr) # Regression: if la.likelihood == :regression From 7af9378017bb2373eed76478dac637debfe52135 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Fri, 21 Jun 2024 17:11:47 +0200 Subject: [PATCH 17/46] changed docstrings in predicting.jl --- src/baselaplace/predicting.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index e58a5bc0..c40f657f 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -24,6 +24,8 @@ Computes the linearized GLM predictive. - `fμ::AbstractArray`: Mean of the predictive distribution. The output shape is column-major as in Flux. - `fvar::AbstractArray`: Variance of the predictive distribution. The output shape is column-major as in Flux. +- `normal_distr` An array of normal distributions approximating the predictive distribution p(y|X) given the input data X. + # Examples ```julia-repl @@ -60,9 +62,10 @@ Computes predictions from Bayesian neural network. - `predict_proba::Bool=true`: If `true` (default), returns probabilities for classification tasks. # Returns - +For classification tasks: - `fμ::AbstractArray`: Mean of the predictive distribution if link function is set to `:plugin`, otherwise the probit approximation. The output shape is column-major as in Flux. -- `fvar::AbstractArray`: If regression, it also returns the variance of the predictive distribution. The output shape is column-major as in Flux. +For regression tasks: +- `normal_distr::Distributions.Normal`:the array of Normal distributions computed by glm_predictive_distribution. The output shape is column-major as in Flux. # Examples From 2c422362456e727c09dcc3354168e9fd620a98c9 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sat, 22 Jun 2024 10:46:04 +0200 Subject: [PATCH 18/46] fixed glm_predictive_distribution --- CHANGELOG.md | 5 +++++ src/baselaplace/predicting.jl | 6 +++--- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6ddcbae2..24cb8d30 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,11 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), *Note*: We try to adhere to these practices as of version [v0.2.1]. +## Version [0.3.1] - 2024-06-22 + +### Changed + +- Changed `glm_predictive_distribution` so that return a tuple(Normal distribution,fμ, fvar) rather than the tuple (mean,variance). [#90] ## Version [0.3.0] - 2024-06-21 diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index c40f657f..5d573339 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -46,7 +46,7 @@ function glm_predictive_distribution(la::AbstractLaplace, X::AbstractArray) fstd = sqrt.(fvar) normal_distr = [ Distributions.Normal(fμ[i], fstd[i]) for i in 1:size(fμ, 1)] - return normal_distr + return (normal_distr,fμ,fvar ) end """ @@ -83,8 +83,8 @@ predict(la, hcat(x...)) function predict( la::AbstractLaplace, X::AbstractArray; link_approx=:probit, predict_proba::Bool=true ) - normal_distr = glm_predictive_distribution(la, X) - fμ, fvar = mean.(normal_distr), var.(normal_distr) + normal_distr,fμ, fvar = glm_predictive_distribution(la, X) + #fμ, fvar = mean.(normal_distr), var.(normal_distr) # Regression: if la.likelihood == :regression From 9d67ddc50201d813143054a1dfc0f4a390469407 Mon Sep 17 00:00:00 2001 From: pasquale c <343guiltyspark@outlook.it> Date: Sat, 22 Jun 2024 11:41:28 +0200 Subject: [PATCH 19/46] Update src/utils.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/utils.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 883db8fb..7e888393 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -147,11 +147,9 @@ The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.0026 sampled_distributions: an array of sampled distributions stacked column-wise where in the first row there is the probability for the target class y_1=1 and in the second row y_0=0. """ -function sharpness_classification(y_binary,sampled_distributions) +function sharpness_classification(y_binary, sampled_distributions) + class_one = sampled_distributions[1, findall(y_binary .== 1)] + class_zero = sampled_distributions[2, findall(y_binary .== 0)] - class_one = sampled_distributions[1,findall(y_binary .== 1)] - class_zero = sampled_distributions[2,findall(y_binary .== 0)] - - return mean(class_one),mean(class_zero) - + return mean(class_one), mean(class_zero) end \ No newline at end of file From 9f075837ad3624d1ae94a4d99d5de4aae4dcfb3c Mon Sep 17 00:00:00 2001 From: pasquale c <343guiltyspark@outlook.it> Date: Sat, 22 Jun 2024 11:41:53 +0200 Subject: [PATCH 20/46] Update src/utils.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/utils.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 7e888393..c41400c9 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -131,9 +131,6 @@ function empirical_frequency_binary_classification(y_binary,sampled_distribution end - - - """ sharpness-classification(array) From f81d226e3a570dcf89a6cc4f8f0c9f5b83de2a11 Mon Sep 17 00:00:00 2001 From: pasquale c <343guiltyspark@outlook.it> Date: Sat, 22 Jun 2024 11:43:19 +0200 Subject: [PATCH 21/46] Update src/utils.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/utils.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index c41400c9..43a3e2ad 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -127,8 +127,7 @@ function empirical_frequency_binary_classification(y_binary,sampled_distribution end - return (total_pj_per_intervalj,emp_avg,pred_avg) - + return (total_pj_per_intervalj, emp_avg, pred_avg) end """ From 6cdc50313d724acf49159ff97b63a006f1d002ad Mon Sep 17 00:00:00 2001 From: pasquale c <343guiltyspark@outlook.it> Date: Sat, 22 Jun 2024 11:44:33 +0200 Subject: [PATCH 22/46] Update src/baselaplace/predicting.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/baselaplace/predicting.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index 5d573339..c6fc4ce0 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -44,9 +44,8 @@ function glm_predictive_distribution(la::AbstractLaplace, X::AbstractArray) fvar = functional_variance(la, 𝐉) fvar = reshape(fvar, size(fμ)...) fstd = sqrt.(fvar) - normal_distr = [ - Distributions.Normal(fμ[i], fstd[i]) for i in 1:size(fμ, 1)] - return (normal_distr,fμ,fvar ) + normal_distr = [Distributions.Normal(fμ[i], fstd[i]) for i in 1:size(fμ, 1)] + return (normal_distr, fμ, fvar) end """ From 89bb19b421b6059423735613301f21efdd6bc0ab Mon Sep 17 00:00:00 2001 From: pasquale c <343guiltyspark@outlook.it> Date: Sat, 22 Jun 2024 11:44:48 +0200 Subject: [PATCH 23/46] Update src/baselaplace/predicting.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/baselaplace/predicting.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index c6fc4ce0..d0c66ebd 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -82,7 +82,7 @@ predict(la, hcat(x...)) function predict( la::AbstractLaplace, X::AbstractArray; link_approx=:probit, predict_proba::Bool=true ) - normal_distr,fμ, fvar = glm_predictive_distribution(la, X) + normal_distr, fμ, fvar = glm_predictive_distribution(la, X) #fμ, fvar = mean.(normal_distr), var.(normal_distr) # Regression: From 6fe01a2169c1c268aa1e73ee926621a58d78669f Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sat, 22 Jun 2024 12:06:12 +0200 Subject: [PATCH 24/46] JuliaFormatter --- src/utils.jl | 49 ++++++++++++++++++------------------------------- 1 file changed, 18 insertions(+), 31 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 43a3e2ad..75444b19 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -52,15 +52,15 @@ cumulative distribution function of the predicted distribution targeting y_t. Th Y_val: a vector of values y_t array: an array of sampled distributions F(x_t) stacked column-wise. """ -function empirical_frequency(Y_cal,sampled_distributions) - - - quantiles= collect(0:0.05:1) - quantiles_matrix = hcat([quantile(samples, quantiles) for samples in sampled_distributions]...) - n_rows = size(bounds_quantiles_matrix,1) +function empirical_frequency(Y_cal, sampled_distributions) + quantiles = collect(0:0.05:1) + quantiles_matrix = hcat( + [quantile(samples, quantiles) for samples in sampled_distributions]... + ) + n_rows = size(bounds_quantiles_matrix, 1) counts = [] - for i in 1:n_rows + for i in 1:n_rows push!(counts, sum(Y_cal .<= quantiles_matrix[i, :]) / length(Y_cal)) end return counts @@ -77,14 +77,10 @@ The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.0026 sampled_distributions: an array of sampled distributions F(x_t) stacked column-wise. """ function sharpness(sampled_distributions) - sharpness= mean(var.(sampled_distributions)) + sharpness = mean(var.(sampled_distributions)) return sharpness - end - - - """ empirical_frequency-classification(Y_val,array) @@ -99,32 +95,24 @@ The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.0026 sampled_distributions: an array of sampled distributions stacked column-wise where in the first row there is the probability for the target class y_1=1 and in the second row y_0=0. """ -function empirical_frequency_binary_classification(y_binary,sampled_distributions) - - pred_avg= collect(range(0,step=0.1,stop=0.9)) +function empirical_frequency_binary_classification(y_binary, sampled_distributions) + pred_avg = collect(range(0; step=0.1, stop=0.9)) emp_avg = [] total_pj_per_intervalj = [] class_probs = sampled_distributions[1, :] for j in 1:10 - j_float = j / 10.0 -0.1 - push!(total_pj_per_intervalj,sum( j_float. j_float < x j_float < x < j_float + 0.1, class_probs) + push!(emp_avg, 1 / total_pj_per_intervalj[j] * sum(y_binary[indices])) println(" numero $j") - pred_avg[j] = 1/total_pj_per_intervalj[j] * sum(sampled_distributions[1,indices]) + pred_avg[j] = + 1 / total_pj_per_intervalj[j] * sum(sampled_distributions[1, indices]) end - end return (total_pj_per_intervalj, emp_avg, pred_avg) @@ -146,6 +134,5 @@ The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.0026 function sharpness_classification(y_binary, sampled_distributions) class_one = sampled_distributions[1, findall(y_binary .== 1)] class_zero = sampled_distributions[2, findall(y_binary .== 0)] - return mean(class_one), mean(class_zero) -end \ No newline at end of file +end From 0bba48856c381f4495cf238927927c9cc1f0197c Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sun, 23 Jun 2024 23:39:51 +0200 Subject: [PATCH 25/46] fixed docstrings --- src/utils.jl | 35 ++++++++++++++++++++--------------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 75444b19..0d781bc5 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -42,15 +42,17 @@ function outdim(model::Chain)::Number end """ - empirical_frequency(Y_val,array) + empirical_frequency(Y_cal, sampled_distributions) FOR REGRESSION MODELS. Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the empirical frequency phat_j = {y_t|F_t(y_t)<= p_j, t= 1,....,T}/T, where T is the number of calibration points, p_j is the confidence level and F_t is the -cumulative distribution function of the predicted distribution targeting y_t. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 +cumulative distribution function of the predicted distribution targeting y_t. + The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 + Arguments: - Y_val: a vector of values y_t - array: an array of sampled distributions F(x_t) stacked column-wise. + -Y_cal: a vector of values y_t + -sampled_distributions: an array of sampled distributions F(x_t) stacked column-wise. """ function empirical_frequency(Y_cal, sampled_distributions) quantiles = collect(0:0.05:1) @@ -67,14 +69,15 @@ function empirical_frequency(Y_cal, sampled_distributions) end """ - sharpness(array) + sharpness(sampled_distributions) FOR REGRESSION MODELS. Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the sharpness of the predicted distributions, i.e., the average of the variances var(F_t) predicted by the forecaster for each x_t The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 + Arguments: - sampled_distributions: an array of sampled distributions F(x_t) stacked column-wise. + -sampled_distributions: an array of sampled distributions F(x_t) stacked column-wise. """ function sharpness(sampled_distributions) sharpness = mean(var.(sampled_distributions)) @@ -82,18 +85,19 @@ function sharpness(sampled_distributions) end """ - empirical_frequency-classification(Y_val,array) + empirical_frequency-classification(y_binary, sampled_distributions) FOR BINARY CLASSIFICATION MODELS. Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} let p_t= H(x_t)∈[0,1] be the forecasted probability. We group the p_t into intervals I-j for j= 1,2,...,m that form a partition of [0,1]. The function computes the observed average p_j= T^-1_j ∑_{t:p_t ∈ I_j} y_j in each interval I_j. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 + Arguments: - y_binary: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. + -y_binary: the array of outputs y_t numerically coded: 1 for the target class, 0 for the negative result. - sampled_distributions: an array of sampled distributions stacked column-wise where in the first row - there is the probability for the target class y_1=1 and in the second row y_0=0. + -sampled_distributions: an array of sampled distributions stacked column-wise so that in the first row + there is the probability for the target class y_1 and in the second row the probability for the null class y_0. """ function empirical_frequency_binary_classification(y_binary, sampled_distributions) pred_avg = collect(range(0; step=0.1, stop=0.9)) @@ -119,17 +123,18 @@ function empirical_frequency_binary_classification(y_binary, sampled_distributio end """ - sharpness-classification(array) + sharpness-classification(y_binary,sampled_distributions) FOR BINARY CLASSIFICATION MODELS. We can also assess sharpness by looking at the distribution of model predictions.When forecasts are sharp, -most predictions are close to 0 or 1; unsharp forecasters make predictions closer to 0.5. +most predictions are close to 0 or 1; not sharp forecasters make predictions closer to 0.5. The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 Arguments: - y_binary: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. - sampled_distributions: an array of sampled distributions stacked column-wise where in the first row - there is the probability for the target class y_1=1 and in the second row y_0=0. + + -y_binary: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. + -sampled_distributions: an array of sampled distributions stacked column-wise so that in the first row + there is the probability for the target class y_1 and in the second row the probability for the null class y_0. """ function sharpness_classification(y_binary, sampled_distributions) class_one = sampled_distributions[1, findall(y_binary .== 1)] From 8311de34a70cc375eeb7b70637d01422c8f09b06 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sun, 23 Jun 2024 23:43:21 +0200 Subject: [PATCH 26/46] made docstrings a lil bit shorter --- src/utils.jl | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 0d781bc5..e3d7ae36 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -48,7 +48,7 @@ FOR REGRESSION MODELS. Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the empirical frequency phat_j = {y_t|F_t(y_t)<= p_j, t= 1,....,T}/T, where T is the number of calibration points, p_j is the confidence level and F_t is the cumulative distribution function of the predicted distribution targeting y_t. - The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 +Source: https://arxiv.org/abs/1807.00263 Arguments: -Y_cal: a vector of values y_t @@ -74,7 +74,7 @@ end FOR REGRESSION MODELS. Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the sharpness of the predicted distributions, i.e., the average of the variances var(F_t) predicted by the forecaster for each x_t -The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 +Source: https://arxiv.org/abs/1807.00263 Arguments: -sampled_distributions: an array of sampled distributions F(x_t) stacked column-wise. @@ -91,7 +91,7 @@ FOR BINARY CLASSIFICATION MODELS. Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} let p_t= H(x_t)∈[0,1] be the forecasted probability. We group the p_t into intervals I-j for j= 1,2,...,m that form a partition of [0,1]. The function computes the observed average p_j= T^-1_j ∑_{t:p_t ∈ I_j} y_j in each interval I_j. -The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 +Source: https://arxiv.org/abs/1807.00263 Arguments: -y_binary: the array of outputs y_t numerically coded: 1 for the target class, 0 for the negative result. @@ -126,12 +126,11 @@ end sharpness-classification(y_binary,sampled_distributions) FOR BINARY CLASSIFICATION MODELS. -We can also assess sharpness by looking at the distribution of model predictions.When forecasts are sharp, +Assess the sharpness of the model by looking at the distribution of model predictions. When forecasts are sharp, most predictions are close to 0 or 1; not sharp forecasters make predictions closer to 0.5. -The function was suggested by Kuleshov(2018) in https://arxiv.org/abs/1807.00263 +Source: https://arxiv.org/abs/1807.00263 Arguments: - -y_binary: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. -sampled_distributions: an array of sampled distributions stacked column-wise so that in the first row there is the probability for the target class y_1 and in the second row the probability for the null class y_0. From 78373334f8e40a745630e94273da8846d6c5d439 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Mon, 24 Jun 2024 02:17:40 +0200 Subject: [PATCH 27/46] docstrings again (added output) --- src/utils.jl | 38 ++++++++++++++++++++++++-------------- 1 file changed, 24 insertions(+), 14 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index e3d7ae36..16dc236b 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -50,9 +50,11 @@ phat_j = {y_t|F_t(y_t)<= p_j, t= 1,....,T}/T, where T is the number of calibrati cumulative distribution function of the predicted distribution targeting y_t. Source: https://arxiv.org/abs/1807.00263 - Arguments: - -Y_cal: a vector of values y_t - -sampled_distributions: an array of sampled distributions F(x_t) stacked column-wise. +Inputs: + - Y_cal: a vector of values y_t + - sampled_distributions: an array of sampled distributions F(x_t) stacked column-wise. +Outputs: + - counts: an array cointaining the empirical frequencies for each quantile interval. """ function empirical_frequency(Y_cal, sampled_distributions) quantiles = collect(0:0.05:1) @@ -76,8 +78,10 @@ Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predi sharpness of the predicted distributions, i.e., the average of the variances var(F_t) predicted by the forecaster for each x_t Source: https://arxiv.org/abs/1807.00263 - Arguments: - -sampled_distributions: an array of sampled distributions F(x_t) stacked column-wise. +Inputs: + - sampled_distributions: an array of sampled distributions F(x_t) stacked column-wise. +Outputs: + - sharpness: a scalar that measure the level of sharpness of the regressor """ function sharpness(sampled_distributions) sharpness = mean(var.(sampled_distributions)) @@ -93,11 +97,14 @@ We group the p_t into intervals I-j for j= 1,2,...,m that form a partition of [0 the observed average p_j= T^-1_j ∑_{t:p_t ∈ I_j} y_j in each interval I_j. Source: https://arxiv.org/abs/1807.00263 - Arguments: - -y_binary: the array of outputs y_t numerically coded: 1 for the target class, 0 for the negative result. - - -sampled_distributions: an array of sampled distributions stacked column-wise so that in the first row +Inputs: + - y_binary: the array of outputs y_t numerically coded: 1 for the target class, 0 for the null class. + - sampled_distributions: an array of sampled distributions stacked column-wise so that in the first row there is the probability for the target class y_1 and in the second row the probability for the null class y_0. +Outputs: + - total_pj_per_intervalj: + - emp_avg: + - pred_avg: """ function empirical_frequency_binary_classification(y_binary, sampled_distributions) pred_avg = collect(range(0; step=0.1, stop=0.9)) @@ -130,13 +137,16 @@ Assess the sharpness of the model by looking at the distribution of model predi most predictions are close to 0 or 1; not sharp forecasters make predictions closer to 0.5. Source: https://arxiv.org/abs/1807.00263 - Arguments: +Inputs: -y_binary: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. -sampled_distributions: an array of sampled distributions stacked column-wise so that in the first row - there is the probability for the target class y_1 and in the second row the probability for the null class y_0. + there is the probability for the target class and in the second row the probability for the null class. + Outputs: + - mean_class_one: a scalar that measure the average prediction for the target class + - mean_class_zero: a scalar that measure the average prediction for the null class """ function sharpness_classification(y_binary, sampled_distributions) - class_one = sampled_distributions[1, findall(y_binary .== 1)] - class_zero = sampled_distributions[2, findall(y_binary .== 0)] - return mean(class_one), mean(class_zero) + mean_class_one = mean(sampled_distributions[1, findall(y_binary .== 1)]) + mean_class_zero= mean(sampled_distributions[2, findall(y_binary .== 0)]) + return mean_class_one, mean_class_zero end From b0518b20452296f8e0e9820f3d92295fd693c389 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Mon, 24 Jun 2024 10:14:18 +0200 Subject: [PATCH 28/46] fixed binary classification case, exported function from utils. --- src/LaplaceRedux.jl | 1 + src/utils.jl | 53 +++++++++++++++++++++++++++++---------------- 2 files changed, 35 insertions(+), 19 deletions(-) diff --git a/src/LaplaceRedux.jl b/src/LaplaceRedux.jl index 6ac5e95d..c87a1f4f 100644 --- a/src/LaplaceRedux.jl +++ b/src/LaplaceRedux.jl @@ -1,6 +1,7 @@ module LaplaceRedux include("utils.jl") +export empirical_frequency_binary_classification, sharpness_classification, empirical_frequency_regression, sharpness_regression include("data/Data.jl") using .Data diff --git a/src/utils.jl b/src/utils.jl index 16dc236b..4d1c3911 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -56,7 +56,7 @@ Inputs: Outputs: - counts: an array cointaining the empirical frequencies for each quantile interval. """ -function empirical_frequency(Y_cal, sampled_distributions) +function empirical_frequency_regression(Y_cal, sampled_distributions) quantiles = collect(0:0.05:1) quantiles_matrix = hcat( [quantile(samples, quantiles) for samples in sampled_distributions]... @@ -83,7 +83,7 @@ Inputs: Outputs: - sharpness: a scalar that measure the level of sharpness of the regressor """ -function sharpness(sampled_distributions) +function sharpness_regression(sampled_distributions) sharpness = mean(var.(sampled_distributions)) return sharpness end @@ -102,31 +102,46 @@ Inputs: - sampled_distributions: an array of sampled distributions stacked column-wise so that in the first row there is the probability for the target class y_1 and in the second row the probability for the null class y_0. Outputs: - - total_pj_per_intervalj: - - emp_avg: - - pred_avg: + - num_p_per_interval: array with the number of probabilities falling within interval + - emp_avg: array with the observed empirical average per interval + - bin_centers: array with the centers of the bins + """ function empirical_frequency_binary_classification(y_binary, sampled_distributions) - pred_avg = collect(range(0; step=0.1, stop=0.9)) + + # Number of bins + n_bins=20 + #intervals boundaries + int_bds = collect(range(0, stop=1, length=n_bins + 1)) + #bin centers + bin_centers = [(int_bds[i] + int_bds[i+1]) / 2 for i in 1:length(int_bds)-1] + #initialize list for empirical averages per interval emp_avg = [] - total_pj_per_intervalj = [] + #initialize list for predicted averages per interval + pred_avg=[] + # initialize list of number of probabilities falling within each intervals + num_p_per_interval = [] + #list of the predicted probabilities for the target class class_probs = sampled_distributions[1, :] - - for j in 1:10 - j_float = j / 10.0 - 0.1 - push!(total_pj_per_intervalj, sum(j_float .< class_probs .< j_float + 0.1)) - if total_pj_per_intervalj[j] == 0 + # iterate over the bins + for j in 1:n_bins + push!(num_p_per_interval, sum( int_bds[j]. j_float < x < j_float + 0.1, class_probs) - push!(emp_avg, 1 / total_pj_per_intervalj[j] * sum(y_binary[indices])) - println(" numero $j") - pred_avg[j] = - 1 / total_pj_per_intervalj[j] * sum(sampled_distributions[1, indices]) + # find the indices fo all istances for which class_probs fall withing the j-th interval + indices = findall(x -> int_bds[j] < x Date: Mon, 24 Jun 2024 10:22:30 +0200 Subject: [PATCH 29/46] juliaformatter --- src/LaplaceRedux.jl | 3 ++- src/utils.jl | 17 ++++++++--------- 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/src/LaplaceRedux.jl b/src/LaplaceRedux.jl index c87a1f4f..619a9d3b 100644 --- a/src/LaplaceRedux.jl +++ b/src/LaplaceRedux.jl @@ -1,7 +1,8 @@ module LaplaceRedux include("utils.jl") -export empirical_frequency_binary_classification, sharpness_classification, empirical_frequency_regression, sharpness_regression +export empirical_frequency_binary_classification, + sharpness_classification, empirical_frequency_regression, sharpness_regression include("data/Data.jl") using .Data diff --git a/src/utils.jl b/src/utils.jl index 4d1c3911..4935bceb 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -110,34 +110,33 @@ Outputs: function empirical_frequency_binary_classification(y_binary, sampled_distributions) # Number of bins - n_bins=20 + n_bins = 20 #intervals boundaries - int_bds = collect(range(0, stop=1, length=n_bins + 1)) + int_bds = collect(range(0; stop=1, length=n_bins + 1)) #bin centers - bin_centers = [(int_bds[i] + int_bds[i+1]) / 2 for i in 1:length(int_bds)-1] + bin_centers = [(int_bds[i] + int_bds[i + 1]) / 2 for i in 1:(length(int_bds) - 1)] #initialize list for empirical averages per interval emp_avg = [] #initialize list for predicted averages per interval - pred_avg=[] + pred_avg = [] # initialize list of number of probabilities falling within each intervals num_p_per_interval = [] #list of the predicted probabilities for the target class class_probs = sampled_distributions[1, :] # iterate over the bins for j in 1:n_bins - push!(num_p_per_interval, sum( int_bds[j]. int_bds[j] < x int_bds[j] < x < int_bds[j + 1], class_probs) #compute the empirical average and saved it in emp_avg in the j-th position push!(emp_avg, 1 / num_p_per_interval[j] * sum(y_binary[indices])) #TO DO: maybe substitute to bin_Centers? - push!(pred_avg, 1 / num_p_per_interval[j] * sum(class_probs[1,indices])) - + push!(pred_avg, 1 / num_p_per_interval[j] * sum(class_probs[1, indices])) end end #return the tuple @@ -162,6 +161,6 @@ Inputs: """ function sharpness_classification(y_binary, sampled_distributions) mean_class_one = mean(sampled_distributions[1, findall(y_binary .== 1)]) - mean_class_zero= mean(sampled_distributions[2, findall(y_binary .== 0)]) + mean_class_zero = mean(sampled_distributions[2, findall(y_binary .== 0)]) return mean_class_one, mean_class_zero end From 203513d29b5ca3a2c4b8e4b014d254c5f78acf2f Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sat, 29 Jun 2024 14:10:36 +0200 Subject: [PATCH 30/46] add n_bins as argument to functions --- src/utils.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 4935bceb..14627b8a 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -51,13 +51,14 @@ cumulative distribution function of the predicted distribution targeting y_t. Source: https://arxiv.org/abs/1807.00263 Inputs: - - Y_cal: a vector of values y_t - - sampled_distributions: an array of sampled distributions F(x_t) stacked column-wise. + - 'Y_cal': a vector of values y_t + - 'sampled_distributions': an array of sampled distributions F(x_t) stacked column-wise. + - 'n_bins': number of equally spaced bins to use. Outputs: - counts: an array cointaining the empirical frequencies for each quantile interval. """ -function empirical_frequency_regression(Y_cal, sampled_distributions) - quantiles = collect(0:0.05:1) +function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins) + quantiles = collect(range(0; stop=1, length=n_bins + 1)) quantiles_matrix = hcat( [quantile(samples, quantiles) for samples in sampled_distributions]... ) @@ -101,16 +102,17 @@ Inputs: - y_binary: the array of outputs y_t numerically coded: 1 for the target class, 0 for the null class. - sampled_distributions: an array of sampled distributions stacked column-wise so that in the first row there is the probability for the target class y_1 and in the second row the probability for the null class y_0. + - 'n_bins': number of equally spaced bins to use. Outputs: - num_p_per_interval: array with the number of probabilities falling within interval - emp_avg: array with the observed empirical average per interval - bin_centers: array with the centers of the bins """ -function empirical_frequency_binary_classification(y_binary, sampled_distributions) +function empirical_frequency_binary_classification(y_binary, sampled_distributions, n_bins) # Number of bins - n_bins = 20 + n_bins = n_bins #intervals boundaries int_bds = collect(range(0; stop=1, length=n_bins + 1)) #bin centers From dce9bdb1e74d1f4819e40c4e3851e583b62a06a1 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sat, 29 Jun 2024 14:57:34 +0200 Subject: [PATCH 31/46] ops forgot default value --- src/utils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 14627b8a..af86a293 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -57,7 +57,7 @@ Inputs: Outputs: - counts: an array cointaining the empirical frequencies for each quantile interval. """ -function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins) +function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins = 20) quantiles = collect(range(0; stop=1, length=n_bins + 1)) quantiles_matrix = hcat( [quantile(samples, quantiles) for samples in sampled_distributions]... @@ -109,7 +109,7 @@ Outputs: - bin_centers: array with the centers of the bins """ -function empirical_frequency_binary_classification(y_binary, sampled_distributions, n_bins) +function empirical_frequency_binary_classification(y_binary, sampled_distributions, n_bins = 20) # Number of bins n_bins = n_bins From b906c3bbc40adad24267f73ee1cc0884bd0ce996 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sat, 29 Jun 2024 14:57:34 +0200 Subject: [PATCH 32/46] ops forgot default value and removed a line --- src/utils.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index 14627b8a..f8d9b973 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -57,7 +57,7 @@ Inputs: Outputs: - counts: an array cointaining the empirical frequencies for each quantile interval. """ -function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins) +function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins = 20) quantiles = collect(range(0; stop=1, length=n_bins + 1)) quantiles_matrix = hcat( [quantile(samples, quantiles) for samples in sampled_distributions]... @@ -109,10 +109,7 @@ Outputs: - bin_centers: array with the centers of the bins """ -function empirical_frequency_binary_classification(y_binary, sampled_distributions, n_bins) - - # Number of bins - n_bins = n_bins +function empirical_frequency_binary_classification(y_binary, sampled_distributions, n_bins = 20) #intervals boundaries int_bds = collect(range(0; stop=1, length=n_bins + 1)) #bin centers From 32586189c70dc2231c6b7df8d712b5a527188df5 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sat, 29 Jun 2024 15:07:09 +0200 Subject: [PATCH 33/46] juliaformatter---- --- src/utils.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index f8d9b973..b4f08f6b 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -57,7 +57,7 @@ Inputs: Outputs: - counts: an array cointaining the empirical frequencies for each quantile interval. """ -function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins = 20) +function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=20) quantiles = collect(range(0; stop=1, length=n_bins + 1)) quantiles_matrix = hcat( [quantile(samples, quantiles) for samples in sampled_distributions]... @@ -109,7 +109,9 @@ Outputs: - bin_centers: array with the centers of the bins """ -function empirical_frequency_binary_classification(y_binary, sampled_distributions, n_bins = 20) +function empirical_frequency_binary_classification( + y_binary, sampled_distributions, n_bins=20 +) #intervals boundaries int_bds = collect(range(0; stop=1, length=n_bins + 1)) #bin centers From c86dc25b790d11f02405c8c6ac44cb107c4800fd Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Sun, 30 Jun 2024 23:04:19 +0200 Subject: [PATCH 34/46] fixed small error in pred_avg --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index b4f08f6b..d2377a7f 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -137,7 +137,7 @@ function empirical_frequency_binary_classification( #compute the empirical average and saved it in emp_avg in the j-th position push!(emp_avg, 1 / num_p_per_interval[j] * sum(y_binary[indices])) #TO DO: maybe substitute to bin_Centers? - push!(pred_avg, 1 / num_p_per_interval[j] * sum(class_probs[1, indices])) + push!(pred_avg, 1 / num_p_per_interval[j] * sum(class_probs[ indices])) end end #return the tuple From 3d2ebd67637c7b398cfd394079ee8091c1d39804 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Mon, 1 Jul 2024 00:45:51 +0200 Subject: [PATCH 35/46] fixed error in empirical_frequency_regression --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index d2377a7f..afac35e6 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -62,7 +62,7 @@ function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=20) quantiles_matrix = hcat( [quantile(samples, quantiles) for samples in sampled_distributions]... ) - n_rows = size(bounds_quantiles_matrix, 1) + n_rows = size(quantiles_matrix, 1) counts = [] for i in 1:n_rows From 4ab04f6bcb5637c8ce09f859a502d06d56eb3e35 Mon Sep 17 00:00:00 2001 From: pasquale c <343guiltyspark@outlook.it> Date: Mon, 1 Jul 2024 01:05:27 +0200 Subject: [PATCH 36/46] Update src/utils.jl yeah basically juliaformatter being a pain for nothing Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index afac35e6..0583d30d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -137,7 +137,7 @@ function empirical_frequency_binary_classification( #compute the empirical average and saved it in emp_avg in the j-th position push!(emp_avg, 1 / num_p_per_interval[j] * sum(y_binary[indices])) #TO DO: maybe substitute to bin_Centers? - push!(pred_avg, 1 / num_p_per_interval[j] * sum(class_probs[ indices])) + push!(pred_avg, 1 / num_p_per_interval[j] * sum(class_probs[indices])) end end #return the tuple From 267b8f4867662b8bb509e43c8b8641991081afdc Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Tue, 2 Jul 2024 03:19:39 +0200 Subject: [PATCH 37/46] docstrings fixes and predict update --- src/baselaplace/predicting.jl | 15 +++-- src/utils.jl | 111 ++++++++++++++++++---------------- 2 files changed, 69 insertions(+), 57 deletions(-) diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index d0c66ebd..80514428 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -61,10 +61,12 @@ Computes predictions from Bayesian neural network. - `predict_proba::Bool=true`: If `true` (default), returns probabilities for classification tasks. # Returns -For classification tasks: -- `fμ::AbstractArray`: Mean of the predictive distribution if link function is set to `:plugin`, otherwise the probit approximation. The output shape is column-major as in Flux. +For classification tasks, LaplaceRedux provides different options: + -`normal_distr::Distributions.Normal`:the array of Normal distributions computed by glm_predictive_distribution If the `link_approx` is set to :distribution + -`fμ::AbstractArray` Mean of the normal distribution if link_approx is set to :plugin + -`fμ::AbstractArray` The probit approximation if link_approx is set to :probit For regression tasks: -- `normal_distr::Distributions.Normal`:the array of Normal distributions computed by glm_predictive_distribution. The output shape is column-major as in Flux. +- `normal_distr::Distributions.Normal`:the array of Normal distributions computed by glm_predictive_distribution. # Examples @@ -93,6 +95,11 @@ function predict( # Classification: if la.likelihood == :classification + # Probit approximation + if link_approx == :distribution + z = normal_distr + end + # Probit approximation if link_approx == :probit z = probit(fμ, fvar) @@ -103,7 +110,7 @@ function predict( end # Sigmoid/Softmax - if predict_proba + if (predict_proba && link_approx != :distribution) if la.posterior.n_out == 1 p = Flux.sigmoid(z) else diff --git a/src/utils.jl b/src/utils.jl index 0583d30d..f4387bd7 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -41,21 +41,24 @@ function outdim(model::Chain)::Number return [size(p) for p in Flux.params(model)][end][1] end -""" +@doc raw""" empirical_frequency(Y_cal, sampled_distributions) -FOR REGRESSION MODELS. -Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the empirical frequency -phat_j = {y_t|F_t(y_t)<= p_j, t= 1,....,T}/T, where T is the number of calibration points, p_j is the confidence level and F_t is the -cumulative distribution function of the predicted distribution targeting y_t. -Source: https://arxiv.org/abs/1807.00263 - -Inputs: - - 'Y_cal': a vector of values y_t - - 'sampled_distributions': an array of sampled distributions F(x_t) stacked column-wise. - - 'n_bins': number of equally spaced bins to use. -Outputs: - - counts: an array cointaining the empirical frequencies for each quantile interval. +FOR REGRESSION MODELS. \ +Given a calibration dataset ``(x_t, y_t)`` for ``i ∈ {1,...,T}`` and an array of predicted distributions, the function calculates the empirical frequency +```math +p^hat_j = {y_t|F_t(y_t)<= p_j, t= 1,....,T}/T, +``` +where ``T`` is the number of calibration points, ``p_j`` is the confidence level and ``F_t`` is the +cumulative distribution function of the predicted distribution targeting ``y_t``. \ +Source: [Kuleshov, Fenner, Ermon 2018](https://arxiv.org/abs/1807.00263) + +Inputs: \ + - `Y_cal`: a vector of values ``y_t``\ + - `sampled_distributions`: an array of sampled distributions ``F(x_t)`` stacked column-wise.\ + - `n_bins`: number of equally spaced bins to use.\ +Outputs:\ + - `counts`: an array cointaining the empirical frequencies for each quantile interval. """ function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=20) quantiles = collect(range(0; stop=1, length=n_bins + 1)) @@ -71,42 +74,43 @@ function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=20) return counts end -""" +@doc raw""" sharpness(sampled_distributions) -FOR REGRESSION MODELS. -Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} and an array of predicted distributions, the function calculates the -sharpness of the predicted distributions, i.e., the average of the variances var(F_t) predicted by the forecaster for each x_t -Source: https://arxiv.org/abs/1807.00263 +FOR REGRESSION MODELS. \ +Given a calibration dataset ``(x_t, y_t)`` for ``i ∈ {1,...,T}`` and an array of predicted distributions, the function calculates the +sharpness of the predicted distributions, i.e., the average of the variances ``\sigma^2(F_t)`` predicted by the forecaster for each ``x_t``. \ +source: [Kuleshov, Fenner, Ermon 2018](https://arxiv.org/abs/1807.00263) -Inputs: - - sampled_distributions: an array of sampled distributions F(x_t) stacked column-wise. -Outputs: - - sharpness: a scalar that measure the level of sharpness of the regressor +Inputs: \ + - `sampled_distributions`: an array of sampled distributions ``F(x_t)`` stacked column-wise. \ +Outputs: \ + - `sharpness`: a scalar that measure the level of sharpness of the regressor """ function sharpness_regression(sampled_distributions) sharpness = mean(var.(sampled_distributions)) return sharpness end -""" +@doc raw""" empirical_frequency-classification(y_binary, sampled_distributions) -FOR BINARY CLASSIFICATION MODELS. -Given a calibration dataset (x_t, y_t) for i ∈ {1,...,T} let p_t= H(x_t)∈[0,1] be the forecasted probability. -We group the p_t into intervals I-j for j= 1,2,...,m that form a partition of [0,1]. The function computes -the observed average p_j= T^-1_j ∑_{t:p_t ∈ I_j} y_j in each interval I_j. -Source: https://arxiv.org/abs/1807.00263 - -Inputs: - - y_binary: the array of outputs y_t numerically coded: 1 for the target class, 0 for the null class. - - sampled_distributions: an array of sampled distributions stacked column-wise so that in the first row - there is the probability for the target class y_1 and in the second row the probability for the null class y_0. - - 'n_bins': number of equally spaced bins to use. -Outputs: - - num_p_per_interval: array with the number of probabilities falling within interval - - emp_avg: array with the observed empirical average per interval - - bin_centers: array with the centers of the bins +FOR BINARY CLASSIFICATION MODELS.\ +Given a calibration dataset ``(x_t, y_t)`` for ``i ∈ {1,...,T}`` let ``p_t= H(x_t)∈[0,1]`` be the forecasted probability. \ +We group the ``p_t`` into intervals ``I_j`` for ``j= 1,2,...,m`` that form a partition of [0,1]. +The function computes the observed average ``p_j= T^-1_j ∑_{t:p_t ∈ I_j} y_j`` in each interval ``I_j``. \ +Source: [Kuleshov, Fenner, Ermon 2018](https://arxiv.org/abs/1807.00263) + +Inputs: \ + - `y_binary`: the array of outputs ``y_t`` numerically coded: 1 for the target class, 0 for the null class. \ + - `sampled_distributions`: an array of sampled distributions stacked column-wise so that in the first row + there is the probability for the target class ``y_1`` and in the second row the probability for the null class ``y_0``. \ + - `n_bins`: number of equally spaced bins to use. + +Outputs: \ + - `num_p_per_interval`: array with the number of probabilities falling within interval. \ + - `emp_avg`: array with the observed empirical average per interval. \ + - `bin_centers`: array with the centers of the bins. """ function empirical_frequency_binary_classification( @@ -144,21 +148,22 @@ function empirical_frequency_binary_classification( return (num_p_per_interval, emp_avg, bin_centers) end -""" - sharpness-classification(y_binary,sampled_distributions) - -FOR BINARY CLASSIFICATION MODELS. -Assess the sharpness of the model by looking at the distribution of model predictions. When forecasts are sharp, -most predictions are close to 0 or 1; not sharp forecasters make predictions closer to 0.5. -Source: https://arxiv.org/abs/1807.00263 - -Inputs: - -y_binary: the array of outputs y_t numerically coded . 1 for the target class, 0 for the negative result. - -sampled_distributions: an array of sampled distributions stacked column-wise so that in the first row - there is the probability for the target class and in the second row the probability for the null class. - Outputs: - - mean_class_one: a scalar that measure the average prediction for the target class - - mean_class_zero: a scalar that measure the average prediction for the null class +@doc raw""" + sharpness_classification(y_binary,sampled_distributions) + +FOR BINARY CLASSIFICATION MODELS. \ +Assess the sharpness of the model by looking at the distribution of model predictions. +When forecasts are sharp, most predictions are close to either 0 or 1 \ +Source: [Kuleshov, Fenner, Ermon 2018](https://arxiv.org/abs/1807.00263) + +Inputs: \ + - `y_binary` : the array of outputs ``y_t`` numerically coded: 1 for the target class, 0 for the negative result. \ + - `sampled_distributions` : an array of sampled distributions stacked column-wise so that in the first row there is the probability for the target class ``y_1`` and in the second row the probability for the null class ``y_0``. \ + +Outputs: \ + - `mean_class_one` : a scalar that measure the average prediction for the target class \ + - `mean_class_zero` : a scalar that measure the average prediction for the null class + """ function sharpness_classification(y_binary, sampled_distributions) mean_class_one = mean(sampled_distributions[1, findall(y_binary .== 1)]) From d188daf2ec4f5f7047db59a10f760dddd883816d Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Tue, 2 Jul 2024 14:37:55 +0200 Subject: [PATCH 38/46] fixed typos --- src/utils.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/utils.jl b/src/utils.jl index f4387bd7..5c94b3a5 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -42,7 +42,7 @@ function outdim(model::Chain)::Number end @doc raw""" - empirical_frequency(Y_cal, sampled_distributions) + empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=20) FOR REGRESSION MODELS. \ Given a calibration dataset ``(x_t, y_t)`` for ``i ∈ {1,...,T}`` and an array of predicted distributions, the function calculates the empirical frequency @@ -75,7 +75,7 @@ function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=20) end @doc raw""" - sharpness(sampled_distributions) + sharpness_regression(sampled_distributions) FOR REGRESSION MODELS. \ Given a calibration dataset ``(x_t, y_t)`` for ``i ∈ {1,...,T}`` and an array of predicted distributions, the function calculates the @@ -93,7 +93,7 @@ function sharpness_regression(sampled_distributions) end @doc raw""" - empirical_frequency-classification(y_binary, sampled_distributions) + empirical_frequency_classification(y_binary, sampled_distributions) FOR BINARY CLASSIFICATION MODELS.\ Given a calibration dataset ``(x_t, y_t)`` for ``i ∈ {1,...,T}`` let ``p_t= H(x_t)∈[0,1]`` be the forecasted probability. \ From 270b70ae12893b2517e16b086e08ea7eb7b8d0c5 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Tue, 2 Jul 2024 17:34:37 +0200 Subject: [PATCH 39/46] moved sharpness functions units tests in calibration.jl. changed runtests.jl to add the new file. --- src/utils.jl | 16 ++++++++++++---- test/calibration.jl | 25 +++++++++++++++++++++++++ test/runtests.jl | 3 +++ 3 files changed, 40 insertions(+), 4 deletions(-) create mode 100644 test/calibration.jl diff --git a/src/utils.jl b/src/utils.jl index 5c94b3a5..b832154e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -60,7 +60,12 @@ Inputs: \ Outputs:\ - `counts`: an array cointaining the empirical frequencies for each quantile interval. """ -function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=20) +function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins::Int=20) + if n_bins <= 0 + throw(ArgumentError("n_bins must be a positive integer")) + elseif all(x -> x == 0 || x == 1, y_binary) + throw(ArgumentError("y_binary must be an array of 0 and 1")) + end quantiles = collect(range(0; stop=1, length=n_bins + 1)) quantiles_matrix = hcat( [quantile(samples, quantiles) for samples in sampled_distributions]... @@ -113,9 +118,12 @@ Outputs: \ - `bin_centers`: array with the centers of the bins. """ -function empirical_frequency_binary_classification( - y_binary, sampled_distributions, n_bins=20 -) +function empirical_frequency_binary_classification(y_binary, sampled_distributions, n_bins::Int=20) + if n_bins <= 0 + throw(ArgumentError("n_bins must be a positive integer")) + elseif all(x -> x == 0 || x == 1, y_binary) + throw(ArgumentError("y_binary must be an array of 0 and 1")) + end #intervals boundaries int_bds = collect(range(0; stop=1, length=n_bins + 1)) #bin centers diff --git a/test/calibration.jl b/test/calibration.jl new file mode 100644 index 00000000..d7c669b4 --- /dev/null +++ b/test/calibration.jl @@ -0,0 +1,25 @@ +using Statistics +using LaplaceRedux + +@testset "sharpness_classification tests" begin + y_binary = [0, 1, 0, 1, 1, 0, 1, 0] + sampled_distributions = [ + 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8; + 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 + ] + mean_class_one, mean_class_zero = sharpness_classification(y_binary, sampled_distributions) + @test mean_class_one ≈ mean(sampled_distributions[1,[2,4,5,7]]) + @test mean_class_zero ≈ mean(sampled_distributions[2,[1,3,6,8]]) + +end + + + +# Test for `sharpness_regression` function +@testset "sharpness_regression tests" begin + sampled_distributions = [[0.1, 0.2, 0.3, 0.7, 0.6], [0.2, 0.3, 0.4, 0.3 , 0.5 ], [0.3, 0.4, 0.5, 0.9, 0.2]] + mean_variance = mean(map(var, sampled_distributions)) + sharpness = sharpness_regression(sampled_distributions) + + @test sharpness ≈ mean_variance +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 5fa82a88..0459cf5a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -17,6 +17,9 @@ using Test @testset "Laplace" begin include("laplace.jl") end + @testset "Calibration Plots" begin + include("calibration.jl") + end if VERSION >= v"1.8.0" @testset "PyTorch Comparisons" begin From 33200635845e31456c29a1a787f0e3c47fe1dcce Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Tue, 2 Jul 2024 18:09:17 +0200 Subject: [PATCH 40/46] more sharpness unit tests --- test/calibration.jl | 41 ++++++++++++++++++++++++++++++++++++++++- 1 file changed, 40 insertions(+), 1 deletion(-) diff --git a/test/calibration.jl b/test/calibration.jl index d7c669b4..f7a63e2a 100644 --- a/test/calibration.jl +++ b/test/calibration.jl @@ -2,6 +2,17 @@ using Statistics using LaplaceRedux @testset "sharpness_classification tests" begin + + + # Test 1: Check that the function runs without errors and returns two scalars for a simple case + y_binary = [1, 0, 1, 0, 1] + sampled_distributions = [0.9 0.1 0.8 0.2 0.7; 0.1 0.9 0.2 0.8 0.3] # Sampled probabilities + mean_class_one, mean_class_zero = sharpness_classification(y_binary, sampled_distributions) + @test typeof(mean_class_one) <: Real # Check if mean_class_one is a scalar + @test typeof(mean_class_zero) <: Real # Check if mean_class_zero is a scalar + + + # Test 2: Check the function with a known input y_binary = [0, 1, 0, 1, 1, 0, 1, 0] sampled_distributions = [ 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8; @@ -10,6 +21,21 @@ using LaplaceRedux mean_class_one, mean_class_zero = sharpness_classification(y_binary, sampled_distributions) @test mean_class_one ≈ mean(sampled_distributions[1,[2,4,5,7]]) @test mean_class_zero ≈ mean(sampled_distributions[2,[1,3,6,8]]) + + # Test 3: Edge case with all ones in y_binary + y_binary_all_ones = [1, 1, 1] + sampled_distributions_all_ones = [0.8 0.9 0.7; 0.2 0.1 0.3] + mean_class_one_all_ones, mean_class_zero_all_ones = sharpness_classification(y_binary_all_ones, sampled_distributions_all_ones) + @test mean_class_one_all_ones == mean([0.8, 0.9, 0.7]) + @test isnan(mean_class_zero_all_ones) # Since there are no zeros in y_binary, the mean should be NaN + + # Test 4: Edge case with all zeros in y_binary + y_binary_all_zeros = [0, 0, 0] + sampled_distributions_all_zeros = [0.1 0.2 0.3; 0.9 0.8 0.7] + mean_class_one_all_zeros, mean_class_zero_all_zeros = sharpness_classification(y_binary_all_zeros, sampled_distributions_all_zeros) + @test mean_class_zero_all_zeros == mean([0.9, 0.8, 0.7]) + @test isnan(mean_class_one_all_zeros) # Since there are no ones in y_binary, the mean should be NaN + end @@ -17,9 +43,22 @@ end # Test for `sharpness_regression` function @testset "sharpness_regression tests" begin + + # Test 1: Check that the function runs without errors and returns a scalar for a simple case + sampled_distributions = [randn(100) for _ in 1:10] # Create 10 distributions, each with 100 samples + sharpness = sharpness_regression(sampled_distributions) + @test typeof(sharpness) <: Real # Check if the output is a scalar + + # Test 2: Check the function with a known input sampled_distributions = [[0.1, 0.2, 0.3, 0.7, 0.6], [0.2, 0.3, 0.4, 0.3 , 0.5 ], [0.3, 0.4, 0.5, 0.9, 0.2]] mean_variance = mean(map(var, sampled_distributions)) sharpness = sharpness_regression(sampled_distributions) - @test sharpness ≈ mean_variance + + # Test 3: Edge case with identical distributions + sampled_distributions_identical = [ones(100) for _ in 1:10] # Identical distributions, zero variance + sharpness_identical = sharpness_regression(sampled_distributions_identical) + @test sharpness_identical == 0.0 # Sharpness should be zero for identical distributions + + end \ No newline at end of file From 3750dbe5716d0ee5bccf2c5bea00a20d050e8da9 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Wed, 3 Jul 2024 00:14:46 +0200 Subject: [PATCH 41/46] fixes and more unit tests --- src/baselaplace/predicting.jl | 2 +- src/utils.jl | 16 +++++------ test/Project.toml | 1 + test/calibration.jl | 52 +++++++++++++++++++++++++++++++++++ 4 files changed, 62 insertions(+), 9 deletions(-) diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index 80514428..e0d657bc 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -44,7 +44,7 @@ function glm_predictive_distribution(la::AbstractLaplace, X::AbstractArray) fvar = functional_variance(la, 𝐉) fvar = reshape(fvar, size(fμ)...) fstd = sqrt.(fvar) - normal_distr = [Distributions.Normal(fμ[i], fstd[i]) for i in 1:size(fμ, 1)] + normal_distr = [Distributions.Normal(fμ[i], fstd[i]) for i in 1:size(fμ, 2)] return (normal_distr, fμ, fvar) end diff --git a/src/utils.jl b/src/utils.jl index b832154e..d304966c 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -60,18 +60,17 @@ Inputs: \ Outputs:\ - `counts`: an array cointaining the empirical frequencies for each quantile interval. """ -function empirical_frequency_regression(Y_cal, sampled_distributions, n_bins::Int=20) +function empirical_frequency_regression(Y_cal, sampled_distributions; n_bins::Int=20) if n_bins <= 0 throw(ArgumentError("n_bins must be a positive integer")) - elseif all(x -> x == 0 || x == 1, y_binary) - throw(ArgumentError("y_binary must be an array of 0 and 1")) end - quantiles = collect(range(0; stop=1, length=n_bins + 1)) + n_edges = n_bins + 1 + quantiles = collect(range(0; stop=1, length=n_edges)) quantiles_matrix = hcat( [quantile(samples, quantiles) for samples in sampled_distributions]... ) n_rows = size(quantiles_matrix, 1) - counts = [] + counts = Float64[] for i in 1:n_rows push!(counts, sum(Y_cal .<= quantiles_matrix[i, :]) / length(Y_cal)) @@ -118,14 +117,15 @@ Outputs: \ - `bin_centers`: array with the centers of the bins. """ -function empirical_frequency_binary_classification(y_binary, sampled_distributions, n_bins::Int=20) +function empirical_frequency_binary_classification(y_binary, sampled_distributions; n_bins::Int=20) if n_bins <= 0 throw(ArgumentError("n_bins must be a positive integer")) - elseif all(x -> x == 0 || x == 1, y_binary) + elseif !all(x -> x == 0 || x == 1, y_binary) throw(ArgumentError("y_binary must be an array of 0 and 1")) end #intervals boundaries - int_bds = collect(range(0; stop=1, length=n_bins + 1)) + n_edges = n_bins + 1 + int_bds = collect(range(0; stop=1, length=n_edges)) #bin centers bin_centers = [(int_bds[i] + int_bds[i + 1]) / 2 for i in 1:(length(int_bds) - 1)] #initialize list for empirical averages per interval diff --git a/test/Project.toml b/test/Project.toml index 7af07ce3..e0aee2a2 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -3,6 +3,7 @@ Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/test/calibration.jl b/test/calibration.jl index f7a63e2a..a033ab87 100644 --- a/test/calibration.jl +++ b/test/calibration.jl @@ -1,5 +1,6 @@ using Statistics using LaplaceRedux +using Distributions @testset "sharpness_classification tests" begin @@ -61,4 +62,55 @@ end @test sharpness_identical == 0.0 # Sharpness should be zero for identical distributions +end + + + +# Test for `empirical_frequency_regression` function +@testset "empirical_frequency_regression tests" begin + # Test 1: Check that the function runs without errors and returns an array for a simple case + Y_cal = [0.5, 1.5, 2.5, 3.5, 4.5] + sampled_distributions = [rand(Distributions.Normal(1, 1.0),6) for _ in 1:5] + counts = empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=10) + @test typeof(counts) == Array{Float64, 1} # Check if the output is an array of Float64 + + # Test 2: Check the function with a known input + #to do + + # Test 3: Invalid n_bins input + Y_cal = [0.5, 1.5, 2.5, 3.5, 4.5] + sampled_distributions = [rand(Distributions.Normal(1, 1.0),6) for _ in 1:5] + @test_throws ArgumentError empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=0) + +end + + + +# Test for `empirical_frequency_binary_classification` function +@testset "empirical_frequency_binary_classification tests" begin + # Test 1: Check that the function runs without errors and returns an array for a simple case + y_binary = rand(0:1, 10) + sampled_distributions = rand(Normal(0.5, 0.1), 10, 6) + n_bins = 4 + num_p_per_interval, emp_avg, bin_centers = empirical_frequency_binary_classification(y_binary, sampled_distributions, n_bins=n_bins) + @test length(num_p_per_interval) == n_bins + @test length(emp_avg) == n_bins + @test length(bin_centers) == n_bins + + # Test 2: Check the function with a known input + + #to do + + + + # Test 4: Invalid Y_cal input + Y_cal = [0, 1, 0, 1.2, 4] + sampled_distributions = rand(Normal(0.5, 0.1), 5, 6) + @test_throws ArgumentError empirical_frequency_binary_classification(Y_cal, sampled_distributions, n_bins=10) + + + # Test 5: Invalid n_bins input + Y_cal = rand(0:1, 5) + sampled_distributions = rand(Normal(0.5, 0.1), 5, 6) + @test_throws ArgumentError empirical_frequency_binary_classification(Y_cal, sampled_distributions, n_bins=0) end \ No newline at end of file From 39d4bdc5273a9291041ae98b0cec0f0a41c33130 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Wed, 3 Jul 2024 03:22:33 +0200 Subject: [PATCH 42/46] small stuff --- src/baselaplace/predicting.jl | 1 + test/calibration.jl | 5 +++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index e0d657bc..2f7c63ef 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -98,6 +98,7 @@ function predict( # Probit approximation if link_approx == :distribution z = normal_distr + #to complete e transform in a distribution.categorical .... end # Probit approximation diff --git a/test/calibration.jl b/test/calibration.jl index a033ab87..5b7f2971 100644 --- a/test/calibration.jl +++ b/test/calibration.jl @@ -73,6 +73,7 @@ end sampled_distributions = [rand(Distributions.Normal(1, 1.0),6) for _ in 1:5] counts = empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=10) @test typeof(counts) == Array{Float64, 1} # Check if the output is an array of Float64 + @test length(counts) == 10 # Test 2: Check the function with a known input #to do @@ -103,13 +104,13 @@ end - # Test 4: Invalid Y_cal input + # Test 3: Invalid Y_cal input Y_cal = [0, 1, 0, 1.2, 4] sampled_distributions = rand(Normal(0.5, 0.1), 5, 6) @test_throws ArgumentError empirical_frequency_binary_classification(Y_cal, sampled_distributions, n_bins=10) - # Test 5: Invalid n_bins input + # Test 4: Invalid n_bins input Y_cal = rand(0:1, 5) sampled_distributions = rand(Normal(0.5, 0.1), 5, 6) @test_throws ArgumentError empirical_frequency_binary_classification(Y_cal, sampled_distributions, n_bins=0) From 56c3b662e08589e1d134cea5fc5bb30764f4cd36 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Wed, 3 Jul 2024 06:50:45 +0200 Subject: [PATCH 43/46] fix. there is still an issue with the shape of the input to use. --- docs/src/tutorials/calibration_plots.md | 16 ++++++++++++++++ test/calibration.jl | 11 ++++++----- 2 files changed, 22 insertions(+), 5 deletions(-) create mode 100644 docs/src/tutorials/calibration_plots.md diff --git a/docs/src/tutorials/calibration_plots.md b/docs/src/tutorials/calibration_plots.md new file mode 100644 index 00000000..ef56d913 --- /dev/null +++ b/docs/src/tutorials/calibration_plots.md @@ -0,0 +1,16 @@ + +``` @meta +CurrentModule = LaplaceRedux +``` + + +# Calibration Plots + +## Libraries + +``` julia +using Pkg; Pkg.activate("docs") +# Import libraries +using Flux, Plots, TaijaPlotting, Random, Statistics, LaplaceRedux, LinearAlgebra +theme(:lime) +``` \ No newline at end of file diff --git a/test/calibration.jl b/test/calibration.jl index 5b7f2971..6db564cd 100644 --- a/test/calibration.jl +++ b/test/calibration.jl @@ -70,10 +70,11 @@ end @testset "empirical_frequency_regression tests" begin # Test 1: Check that the function runs without errors and returns an array for a simple case Y_cal = [0.5, 1.5, 2.5, 3.5, 4.5] + n_bins = 10 sampled_distributions = [rand(Distributions.Normal(1, 1.0),6) for _ in 1:5] - counts = empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=10) + counts = empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=n_bins) @test typeof(counts) == Array{Float64, 1} # Check if the output is an array of Float64 - @test length(counts) == 10 + @test length(counts) == n_bins + 1 # Test 2: Check the function with a known input #to do @@ -91,7 +92,7 @@ end @testset "empirical_frequency_binary_classification tests" begin # Test 1: Check that the function runs without errors and returns an array for a simple case y_binary = rand(0:1, 10) - sampled_distributions = rand(Normal(0.5, 0.1), 10, 6) + sampled_distributions = rand(2,10) n_bins = 4 num_p_per_interval, emp_avg, bin_centers = empirical_frequency_binary_classification(y_binary, sampled_distributions, n_bins=n_bins) @test length(num_p_per_interval) == n_bins @@ -106,12 +107,12 @@ end # Test 3: Invalid Y_cal input Y_cal = [0, 1, 0, 1.2, 4] - sampled_distributions = rand(Normal(0.5, 0.1), 5, 6) + sampled_distributions = rand(2,5) @test_throws ArgumentError empirical_frequency_binary_classification(Y_cal, sampled_distributions, n_bins=10) # Test 4: Invalid n_bins input Y_cal = rand(0:1, 5) - sampled_distributions = rand(Normal(0.5, 0.1), 5, 6) + sampled_distributions = rand(2,5) @test_throws ArgumentError empirical_frequency_binary_classification(Y_cal, sampled_distributions, n_bins=0) end \ No newline at end of file From 908c804e56b10eeb664c64dc8c0f6aa66ccf11aa Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Thu, 4 Jul 2024 05:51:33 +0200 Subject: [PATCH 44/46] fixed logit.md ,moved functions to new file, removed changes to predict ( will address it in a different issue) fixed docstring for calibration.jl --- docs/src/tutorials/logit.md | 1 + src/LaplaceRedux.jl | 5 +- src/baselaplace/predicting.jl | 9 +-- src/calibration_functions.jl | 141 ++++++++++++++++++++++++++++++++++ src/utils.jl | 138 --------------------------------- test/calibration.jl | 80 +++++++++---------- 6 files changed, 187 insertions(+), 187 deletions(-) create mode 100644 src/calibration_functions.jl diff --git a/docs/src/tutorials/logit.md b/docs/src/tutorials/logit.md index e5185187..bee5f5d5 100644 --- a/docs/src/tutorials/logit.md +++ b/docs/src/tutorials/logit.md @@ -81,3 +81,4 @@ p_untuned = plot(la_untuned, X, ys; title="LA - raw (λ=$(unique(diag(la_untuned p_laplace = plot(la, X, ys; title="LA - tuned (λ=$(round(unique(diag(la.prior.P₀))[1],digits=2)))", clim=(0,1), zoom=zoom) plot(p_plugin, p_untuned, p_laplace, layout=(1,3), size=(1700,400)) ``` +![](logit_files/figure-commonmark/cell-output-1.svg) \ No newline at end of file diff --git a/src/LaplaceRedux.jl b/src/LaplaceRedux.jl index 619a9d3b..63f969f1 100644 --- a/src/LaplaceRedux.jl +++ b/src/LaplaceRedux.jl @@ -1,9 +1,10 @@ module LaplaceRedux - -include("utils.jl") +include("calibration_functions.jl") export empirical_frequency_binary_classification, sharpness_classification, empirical_frequency_regression, sharpness_regression +include("utils.jl") + include("data/Data.jl") using .Data diff --git a/src/baselaplace/predicting.jl b/src/baselaplace/predicting.jl index 2f7c63ef..9891e0fb 100644 --- a/src/baselaplace/predicting.jl +++ b/src/baselaplace/predicting.jl @@ -62,7 +62,6 @@ Computes predictions from Bayesian neural network. # Returns For classification tasks, LaplaceRedux provides different options: - -`normal_distr::Distributions.Normal`:the array of Normal distributions computed by glm_predictive_distribution If the `link_approx` is set to :distribution -`fμ::AbstractArray` Mean of the normal distribution if link_approx is set to :plugin -`fμ::AbstractArray` The probit approximation if link_approx is set to :probit For regression tasks: @@ -95,12 +94,6 @@ function predict( # Classification: if la.likelihood == :classification - # Probit approximation - if link_approx == :distribution - z = normal_distr - #to complete e transform in a distribution.categorical .... - end - # Probit approximation if link_approx == :probit z = probit(fμ, fvar) @@ -111,7 +104,7 @@ function predict( end # Sigmoid/Softmax - if (predict_proba && link_approx != :distribution) + if (predict_proba) if la.posterior.n_out == 1 p = Flux.sigmoid(z) else diff --git a/src/calibration_functions.jl b/src/calibration_functions.jl new file mode 100644 index 00000000..fc7065da --- /dev/null +++ b/src/calibration_functions.jl @@ -0,0 +1,141 @@ +using Statistics +@doc raw""" + empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=20) + +FOR REGRESSION MODELS. \ +Given a calibration dataset ``(x_t, y_t)`` for ``i ∈ {1,...,T}`` and an array of predicted distributions, the function calculates the empirical frequency +```math +p^hat_j = {y_t|F_t(y_t)<= p_j, t= 1,....,T}/T, +``` +where ``T`` is the number of calibration points, ``p_j`` is the confidence level and ``F_t`` is the +cumulative distribution function of the predicted distribution targeting ``y_t``. \ +Source: [Kuleshov, Fenner, Ermon 2018](https://arxiv.org/abs/1807.00263) + +Inputs: \ + - `Y_cal`: a vector of values ``y_t``\ + - `sampled_distributions`:a Vector{Vector{Float64}} of sampled distributions ``F(x_t)`` stacked row-wise.\ + For example [rand(distr,50) for distr in LaplaceRedux.predict(la,X)] + - `n_bins`: number of equally spaced bins to use.\ +Outputs:\ + - `counts`: an array cointaining the empirical frequencies for each quantile interval. +""" +function empirical_frequency_regression(Y_cal, sampled_distributions; n_bins::Int=20) + if n_bins <= 0 + throw(ArgumentError("n_bins must be a positive integer")) + end + n_edges = n_bins + 1 + quantiles = collect(range(0; stop=1, length=n_edges)) + quantiles_matrix = hcat( + [quantile(samples, quantiles) for samples in sampled_distributions]... + ) + n_rows = size(quantiles_matrix, 1) + counts = Float64[] + + for i in 1:n_rows + push!(counts, sum(Y_cal .<= quantiles_matrix[i, :]) / length(Y_cal)) + end + return counts +end + +@doc raw""" + sharpness_regression(sampled_distributions) + +FOR REGRESSION MODELS. \ +Given a calibration dataset ``(x_t, y_t)`` for ``i ∈ {1,...,T}`` and an array of predicted distributions, the function calculates the +sharpness of the predicted distributions, i.e., the average of the variances ``\sigma^2(F_t)`` predicted by the forecaster for each ``x_t``. \ +source: [Kuleshov, Fenner, Ermon 2018](https://arxiv.org/abs/1807.00263) + +Inputs: \ + - `sampled_distributions`: an array of sampled distributions ``F(x_t)`` stacked column-wise. \ +Outputs: \ + - `sharpness`: a scalar that measure the level of sharpness of the regressor +""" +function sharpness_regression(sampled_distributions) + sharpness = mean(var.(sampled_distributions)) + return sharpness +end + +@doc raw""" + empirical_frequency_classification(y_binary, sampled_distributions) + +FOR BINARY CLASSIFICATION MODELS.\ +Given a calibration dataset ``(x_t, y_t)`` for ``i ∈ {1,...,T}`` let ``p_t= H(x_t)∈[0,1]`` be the forecasted probability. \ +We group the ``p_t`` into intervals ``I_j`` for ``j= 1,2,...,m`` that form a partition of [0,1]. +The function computes the observed average ``p_j= T^-1_j ∑_{t:p_t ∈ I_j} y_j`` in each interval ``I_j``. \ +Source: [Kuleshov, Fenner, Ermon 2018](https://arxiv.org/abs/1807.00263) + +Inputs: \ + - `y_binary`: the array of outputs ``y_t`` numerically coded: 1 for the target class, 0 for the null class. \ + - `sampled_distributions`: an array of sampled distributions stacked column-wise so that in the first row + there is the probability for the target class ``y_1`` and in the second row the probability for the null class ``y_0``. \ + - `n_bins`: number of equally spaced bins to use. + +Outputs: \ + - `num_p_per_interval`: array with the number of probabilities falling within interval. \ + - `emp_avg`: array with the observed empirical average per interval. \ + - `bin_centers`: array with the centers of the bins. + +""" +function empirical_frequency_binary_classification( + y_binary, sampled_distributions; n_bins::Int=20 +) + if n_bins <= 0 + throw(ArgumentError("n_bins must be a positive integer")) + elseif !all(x -> x == 0 || x == 1, y_binary) + throw(ArgumentError("y_binary must be an array of 0 and 1")) + end + #intervals boundaries + n_edges = n_bins + 1 + int_bds = collect(range(0; stop=1, length=n_edges)) + #bin centers + bin_centers = [(int_bds[i] + int_bds[i + 1]) / 2 for i in 1:(length(int_bds) - 1)] + #initialize list for empirical averages per interval + emp_avg = [] + #initialize list for predicted averages per interval + pred_avg = [] + # initialize list of number of probabilities falling within each intervals + num_p_per_interval = [] + #list of the predicted probabilities for the target class + class_probs = sampled_distributions[1, :] + # iterate over the bins + for j in 1:n_bins + push!(num_p_per_interval, sum(int_bds[j] .< class_probs .< int_bds[j + 1])) + if num_p_per_interval[j] == 0 + push!(emp_avg, 0) + push!(pred_avg, bin_centers[j]) + + else + # find the indices fo all istances for which class_probs fall withing the j-th interval + indices = findall(x -> int_bds[j] < x < int_bds[j + 1], class_probs) + #compute the empirical average and saved it in emp_avg in the j-th position + push!(emp_avg, 1 / num_p_per_interval[j] * sum(y_binary[indices])) + #TO DO: maybe substitute to bin_Centers? + push!(pred_avg, 1 / num_p_per_interval[j] * sum(class_probs[indices])) + end + end + #return the tuple + return (num_p_per_interval, emp_avg, bin_centers) +end + +@doc raw""" + sharpness_classification(y_binary,sampled_distributions) + +FOR BINARY CLASSIFICATION MODELS. \ +Assess the sharpness of the model by looking at the distribution of model predictions. +When forecasts are sharp, most predictions are close to either 0 or 1 \ +Source: [Kuleshov, Fenner, Ermon 2018](https://arxiv.org/abs/1807.00263) + +Inputs: \ + - `y_binary` : the array of outputs ``y_t`` numerically coded: 1 for the target class, 0 for the negative result. \ + - `sampled_distributions` : an array of sampled distributions stacked column-wise so that in the first row there is the probability for the target class ``y_1`` and in the second row the probability for the null class ``y_0``. \ + +Outputs: \ + - `mean_class_one` : a scalar that measure the average prediction for the target class \ + - `mean_class_zero` : a scalar that measure the average prediction for the null class + +""" +function sharpness_classification(y_binary, sampled_distributions) + mean_class_one = mean(sampled_distributions[1, findall(y_binary .== 1)]) + mean_class_zero = mean(sampled_distributions[2, findall(y_binary .== 0)]) + return mean_class_one, mean_class_zero +end diff --git a/src/utils.jl b/src/utils.jl index d304966c..4d2a9194 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -40,141 +40,3 @@ corresponding to the number of neurons on the last layer of the NN. function outdim(model::Chain)::Number return [size(p) for p in Flux.params(model)][end][1] end - -@doc raw""" - empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=20) - -FOR REGRESSION MODELS. \ -Given a calibration dataset ``(x_t, y_t)`` for ``i ∈ {1,...,T}`` and an array of predicted distributions, the function calculates the empirical frequency -```math -p^hat_j = {y_t|F_t(y_t)<= p_j, t= 1,....,T}/T, -``` -where ``T`` is the number of calibration points, ``p_j`` is the confidence level and ``F_t`` is the -cumulative distribution function of the predicted distribution targeting ``y_t``. \ -Source: [Kuleshov, Fenner, Ermon 2018](https://arxiv.org/abs/1807.00263) - -Inputs: \ - - `Y_cal`: a vector of values ``y_t``\ - - `sampled_distributions`: an array of sampled distributions ``F(x_t)`` stacked column-wise.\ - - `n_bins`: number of equally spaced bins to use.\ -Outputs:\ - - `counts`: an array cointaining the empirical frequencies for each quantile interval. -""" -function empirical_frequency_regression(Y_cal, sampled_distributions; n_bins::Int=20) - if n_bins <= 0 - throw(ArgumentError("n_bins must be a positive integer")) - end - n_edges = n_bins + 1 - quantiles = collect(range(0; stop=1, length=n_edges)) - quantiles_matrix = hcat( - [quantile(samples, quantiles) for samples in sampled_distributions]... - ) - n_rows = size(quantiles_matrix, 1) - counts = Float64[] - - for i in 1:n_rows - push!(counts, sum(Y_cal .<= quantiles_matrix[i, :]) / length(Y_cal)) - end - return counts -end - -@doc raw""" - sharpness_regression(sampled_distributions) - -FOR REGRESSION MODELS. \ -Given a calibration dataset ``(x_t, y_t)`` for ``i ∈ {1,...,T}`` and an array of predicted distributions, the function calculates the -sharpness of the predicted distributions, i.e., the average of the variances ``\sigma^2(F_t)`` predicted by the forecaster for each ``x_t``. \ -source: [Kuleshov, Fenner, Ermon 2018](https://arxiv.org/abs/1807.00263) - -Inputs: \ - - `sampled_distributions`: an array of sampled distributions ``F(x_t)`` stacked column-wise. \ -Outputs: \ - - `sharpness`: a scalar that measure the level of sharpness of the regressor -""" -function sharpness_regression(sampled_distributions) - sharpness = mean(var.(sampled_distributions)) - return sharpness -end - -@doc raw""" - empirical_frequency_classification(y_binary, sampled_distributions) - -FOR BINARY CLASSIFICATION MODELS.\ -Given a calibration dataset ``(x_t, y_t)`` for ``i ∈ {1,...,T}`` let ``p_t= H(x_t)∈[0,1]`` be the forecasted probability. \ -We group the ``p_t`` into intervals ``I_j`` for ``j= 1,2,...,m`` that form a partition of [0,1]. -The function computes the observed average ``p_j= T^-1_j ∑_{t:p_t ∈ I_j} y_j`` in each interval ``I_j``. \ -Source: [Kuleshov, Fenner, Ermon 2018](https://arxiv.org/abs/1807.00263) - -Inputs: \ - - `y_binary`: the array of outputs ``y_t`` numerically coded: 1 for the target class, 0 for the null class. \ - - `sampled_distributions`: an array of sampled distributions stacked column-wise so that in the first row - there is the probability for the target class ``y_1`` and in the second row the probability for the null class ``y_0``. \ - - `n_bins`: number of equally spaced bins to use. - -Outputs: \ - - `num_p_per_interval`: array with the number of probabilities falling within interval. \ - - `emp_avg`: array with the observed empirical average per interval. \ - - `bin_centers`: array with the centers of the bins. - -""" -function empirical_frequency_binary_classification(y_binary, sampled_distributions; n_bins::Int=20) - if n_bins <= 0 - throw(ArgumentError("n_bins must be a positive integer")) - elseif !all(x -> x == 0 || x == 1, y_binary) - throw(ArgumentError("y_binary must be an array of 0 and 1")) - end - #intervals boundaries - n_edges = n_bins + 1 - int_bds = collect(range(0; stop=1, length=n_edges)) - #bin centers - bin_centers = [(int_bds[i] + int_bds[i + 1]) / 2 for i in 1:(length(int_bds) - 1)] - #initialize list for empirical averages per interval - emp_avg = [] - #initialize list for predicted averages per interval - pred_avg = [] - # initialize list of number of probabilities falling within each intervals - num_p_per_interval = [] - #list of the predicted probabilities for the target class - class_probs = sampled_distributions[1, :] - # iterate over the bins - for j in 1:n_bins - push!(num_p_per_interval, sum(int_bds[j] .< class_probs .< int_bds[j + 1])) - if num_p_per_interval[j] == 0 - push!(emp_avg, 0) - push!(pred_avg, bin_centers[j]) - - else - # find the indices fo all istances for which class_probs fall withing the j-th interval - indices = findall(x -> int_bds[j] < x < int_bds[j + 1], class_probs) - #compute the empirical average and saved it in emp_avg in the j-th position - push!(emp_avg, 1 / num_p_per_interval[j] * sum(y_binary[indices])) - #TO DO: maybe substitute to bin_Centers? - push!(pred_avg, 1 / num_p_per_interval[j] * sum(class_probs[indices])) - end - end - #return the tuple - return (num_p_per_interval, emp_avg, bin_centers) -end - -@doc raw""" - sharpness_classification(y_binary,sampled_distributions) - -FOR BINARY CLASSIFICATION MODELS. \ -Assess the sharpness of the model by looking at the distribution of model predictions. -When forecasts are sharp, most predictions are close to either 0 or 1 \ -Source: [Kuleshov, Fenner, Ermon 2018](https://arxiv.org/abs/1807.00263) - -Inputs: \ - - `y_binary` : the array of outputs ``y_t`` numerically coded: 1 for the target class, 0 for the negative result. \ - - `sampled_distributions` : an array of sampled distributions stacked column-wise so that in the first row there is the probability for the target class ``y_1`` and in the second row the probability for the null class ``y_0``. \ - -Outputs: \ - - `mean_class_one` : a scalar that measure the average prediction for the target class \ - - `mean_class_zero` : a scalar that measure the average prediction for the null class - -""" -function sharpness_classification(y_binary, sampled_distributions) - mean_class_one = mean(sampled_distributions[1, findall(y_binary .== 1)]) - mean_class_zero = mean(sampled_distributions[2, findall(y_binary .== 0)]) - return mean_class_one, mean_class_zero -end diff --git a/test/calibration.jl b/test/calibration.jl index 6db564cd..649f0b43 100644 --- a/test/calibration.jl +++ b/test/calibration.jl @@ -4,44 +4,46 @@ using Distributions @testset "sharpness_classification tests" begin - # Test 1: Check that the function runs without errors and returns two scalars for a simple case y_binary = [1, 0, 1, 0, 1] sampled_distributions = [0.9 0.1 0.8 0.2 0.7; 0.1 0.9 0.2 0.8 0.3] # Sampled probabilities - mean_class_one, mean_class_zero = sharpness_classification(y_binary, sampled_distributions) + mean_class_one, mean_class_zero = sharpness_classification( + y_binary, sampled_distributions + ) @test typeof(mean_class_one) <: Real # Check if mean_class_one is a scalar @test typeof(mean_class_zero) <: Real # Check if mean_class_zero is a scalar - # Test 2: Check the function with a known input y_binary = [0, 1, 0, 1, 1, 0, 1, 0] sampled_distributions = [ - 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8; - 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 + 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 + 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 ] - mean_class_one, mean_class_zero = sharpness_classification(y_binary, sampled_distributions) - @test mean_class_one ≈ mean(sampled_distributions[1,[2,4,5,7]]) - @test mean_class_zero ≈ mean(sampled_distributions[2,[1,3,6,8]]) + mean_class_one, mean_class_zero = sharpness_classification( + y_binary, sampled_distributions + ) + @test mean_class_one ≈ mean(sampled_distributions[1, [2, 4, 5, 7]]) + @test mean_class_zero ≈ mean(sampled_distributions[2, [1, 3, 6, 8]]) # Test 3: Edge case with all ones in y_binary y_binary_all_ones = [1, 1, 1] sampled_distributions_all_ones = [0.8 0.9 0.7; 0.2 0.1 0.3] - mean_class_one_all_ones, mean_class_zero_all_ones = sharpness_classification(y_binary_all_ones, sampled_distributions_all_ones) + mean_class_one_all_ones, mean_class_zero_all_ones = sharpness_classification( + y_binary_all_ones, sampled_distributions_all_ones + ) @test mean_class_one_all_ones == mean([0.8, 0.9, 0.7]) @test isnan(mean_class_zero_all_ones) # Since there are no zeros in y_binary, the mean should be NaN # Test 4: Edge case with all zeros in y_binary y_binary_all_zeros = [0, 0, 0] sampled_distributions_all_zeros = [0.1 0.2 0.3; 0.9 0.8 0.7] - mean_class_one_all_zeros, mean_class_zero_all_zeros = sharpness_classification(y_binary_all_zeros, sampled_distributions_all_zeros) + mean_class_one_all_zeros, mean_class_zero_all_zeros = sharpness_classification( + y_binary_all_zeros, sampled_distributions_all_zeros + ) @test mean_class_zero_all_zeros == mean([0.9, 0.8, 0.7]) @test isnan(mean_class_one_all_zeros) # Since there are no ones in y_binary, the mean should be NaN - - end - - # Test for `sharpness_regression` function @testset "sharpness_regression tests" begin @@ -51,7 +53,9 @@ end @test typeof(sharpness) <: Real # Check if the output is a scalar # Test 2: Check the function with a known input - sampled_distributions = [[0.1, 0.2, 0.3, 0.7, 0.6], [0.2, 0.3, 0.4, 0.3 , 0.5 ], [0.3, 0.4, 0.5, 0.9, 0.2]] + sampled_distributions = [ + [0.1, 0.2, 0.3, 0.7, 0.6], [0.2, 0.3, 0.4, 0.3, 0.5], [0.3, 0.4, 0.5, 0.9, 0.2] + ] mean_variance = mean(map(var, sampled_distributions)) sharpness = sharpness_regression(sampled_distributions) @test sharpness ≈ mean_variance @@ -60,41 +64,38 @@ end sampled_distributions_identical = [ones(100) for _ in 1:10] # Identical distributions, zero variance sharpness_identical = sharpness_regression(sampled_distributions_identical) @test sharpness_identical == 0.0 # Sharpness should be zero for identical distributions - - end - - # Test for `empirical_frequency_regression` function @testset "empirical_frequency_regression tests" begin # Test 1: Check that the function runs without errors and returns an array for a simple case Y_cal = [0.5, 1.5, 2.5, 3.5, 4.5] n_bins = 10 - sampled_distributions = [rand(Distributions.Normal(1, 1.0),6) for _ in 1:5] - counts = empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=n_bins) - @test typeof(counts) == Array{Float64, 1} # Check if the output is an array of Float64 - @test length(counts) == n_bins + 1 + sampled_distributions = [rand(Distributions.Normal(1, 1.0), 6) for _ in 1:5] + counts = empirical_frequency_regression(Y_cal, sampled_distributions; n_bins=n_bins) + @test typeof(counts) == Array{Float64,1} # Check if the output is an array of Float64 + @test length(counts) == n_bins + 1 # Test 2: Check the function with a known input #to do # Test 3: Invalid n_bins input Y_cal = [0.5, 1.5, 2.5, 3.5, 4.5] - sampled_distributions = [rand(Distributions.Normal(1, 1.0),6) for _ in 1:5] - @test_throws ArgumentError empirical_frequency_regression(Y_cal, sampled_distributions, n_bins=0) - + sampled_distributions = [rand(Distributions.Normal(1, 1.0), 6) for _ in 1:5] + @test_throws ArgumentError empirical_frequency_regression( + Y_cal, sampled_distributions, n_bins=0 + ) end - - # Test for `empirical_frequency_binary_classification` function @testset "empirical_frequency_binary_classification tests" begin # Test 1: Check that the function runs without errors and returns an array for a simple case y_binary = rand(0:1, 10) - sampled_distributions = rand(2,10) + sampled_distributions = rand(2, 10) n_bins = 4 - num_p_per_interval, emp_avg, bin_centers = empirical_frequency_binary_classification(y_binary, sampled_distributions, n_bins=n_bins) + num_p_per_interval, emp_avg, bin_centers = empirical_frequency_binary_classification( + y_binary, sampled_distributions; n_bins=n_bins + ) @test length(num_p_per_interval) == n_bins @test length(emp_avg) == n_bins @test length(bin_centers) == n_bins @@ -103,16 +104,17 @@ end #to do - - # Test 3: Invalid Y_cal input - Y_cal = [0, 1, 0, 1.2, 4] - sampled_distributions = rand(2,5) - @test_throws ArgumentError empirical_frequency_binary_classification(Y_cal, sampled_distributions, n_bins=10) - + Y_cal = [0, 1, 0, 1.2, 4] + sampled_distributions = rand(2, 5) + @test_throws ArgumentError empirical_frequency_binary_classification( + Y_cal, sampled_distributions, n_bins=10 + ) # Test 4: Invalid n_bins input Y_cal = rand(0:1, 5) - sampled_distributions = rand(2,5) - @test_throws ArgumentError empirical_frequency_binary_classification(Y_cal, sampled_distributions, n_bins=0) -end \ No newline at end of file + sampled_distributions = rand(2, 5) + @test_throws ArgumentError empirical_frequency_binary_classification( + Y_cal, sampled_distributions, n_bins=0 + ) +end From f4688035e54f418961905e446f673a3641e97d8d Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Thu, 4 Jul 2024 05:53:05 +0200 Subject: [PATCH 45/46] removed calibration_plots.md --- docs/src/tutorials/calibration_plots.md | 16 ---------------- 1 file changed, 16 deletions(-) delete mode 100644 docs/src/tutorials/calibration_plots.md diff --git a/docs/src/tutorials/calibration_plots.md b/docs/src/tutorials/calibration_plots.md deleted file mode 100644 index ef56d913..00000000 --- a/docs/src/tutorials/calibration_plots.md +++ /dev/null @@ -1,16 +0,0 @@ - -``` @meta -CurrentModule = LaplaceRedux -``` - - -# Calibration Plots - -## Libraries - -``` julia -using Pkg; Pkg.activate("docs") -# Import libraries -using Flux, Plots, TaijaPlotting, Random, Statistics, LaplaceRedux, LinearAlgebra -theme(:lime) -``` \ No newline at end of file From 459b2feb380820e16a6bcfab529557e17837ac14 Mon Sep 17 00:00:00 2001 From: "pasquale c." <343guiltyspark@outlook.it> Date: Thu, 4 Jul 2024 08:06:59 +0200 Subject: [PATCH 46/46] test plot --- docs/src/tutorials/calibration.md | 9 +++ docs/src/tutorials/regression.md | 8 +++ .../figure-commonmark/miscalibration.svg | 56 +++++++++++++++++++ 3 files changed, 73 insertions(+) create mode 100644 docs/src/tutorials/calibration.md create mode 100644 docs/src/tutorials/regression_files/figure-commonmark/miscalibration.svg diff --git a/docs/src/tutorials/calibration.md b/docs/src/tutorials/calibration.md new file mode 100644 index 00000000..19db8ba2 --- /dev/null +++ b/docs/src/tutorials/calibration.md @@ -0,0 +1,9 @@ +# Uncertainty Calibration +## The issue of calibrated uncertainty distributions +Bayesian methods offer a general framework for quantifying uncertainty. However, due to model misspecification and the use of approximate inference, Bayesian uncertainty estimates are often inaccurate: for example, a 90% credible interval may not contain the true outcome 90% of the time. A model is considered calibrated when uncertainty estimates, such as Bayesian credible intervals, accurately reflect the true likelihood of outcomes. In other words, a 90% credible interval is calibrated if it contains the true outcome approximately 90% of the time, thereby indicating the reliability and accuracy of the inference method. In other words, a good forecaster must be calibrated. Perfect calibration + + +## Calibration Plots + + +yadda yadda \ No newline at end of file diff --git a/docs/src/tutorials/regression.md b/docs/src/tutorials/regression.md index eb660591..d2cef2fa 100644 --- a/docs/src/tutorials/regression.md +++ b/docs/src/tutorials/regression.md @@ -132,3 +132,11 @@ plot(la, X, y; zoom=-5, size=(400,400)) Scatter: 8.497215713339543 ![](regression_files/figure-commonmark/cell-7-output-5.svg) + + +## Calibration Plot +Once the prior precision has been optimized it is possible to evaluate the quality of the predictive distribution +obtained through a calibration plot [Link text Here](https://link-url-here.org) [add cross link]. + +![](regression_files/figure-commonmark/miscalibration.svg) + diff --git a/docs/src/tutorials/regression_files/figure-commonmark/miscalibration.svg b/docs/src/tutorials/regression_files/figure-commonmark/miscalibration.svg new file mode 100644 index 00000000..9ca6b50d --- /dev/null +++ b/docs/src/tutorials/regression_files/figure-commonmark/miscalibration.svg @@ -0,0 +1,56 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +