Skip to content

Commit

Permalink
Merge #19
Browse files Browse the repository at this point in the history
19: Update CalibrationErrors r=devmotion a=devmotion



Co-authored-by: David Widmann <dev+git@devmotion.de>
  • Loading branch information
bors[bot] and devmotion committed May 7, 2020
2 parents d3563f1 + c2c675d commit b621b3b
Show file tree
Hide file tree
Showing 14 changed files with 339 additions and 295 deletions.
4 changes: 2 additions & 2 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -14,10 +14,10 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"

[compat]
CalibrationErrors = "0.4"
CalibrationErrors = "0.5"
HypothesisTests = "0.8, 0.9, 0.10"
Parameters = "0.12"
StatsBase = "0.32"
StatsBase = "0.32, 0.33"
StatsFuns = "0.8, 0.9"
julia = "1.1"

Expand Down
6 changes: 3 additions & 3 deletions src/CalibrationTests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,15 +12,15 @@ using Random
using CalibrationErrors: CalibrationErrorEstimator, SKCE, unsafe_skce_eval

export ConsistencyTest
export DistributionFreeTest, AsymptoticLinearTest, AsymptoticQuadraticTest
export DistributionFreeSKCETest, AsymptoticBlockSKCETest, AsymptoticSKCETest

# re-export
export pvalue, confint

include("consistency.jl")

include("skce/asymptotic.jl")
include("skce/asymptotic_block.jl")
include("skce/distribution_free.jl")
include("skce/asymptotic_linear.jl")
include("skce/asymptotic_quadratic.jl")

end # module
3 changes: 1 addition & 2 deletions src/consistency.jl
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,7 @@ end

HypothesisTests.default_tail(::ConsistencyTest) = :right

HypothesisTests.testname(::ConsistencyTest) =
"Calibration test based on consistency resampling"
HypothesisTests.testname(::ConsistencyTest) = "Consistency resampling test"

HypothesisTests.pvalue(test::ConsistencyTest; kwargs...) =
pvalue(Random.GLOBAL_RNG, test; kwargs...)
Expand Down
29 changes: 14 additions & 15 deletions src/skce/asymptotic_quadratic.jl → src/skce/asymptotic.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
struct AsymptoticQuadraticTest{K<:Kernel,P,T,E,V} <: HypothesisTests.HypothesisTest
struct AsymptoticSKCETest{K<:Kernel,P,T,E,V} <: HypothesisTests.HypothesisTest
"""Kernel."""
kernel::K
"""Predictions."""
Expand All @@ -11,42 +11,41 @@ struct AsymptoticQuadraticTest{K<:Kernel,P,T,E,V} <: HypothesisTests.HypothesisT
statistic::V
end

AsymptoticQuadraticTest(skce::QuadraticUnbiasedSKCE, data...; kwargs...) =
AsymptoticQuadraticTest(skce.kernel, data...; kwargs...)
AsymptoticSKCETest(skce::UnbiasedSKCE, data...; kwargs...) =
AsymptoticSKCETest(skce.kernel, data...; kwargs...)

AsymptoticQuadraticTest(kernel1::Kernel, kernel2::Kernel, data...; kwargs...) =
AsymptoticQuadraticTest(TensorProductKernel(kernel1, kernel2), data...; kwargs...)
AsymptoticSKCETest(kernel1::Kernel, kernel2::Kernel, data...; kwargs...) =
AsymptoticSKCETest(TensorProduct(kernel1, kernel2), data...; kwargs...)

function AsymptoticQuadraticTest(kernel::Kernel, data...)
function AsymptoticSKCETest(kernel::Kernel, data...)
# obtain the predictions and targets
predictions, targets = CalibrationErrors.predictions_targets(data...)

# compute the calibration error estimate and the test statistic
estimate, statistic = estimate_statistic(kernel, predictions, targets)

AsymptoticQuadraticTest(kernel, predictions, targets, estimate, statistic)
AsymptoticSKCETest(kernel, predictions, targets, estimate, statistic)
end

# HypothesisTests interface

HypothesisTests.default_tail(::AsymptoticQuadraticTest) = :right
HypothesisTests.default_tail(::AsymptoticSKCETest) = :right

