Add beta location scale and dont cheat at tests (#89)

* add beta location scale and dont cheat at tests

* add test for params_sdm

Author: guilhermebodin

src/MLE.jl

@@ -117,15 +117,15 @@ function fit(gas::Model{D, T}, y::Vector{T};
             func = TwiceDifferentiable(psi_tilde -> log_lik(psi_tilde, y, gas_fit, initial_params, unknowns, n), opt_method.initial_points[i])
             opt_result = optimize(func, opt_method, verbose, i)
             update_aux_estimation!(aux_est, func, opt_result)
-            println("initial_point $i of $n_initial_points - $(-opt_result.minimum)")
+            println("initial point $i of $n_initial_points - Log-likelihood: $(-opt_result.minimum)")
         catch err
             println(err)
-            println("initial_point $i diverged")
+            println("initial point $i diverged")
         end
     end
 
     if isempty(aux_est.loglikelihood) 
-        println("No initial_point converged.")
+        println("No initial point converged.")
         return
     end
 
@@ -146,7 +146,7 @@ end
 
 function fit!(gas::Model{D, T}, y::Vector{T};
               initial_params::Matrix{T} = DEFAULT_INITIAL_PARAM,
-              opt_method::AbstractOptimizationMethod = LBFGS(gas, DEFAULT_NUM_SEEDS),
+              opt_method::AbstractOptimizationMethod = NelderMead(gas, DEFAULT_NUM_SEEDS),
               verbose::Int = DEFAULT_VERBOSE) where {D, T}
 
     unknowns = find_unknowns(gas)

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/betalocationscale.jl")
 include("distributions/chi.jl")
 include("distributions/chisq.jl")
 include("distributions/exponential.jl")

src/distributions/beta.jl

@@ -59,6 +59,10 @@ function update_dist(::Type{Beta}, param::Matrix{T}, t::Int) where T
     return Beta(param[t, 1], param[t, 2])
 end 
 
+function params_sdm(d::Beta)
+    return Distributions.params(d)
+end
+
 function num_params(::Type{Beta})
     return 2
 end

src/distributions/betalocationscale.jl

@@ -0,0 +1,73 @@
+"""
+    BetaLocationScale
+
+* Parametrization \\mu, \\sigma, \\alpha, \\beta
+
+* Score
+
+* Fisher Information
+
+* `time_varying_params` map.
+
+* Default link
+"""
+BetaLocationScale
+
+
+# Remember for all the code that param[t, 2] is c - a
+# and param[t, 1] is a so to get c one must make 
+# param[t, 2] + param[t, 1]
+function score!(score_til::Matrix{T}, y::T, ::Type{BetaLocationScale}, 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 fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{BetaLocationScale}, param::Matrix{T}, t::Int) where T
+    return error("Fisher information not implemented for BetaLocationScale distribution.")
+end
+
+function log_likelihood(::Type{BetaLocationScale}, 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{BetaLocationScale}, param::Matrix{T}, t::Int) where T 
+    param_tilde[t, 1] = link(LogLink, param[t, 1], zero(T))
+    param_tilde[t, 2] = link(LogLink, param[t, 2], zero(T))
+    return
+end
+function unlink!(param::Matrix{T}, ::Type{BetaLocationScale}, param_tilde::Matrix{T}, t::Int) where T 
+    param[t, 1] = unlink(LogLink, param_tilde[t, 1], zero(T))
+    param[t, 2] = unlink(LogLink, param_tilde[t, 2], zero(T))
+    return
+end
+function jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{BetaLocationScale}, param::Matrix{T}, t::Int) where T 
+    aux.jac[1] = jacobian_link(LogLink, param[t, 1], zero(T))
+    aux.jac[2] = jacobian_link(LogLink, param[t, 2], zero(T))
+    return
+end
+
+# utils
+function update_dist(::Type{BetaLocationScale}, 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::BetaLocationScale)
+    return (d.μ, d.σ + d.μ, Distributions.params(d.ρ)...)
+end
+
+function num_params(::Type{BetaLocationScale})
+    return 4
+end

src/distributions/chi.jl

@@ -52,6 +52,10 @@ function update_dist(::Type{Chi}, param::Matrix{T}, t::Int) where T
     return Chi(param[t, 1])
 end 
 
+function params_sdm(d::Chi)
+    return Distributions.params(d)
+end
+
 function num_params(::Type{Chi})
     return 1
 end

src/distributions/chisq.jl

@@ -52,6 +52,10 @@ function update_dist(::Type{Chisq}, param::Matrix{T}, t::Int) where T
     return Chisq(param[t, 1])
 end 
 
+function params_sdm(d::Chisq)
+    return Distributions.params(d)
+end
+
 function num_params(::Type{Chisq})
     return 1
 end

src/distributions/common_interface.jl

@@ -1,6 +1,7 @@
 # Currently supported distributions
 const DISTS = [
     Beta,
+    BetaLocationScale,
     Chi,
     Chisq,
     Exponential,
@@ -15,6 +16,7 @@ const DISTS = [
 ]
 
 export Beta,
+    BetaLocationScale,
     Chi,
     Chisq,
     Exponential,
@@ -27,35 +29,59 @@ export Beta,
     TDistLocationScale,
     Weibull
 
-# Query params for generic distribution
-params_sdm(d::Distribution) = Distributions.params(d)
-# Query params for all LocationScale distributions
-params_sdm(d::LocationScale) = (d.μ, d.σ, Distributions.params(d.ρ)...)
-
 function score!(score_til::Matrix{T}, y::T, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T
     return error("score! not implemented for $D distribution")
 end
+
 function score!(score_til::Matrix{T}, y::Int, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T
     return error("score! not implemented for $D distribution")
 end
+
 function fisher_information!(aux::AuxiliaryLinAlg{T}, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T
     return error("fisher_information! not implemented for $D distribution")
 end
+
 function log_likelihood(D::Type{<:Distribution}, y::Vector{T}, param::Matrix{T}, n::Int) where T
     return error("log_likelihood not implemented for $D distribution")
 end
+
 function link!(param_tilde::Matrix{T}, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T 
     return error("link! not implemented for $D distribution")
 end
+
 function unlink!(param::Matrix{T}, D::Type{<:Distribution}, param_tilde::Matrix{T}, t::Int) where T 
     return error("unlink! not implemented for $D distribution")
 end
+
 function jacobian_link!(aux::AuxiliaryLinAlg{T}, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T 
     return error("jacobian_link! not implemented for $D distribution")
 end
+
+"""
+    update_dist(D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T
+
+Create a new distribution from Distributions.jl based on the parametrization used in
+ScoreDrivenModels.jl.
+"""
 function update_dist(D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T
     return error("update_dist not implemented for $D distribution")
+end
+
+"""
+    params_sdm(d::Distribution)
+
+Recover the parametrization used in ScoreDrivenModels.jl based on a Distribution from
+Distributions.jl.
+"""
+function params_sdm(d::Distribution)
+    return error("params_sdm not implemented for $d distribution")
 end 
+
+"""
+    num_params(D::Type{<:Distribution})
+
+Number of parameters of a given distribution.
+"""
 function num_params(D::Type{<:Distribution})
     return error("num_params not implemented for $D distribution")
 end

src/distributions/exponential.jl

@@ -52,6 +52,11 @@ function update_dist(::Type{Exponential}, param::Matrix{T}, t::Int) where T
     return Exponential(1/param[t, 1])
 end 
 
+function params_sdm(d::Exponential)
+    pars = Distributions.params(d)
+    return 1/pars[1]
+end
+
 function num_params(::Type{Exponential})
     return 1
 end

src/distributions/gamma.jl

@@ -59,6 +59,10 @@ function update_dist(::Type{Gamma}, param::Matrix{T}, t::Int) where T
     return Gamma(param[t, 1], param[t, 2])
 end 
 
+function params_sdm(d::Gamma)
+    return Distributions.params(d)
+end
+
 function num_params(::Type{Gamma})
     return 2
 end

src/distributions/logitnormal.jl

@@ -61,6 +61,11 @@ function update_dist(::Type{LogitNormal}, param::Matrix{T}, t::Int) where T
     return LogitNormal(param[t, 1], sqrt(param[t, 2]))
 end 
 
+function params_sdm(d::LogitNormal)
+    pars = Distributions.params(d)
+    return (pars[1], pars[2]^2)
+end
+
 function num_params(::Type{LogitNormal})
     return 2
 end