Add Weibull distribution (#25)

* add Weibull

* return error if scaling is required

* error update

* bug fix

Author: guilhermebodin

src/ScoreDrivenModels.jl

@@ -27,5 +27,6 @@ include("distributions/poisson.jl")
 include("distributions/beta.jl")
 include("distributions/log_normal.jl")
 include("distributions/gamma.jl")
+include("distributions/weibull.jl")
 
 end

src/distributions/common_interface.jl

@@ -4,7 +4,8 @@ const DISTS = [
     Beta;
     Poisson;
     LogNormal;
-    Gamma
+    Gamma;
+    Weibull
 ]
 
 function score(y::T, D::Type{<:Distribution}, param::Vector{T}) where T

src/distributions/poisson.jl

@@ -16,7 +16,7 @@ end
 Proof somewhere
 """
 function log_likelihood(::Type{Poisson}, y::Vector{Int}, param::Vector{Vector{T}}, n::Int) where T
-    loglik = 0.0
+    loglik = zero(T)
     for i in 1:n
         loglik += y[i]*log(param[i][1]) - param[i][1] - lfactorial(y[i])
     end

src/distributions/weibull.jl

@@ -0,0 +1,58 @@
+"""
+Proof somewhere 
+parametrized in \\alpha and \\theta
+"""
+function score(y::T, ::Type{Weibull}, param::Vector{T}) where T
+    return [
+        (1/param[1]) + log(y/param[2]) * (1 - (y/param[2])^param[1]) ;
+        (param[1]/param[2]) * (((y/param[2])^param[1]) - 1)
+    ]
+end
+
+"""
+Proof somewhere
+"""
+function fisher_information(::Type{Weibull}, param::Vector{T}) where T
+    return error("Fisher information not implemented for Weibull distribution.")
+end
+
+"""
+Proof somewhere
+"""
+function log_likelihood(::Type{Weibull}, y::Vector{T}, param::Vector{Vector{T}}, n::Int) where T
+    loglik = zero(T)
+    for i in 1:n
+        loglik += log(param[i][1]) + (param[i][1] - 1) * log(y[i]) - param[i][1] * log(param[i][2]) - (y[i]/param[i][2])^param[i][1]
+    end
+    return -loglik
+end
+
+# Links
+function param_to_param_tilde(::Type{Weibull}, param::Vector{T}) where T 
+    return [
+        param_to_param_tilde(ExponentialLink, param[1], zero(T));
+        param_to_param_tilde(ExponentialLink, param[2], zero(T))
+    ]
+end
+function param_tilde_to_param(::Type{Weibull}, param_tilde::Vector{T}) where T 
+    return [
+        param_tilde_to_param(ExponentialLink, param_tilde[1], zero(T));
+        param_tilde_to_param(ExponentialLink, param_tilde[2], zero(T))
+    ]
+end
+function jacobian_param_tilde(::Type{Weibull}, param_tilde::Vector{T}) where T 
+    return Diagonal([
+        jacobian_param_tilde(ExponentialLink, param_tilde[1], zero(T));
+        jacobian_param_tilde(ExponentialLink, param_tilde[2], zero(T))
+    ])
+end
+
+# utils 
+function update_dist(::Type{Weibull}, param::Vector{T}) where T
+    # normal here is parametrized as sigma^2
+    return Weibull(param[1], param[2])
+end 
+
+function num_params(::Type{Weibull})
+    return 2
+end

src/link_functions.jl

@@ -10,7 +10,11 @@ struct ExponentialLink <: Link end
 
 param_to_param_tilde(::Type{ExponentialLink}, param::T, lower_bound::T) where T = log(param - lower_bound)
 param_tilde_to_param(::Type{ExponentialLink}, param_tilde::T, lower_bound::T) where T = exp(param_tilde) + lower_bound
-jacobian_param_tilde(::Type{ExponentialLink}, param_tilde::T, lower_bound::T) where T = exp(param_tilde) + lower_bound
+jacobian_param_tilde(::Type{ExponentialLink}, param_tilde::T, lower_bound::T) where T = exp(param_tilde)
 
 struct LogitLink <: Link end
 
+const LINKS = [
+    IdentityLink;
+    ExponentialLink
+]

test/runtests.jl

@@ -5,7 +5,8 @@ using Test, Random, HypothesisTests
 const SDM = ScoreDrivenModels
 
 include("test_recursion.jl")
-include("test_estimate.jl")
+include("test_links.jl")
 include("test_distributions.jl")
 include("test_initial_params.jl")
-include("test_diagnostics.jl")
+include("test_diagnostics.jl")
+include("test_estimate.jl")

test/test_links.jl

@@ -0,0 +1,15 @@
+@testset "Links" begin
+    Random.seed!(10)
+    x = rand()
+    atol = 1e-7 
+    rtol = 1e-7
+
+    link = SDM.IdentityLink
+    @test SDM.param_to_param_tilde(link, SDM.param_tilde_to_param(link, x)) ≈ x atol = atol rtol = rtol
+
+    link = SDM.ExponentialLink
+    lb = 1.0
+    @test SDM.param_to_param_tilde(link, SDM.param_tilde_to_param(link, x, lb), lb) ≈ x atol = atol rtol = rtol
+    lb = 0.0
+    @test SDM.param_to_param_tilde(link, SDM.param_tilde_to_param(link, x, lb), lb) ≈ x atol = atol rtol = rtol
+end