HypothesisTests.pvalue(test::AsymptoticQuadraticTest; kwargs...) =
HypothesisTests.pvalue(test::AsymptoticSKCETest; kwargs...) =
pvalue(Random.GLOBAL_RNG, test; kwargs...)
function HypothesisTests.pvalue(rng::AbstractRNG, test::AsymptoticQuadraticTest;
function HypothesisTests.pvalue(rng::AbstractRNG, test::AsymptoticSKCETest;
bootstrap_iters::Int = 1_000)
bootstrap_ccdf(rng, test, bootstrap_iters)
end

HypothesisTests.testname(::AsymptoticQuadraticTest) =
"Asymptotic calibration test based on the quadratic unbiased SKCE estimator"
HypothesisTests.testname(::AsymptoticSKCETest) = "Asymptotic SKCE test"

# parameter of interest: name, value under H0, point estimate
function HypothesisTests.population_param_of_interest(test::AsymptoticQuadraticTest)
function HypothesisTests.population_param_of_interest(test::AsymptoticSKCETest)
"SKCE", zero(test.estimate), test.estimate
end

function HypothesisTests.show_params(io::IO, test::AsymptoticQuadraticTest, ident = "")
function HypothesisTests.show_params(io::IO, test::AsymptoticSKCETest, ident = "")
println(io, ident, "test statistic: $(test.statistic)")
end

Expand Down Expand Up @@ -103,7 +102,7 @@ function estimate_statistic(kernel::Kernel,
end

# estimate the ccdf using bootstrapping
function bootstrap_ccdf(rng::AbstractRNG, test::AsymptoticQuadraticTest,
function bootstrap_ccdf(rng::AbstractRNG, test::AsymptoticSKCETest,
bootstrap_iters::Int)
@unpack kernel, predictions, targets, statistic = test

Expand Down
110 changes: 110 additions & 0 deletions src/skce/asymptotic_block.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,110 @@
struct AsymptoticBlockSKCETest{K<:Kernel,E,S,Z} <: HypothesisTests.ZTest
"""Kernel."""
kernel::K
"""Number of observations per block."""
blocksize::Int
"""Number of blocks of observations."""
nblocks::Int
"""Calibration error estimate (average of evaluations of blocks of observations)."""
estimate::E
"""Standard error of evaluations of blocks of observations."""
stderr::S
"""z-statistic."""
z::Z
end

function AsymptoticBlockSKCETest(skce::BlockUnbiasedSKCE, data...)
return AsymptoticBlockSKCETest(skce.kernel, skce.blocksize, data...)
end

function AsymptoticBlockSKCETest(
κpredictions::Kernel,
κtargets::Kernel,
args...
)
return AsymptoticBlockSKCETest(TensorProduct(κpredictions, κtargets), args...)
end

function AsymptoticBlockSKCETest(kernel::Kernel, blocksize::Int, data...)
# obtain predictions and targets
predictions, targets = CalibrationErrors.predictions_targets(data...)

# obtain number of samples
nsamples = length(predictions)
nsamples blocksize || error("there must be at least ", blocksize, " samples")

# compute number of blocks
nblocks = nsamples ÷ blocksize

# evaluate U-statistic of the first block
istart = 1
iend = blocksize
x = CalibrationErrors.unsafe_unbiasedskce(kernel, predictions, targets, istart, iend)

# initialize the estimate and the sum of squares
estimate = x / 1
S = zero(x)

# for all other blocks
for b in 2:nblocks
# evaluate U-statistic
istart += blocksize
iend += blocksize
x = CalibrationErrors.unsafe_unbiasedskce(kernel, predictions, targets, istart,
iend)

# update the estimate
Δestimate = x - estimate
estimate += Δestimate / b
S += Δestimate * (x - estimate)
end

# compute standard error and z-statistic
stderr = sqrt(S) / nblocks
z = estimate / stderr

return AsymptoticBlockSKCETest(kernel, blocksize, nblocks, estimate, stderr, z)
end

# HypothesisTests interface

HypothesisTests.default_tail(::AsymptoticBlockSKCETest) = :right

## have to specify and check keyword arguments in `pvalue` and `confint` to
## force `tail = :right` due to the default implementation in HypothesisTests

function HypothesisTests.pvalue(test::AsymptoticBlockSKCETest; tail = :right)
if tail === :right
normccdf(test.z)
else
throw(ArgumentError("tail=$(tail) is invalid"))
end
end

# confidence interval by inversion
function StatsBase.confint(test::AsymptoticBlockSKCETest; level = 0.95, tail = :right)
HypothesisTests.check_level(level)

if tail === :right
q = norminvcdf(level)
lowerbound = test.estimate - q * test.stderr
(lowerbound, oftype(lowerbound, Inf))
else
throw(ArgumentError("tail = $(tail) is invalid"))
end
end

HypothesisTests.testname(test::AsymptoticBlockSKCETest) = "Asymptotic block SKCE test"

# parameter of interest: name, value under H0, point estimate
function HypothesisTests.population_param_of_interest(test::AsymptoticBlockSKCETest)
"SKCE", zero(test.estimate), test.estimate
end

function HypothesisTests.show_params(io::IO, test::AsymptoticBlockSKCETest, ident = "")
println(io, ident, "number of observations per block: ", test.blocksize)
println(io, ident, "number of blocks of observations: ", test.nblocks)
println(io, ident, "z-statistic: ", test.z)
println(io, ident, "standard error of evaluations of pairs of observations: ",
test.stderr)
end
100 changes: 0 additions & 100 deletions src/skce/asymptotic_linear.jl

This file was deleted.

23 changes: 11 additions & 12 deletions src/skce/distribution_free.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
struct DistributionFreeTest{E<:SKCE,B,V} <: HypothesisTests.HypothesisTest
struct DistributionFreeSKCETest{E<:SKCE,B,V} <: HypothesisTests.HypothesisTest
"""Calibration estimator."""
estimator::E
"""Uniform upper bound of the terms of the estimator."""
Expand All @@ -9,21 +9,21 @@ struct DistributionFreeTest{E<:SKCE,B,V} <: HypothesisTests.HypothesisTest
estimate::V
end

function DistributionFreeTest(estimator::SKCE, data...; bound = uniformbound(estimator))
function DistributionFreeSKCETest(estimator::SKCE, data...; bound = uniformbound(estimator))
# obtain the predictions and targets
predictions, targets = CalibrationErrors.predictions_targets(data...)

# compute the calibration error estimate
estimate = calibrationerror(estimator, predictions, targets)

DistributionFreeTest(estimator, bound, length(predictions), estimate)
DistributionFreeSKCETest(estimator, bound, length(predictions), estimate)
end

# HypothesisTests interface

HypothesisTests.default_tail(::DistributionFreeTest) = :right
HypothesisTests.default_tail(::DistributionFreeSKCETest) = :right

function HypothesisTests.pvalue(test::DistributionFreeTest{<:BiasedSKCE})
function HypothesisTests.pvalue(test::DistributionFreeSKCETest{<:BiasedSKCE})
@unpack bound, n, estimate = test

s = sqrt(n * estimate / bound) - 1
Expand All @@ -34,21 +34,21 @@ function HypothesisTests.pvalue(test::DistributionFreeTest{<:BiasedSKCE})
exp(- s^2 / 2)
end

function HypothesisTests.pvalue(test::DistributionFreeTest{<:Union{QuadraticUnbiasedSKCE,
LinearUnbiasedSKCE}})
function HypothesisTests.pvalue(test::DistributionFreeSKCETest{<:Union{UnbiasedSKCE,
BlockUnbiasedSKCE}})
@unpack bound, n, estimate = test

exp(- div(n, 2) * estimate^2 / (2 * bound ^ 2))
end

HypothesisTests.testname(::DistributionFreeTest) = "Distribution-free calibration test"
HypothesisTests.testname(::DistributionFreeSKCETest) = "Distribution-free SKCE test"

# parameter of interest: name, value under H0, point estimate
function HypothesisTests.population_param_of_interest(test::DistributionFreeTest)
function HypothesisTests.population_param_of_interest(test::DistributionFreeSKCETest)
nameof(typeof(test.estimator)), zero(test.estimate), test.estimate
end

function HypothesisTests.show_params(io::IO, test::DistributionFreeTest, ident = "")
function HypothesisTests.show_params(io::IO, test::DistributionFreeSKCETest, ident = "")
println(io, ident, "number of observations: $(test.n)")
println(io, ident, "uniform bound of the terms of the estimator: $(test.bound)")
end
Expand All @@ -70,5 +70,4 @@ uniformbound(kernel::ScaledKernel) = first(kernel.σ²) * uniformbound(kernel.ke
uniformbound(kernel::TransformedKernel) = uniformbound(kernel.kernel)

# uniform bounds of the norm of tensor product kernels
uniformbound(kernel::TensorProductKernel) =
uniformbound(kernel.kernel1) * uniformbound(kernel.kernel2)
uniformbound(kernel::TensorProduct) = prod(uniformbound, kernel.kernels)
Loading

0 comments on commit b621b3b

Please sign in to comment.