Add BetaFourParameters (#121)

Author: guilhermebodin

docs/src/manual.md

@@ -95,6 +95,7 @@ this is handled internally.
 | Distribution | Identity scaling | Inverse and inverse square root scalings |
 | :------- |:---------------:|:--------------------:|
 |[`Beta`](@ref)| ✓ | ✓  |
+|[`BetaFourParameters`](@ref)| ✓ |  x  |
 |[`Exponential`](@ref)| ✓ | ✓  |
 |[`Gamma`](@ref)| ✓ | ✓  |
 |[`LogitNormal`](@ref)| ✓ | ✓  |
@@ -109,6 +110,7 @@ this is handled internally.
 
 ```@docs
 ScoreDrivenModels.Beta
+ScoreDrivenModels.BetaFourParameters
 ScoreDrivenModels.Exponential
 ScoreDrivenModels.Gamma
 ScoreDrivenModels.LogitNormal

src/MLE.jl

@@ -155,14 +155,15 @@ end
 function fit!(gas::Model{D, T}, y::Vector{T};
               initial_params::Matrix{T} = DEFAULT_INITIAL_PARAM,
               opt_method::AbstractOptimizationMethod = NelderMead(gas, DEFAULT_NUM_SEEDS),
-              verbose::Int = DEFAULT_VERBOSE) where {D, T}
+              verbose::Int = DEFAULT_VERBOSE,
+              throw_errors::Bool = false) where {D, T}
 
     unknowns = find_unknowns(gas)
     # Check if the model has no unknowns
     n_unknowns = length(unknowns)
     check_model_estimated(n_unknowns) && return gas
 
-    f = fit(gas, y; initial_params = initial_params, opt_method = opt_method, verbose = verbose)
+    f = fit(gas, y; initial_params = initial_params, opt_method = opt_method, verbose = verbose, throw_errors = throw_errors)
     fill_psitilde!(gas, f.coefs, unknowns)
     return f
 end

src/ScoreDrivenModels.jl

@@ -32,6 +32,7 @@ include("opt_methods/NelderMead.jl")
 include("distributions/non_native_dists.jl")
 include("distributions/common_interface.jl")
 include("distributions/beta.jl")
+include("distributions/betafour.jl")
 include("distributions/exponential.jl")
 include("distributions/gamma.jl")
 include("distributions/logitnormal.jl")

src/distributions/betafour.jl

@@ -0,0 +1,86 @@
+"""
+    BetaFourParameters
+* Parametrization a, c, \\alpha, \\beta
+* Score
+* Fisher Information
+* `time_varying_params` map.
+* Default link
+
+Right now when estimating a BetaFourParameters model is recomended to provide fixed parameters a and c, dont 
+with the code:
+
+```julia
+gas_beta_4 = Model([1, 2, 11, 12], [1, 2, 11, 12], BetaFourParameters, 0.0; time_varying_params=[3])
+gas_beta_4.ω[1] = minimum(y) - 0.1*std(y) # parameter a
+gas_beta_4.ω[2] = maximum(y) + 0.1*std(y) # parameter c
+```
+
+This code is a simple and not very accurate heuristic on how to estimate the a and c parameters. If you have estimated 
+using maximum likelihood it will be probably better
+
+The fact occurs because if in the optimization the gradient makes `c < maximum(y)` or `a > minimum(y)` then it 
+is impossible that y comes from the BetaFourParameters distribution leading to DomainErrors
+"""
+BetaFourParameters
+
+function score!(score_til::Matrix{T}, y::T, ::Type{BetaFourParameters}, param::Matrix{T}, t::Int) where T
+    score_til[t, 1] = -(param[t, 3] - 1)/(y - param[t, 1]) + (param[t, 4] + param[t, 3] - 1)/(param[t, 2] - param[t, 1])
+    score_til[t, 2] = (param[t, 4] - 1)/(param[t, 2] - y) - (param[t, 4] + param[t, 3] - 1)/(param[t, 2] - param[t, 1])
+    score_til[t, 3] = log(y - param[t, 1]) - log(param[t, 2] - param[t, 1]) + 
+                        digamma(param[t, 3] + param[t, 4]) - digamma(param[t, 3])
+    score_til[t, 4] = log(param[t, 2] - y) - log(param[t, 2] - param[t, 1]) + 
+                        digamma(param[t, 3] + param[t, 4]) - digamma(param[t, 4])
+    return
+end
+
+function log_likelihood(::Type{BetaFourParameters}, y::Vector{T}, param::Matrix{T}, n::Int) where T
+    loglik = 0.0
+    for t in 1:n
+        loglik += (param[t, 3] - 1)*log(y[t] - param[t, 1]) + (param[t, 4] - 1) * log(param[t, 2] - y[t]) - 
+                  (param[t, 3] + param[t, 4] - 1) * log(param[t, 2] - param[t, 1]) - logbeta(param[t, 3], param[t, 4])
+    end
+    return -loglik
+end
+
+# Links
+function link!(param_tilde::Matrix{T}, ::Type{BetaFourParameters}, param::Matrix{T}, t::Int) where T 
+    param_tilde[t, 1] = link(IdentityLink, param[t, 1])
+    param_tilde[t, 2] = link(IdentityLink, param[t, 2])
+    param_tilde[t, 3] = link(LogLink, param[t, 3], zero(T))
+    param_tilde[t, 4] = link(LogLink, param[t, 4], zero(T))
+    return
+end
+function unlink!(param::Matrix{T}, ::Type{BetaFourParameters}, param_tilde::Matrix{T}, t::Int) where T 
+    param[t, 1] = unlink(IdentityLink, param_tilde[t, 1])
+    param[t, 2] = unlink(IdentityLink, param_tilde[t, 2])
+    param[t, 3] = unlink(LogLink, param_tilde[t, 3], zero(T))
+    param[t, 4] = unlink(LogLink, param_tilde[t, 4], zero(T))
+    return
+end
+function jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{BetaFourParameters}, param::Matrix{T}, t::Int) where T 
+    aux.jac[1] = jacobian_link(IdentityLink, param[t, 1])
+    aux.jac[2] = jacobian_link(IdentityLink, param[t, 2])
+    aux.jac[3] = jacobian_link(LogLink, param[t, 3], zero(T))
+    aux.jac[4] = jacobian_link(LogLink, param[t, 4], zero(T))
+    return
+end
+
+# utils
+function update_dist(::Type{BetaFourParameters}, param::Matrix{T}, t::Int) where T
+    small_threshold!(param[:, 3:4], SMALL_NUM, t)
+    beta = Beta(param[t, 3], param[t, 4])
+    return LocationScale(param[t, 1], param[t, 2] - param[t, 1], beta)
+end 
+
+function params_sdm(d::BetaFourParameters)
+    return (d.μ, d.σ + d.μ, Distributions.params(d.ρ)...)
+end
+
+function num_params(::Type{BetaFourParameters})
+    return 4
+end
+
+# d = Beta(1, 1)
+# d2 = LocationScale(-1, 3, d)
+# maximum(d2)
+# minimum(d2)

src/distributions/common_interface.jl

@@ -1,6 +1,7 @@
 # Currently supported distributions
 const DISTS = [
     Beta,
+    BetaFourParameters,
     Exponential,
     Gamma,
     LogitNormal,
@@ -14,6 +15,7 @@ const DISTS = [
 ]
 
 export Beta,
+    BetaFourParameters,
     Exponential,
     Gamma,
     LogitNormal,

src/distributions/non_native_dists.jl

@@ -1,3 +1,6 @@
 # Define TDistLocationScale as the location scale transformation of the Distributions.jl TDist
 export TDistLocationScale
-TDistLocationScale = LocationScale{Float64,TDist{Float64}}
+TDistLocationScale = LocationScale{Float64,TDist{Float64}}
+
+export BetaFourParameters
+BetaFourParameters = LocationScale{Float64,Beta{Float64}}

test/utils.jl

@@ -1,17 +1,19 @@
-const PARAMETER_REQUIRED_DISTS = [Chisq; Chi; TDist]
-const FI_NOT_IMPLEMENTED = [Weibull; NegativeBinomial]
+const FI_NOT_IMPLEMENTED = [Weibull; BetaFourParameters; NegativeBinomial]
 
 struct FakeDist{T<:Real} <: Distributions.ContinuousUnivariateDistribution
     foo::T
 end
 
 function instantiate_dist(D::Type{<:Distribution})
-    if D == TDistLocationScale
+    if D == ScoreDrivenModels.TDistLocationScale
         params = [1.0 2.0 10]
         return ScoreDrivenModels.update_dist(D, params, 1)
-    elseif D in PARAMETER_REQUIRED_DISTS
+    elseif D == ScoreDrivenModels.TDist
         params = [10][:, :]
         return ScoreDrivenModels.update_dist(D, params, 1)
+    elseif D == ScoreDrivenModels.BetaFourParameters
+        params = [0.0 10.0 2.0 3]
+        return ScoreDrivenModels.update_dist(D, params, 1)
     else
         params = 0.5 * ones(1, ScoreDrivenModels.num_params(D))
         return ScoreDrivenModels.update_dist(D, params, 1)
@@ -26,8 +28,8 @@ function test_score_mean(D::Type{<:Distribution}; n::Int = 10^7, seed::Int = 10,
     n_params = ScoreDrivenModels.num_params(D)
     avg  = zeros(1, n_params)
     obs = rand(dist, n)
+    score_aux = ones(1, n_params)
     for i = 1:n
-        score_aux = ones(1, n_params)
         ScoreDrivenModels.score!(score_aux, obs[i], D, pars, 1)
         avg .+= score_aux
     end
@@ -43,8 +45,8 @@ function test_fisher_information(D::Type{<:Distribution}; n::Int = 10^6, seed::I
     n_params = ScoreDrivenModels.num_params(D)
     var_terms  = zeros(n_params, n_params)
     obs = rand(dist, n)
+    score_aux = ones(1, n_params)
     for i = 1:n
-        score_aux = ones(1, n_params)
         ScoreDrivenModels.score!(score_aux, obs[i], D, pars, 1)
         var_terms .+= score_aux' * score_aux
     end