Add TDist (#88)

guilhermebodin · web-flow · commit f8d7c90c4983 · 2020-02-28T15:16:53.000-03:00
* add TDist

* change location scale name order
diff --git a/docs/src/manual.md b/docs/src/manual.md
@@ -94,11 +94,12 @@ ScoreDrivenModels.Chi
 ScoreDrivenModels.Chisq
 ScoreDrivenModels.Exponential
 ScoreDrivenModels.Gamma
-ScoreDrivenModels.LocationScaleTDist
 ScoreDrivenModels.LogitNormal
 ScoreDrivenModels.LogNormal
 ScoreDrivenModels.Normal
 ScoreDrivenModels.Poisson
+ScoreDrivenModels.TDist
+ScoreDrivenModels.TDistLocationScale
 ScoreDrivenModels.Weibull
 ```
 
diff --git a/examples/cpichg.jl b/examples/cpichg.jl
@@ -3,7 +3,7 @@ using ScoreDrivenModels, DelimitedFiles
 y = readdlm("./test/data/cpichg.csv")[:]
 
 # Define a model where only \mu and \sigma^2 is varying over time
-gas = Model(1, 1, LocationScaleTDist, 0.0, time_varying_params = [1; 2])
+gas = Model(1, 1, TDistLocationScale, 0.0, time_varying_params = [1; 2])
 
 # Estimate the model
 @time f = fit!(gas, y)
diff --git a/src/ScoreDrivenModels.jl b/src/ScoreDrivenModels.jl
@@ -40,7 +40,8 @@ include("distributions/logitnormal.jl")
 include("distributions/lognormal.jl")
 include("distributions/normal.jl")
 include("distributions/poisson.jl")
-include("distributions/locationscalestd.jl")
+include("distributions/tdist.jl")
+include("distributions/tdistlocationscale.jl")
 include("distributions/weibull.jl")
 
 end
diff --git a/src/distributions/beta.jl b/src/distributions/beta.jl
@@ -11,7 +11,7 @@
 
 * Default link
 """
-function Beta end
+Beta
 
 function score!(score_til::Matrix{T}, y::T, ::Type{Beta}, param::Matrix{T}, t::Int) where T
     score_til[t, 1] = log(y) + digamma(param[t, 1] + param[t, 2]) - digamma(param[t, 1])
diff --git a/src/distributions/chi.jl b/src/distributions/chi.jl
@@ -12,7 +12,7 @@ parametrized in k
 
 * Default link
 """
-function Chi end
+Chi
 
 function score!(score_til::Matrix{T}, y::T, ::Type{Chi}, param::Matrix{T}, t::Int) where T
     score_til[t, 1] = log(y) - log(2)/2 - digamma(param[t, 1]/2)/2
diff --git a/src/distributions/chisq.jl b/src/distributions/chisq.jl
@@ -12,7 +12,7 @@ parametrized in k
 
 * Default link
 """
-function Chisq end
+Chisq
 
 function score!(score_til::Matrix{T}, y::T, ::Type{Chisq}, param::Matrix{T}, t::Int) where T
     score_til[t, 1] = (log(y) - log(2) - digamma(param[t, 1]/2))/2
diff --git a/src/distributions/common_interface.jl b/src/distributions/common_interface.jl
@@ -9,7 +9,8 @@ const DISTS = [
     LogNormal,
     Normal,
     Poisson,
-    LocationScaleTDist,
+    TDist,
+    TDistLocationScale,
     Weibull
 ]
 
@@ -22,7 +23,8 @@ export Beta,
     LogNormal,
     Normal,
     Poisson,
-    LocationScaleTDist,
+    TDist,
+    TDistLocationScale,
     Weibull
 
 # Query params for generic distribution
diff --git a/src/distributions/exponential.jl b/src/distributions/exponential.jl
@@ -12,7 +12,7 @@ parametrized in \\lambda
 
 * Default link
 """
-function Exponential end
+Exponential
 
 function score!(score_til::Matrix{T}, y::T, ::Type{Exponential}, param::Matrix{T}, t::Int) where T
     score_til[t, 1] = 1/param[t, 1] - y
diff --git a/src/distributions/gamma.jl b/src/distributions/gamma.jl
@@ -12,7 +12,7 @@ parametrized in \\alpha and \\theta
 
 * Default link
 """
-function Gamma end
+Gamma
 
 function score!(score_til::Matrix{T}, y::T, ::Type{Gamma}, param::Matrix{T}, t::Int) where T
     score_til[t, 1] = log(y) - digamma(param[t, 1]) - log(param[t, 2])
diff --git a/src/distributions/logitnormal.jl b/src/distributions/logitnormal.jl
@@ -12,7 +12,7 @@ parametrized in \\mu and \\sigma^2
 
 * Default link
 """
-function LogitNormal end
+LogitNormal
 
 logit(x) = log(x) - log(1-x)
 
diff --git a/src/distributions/lognormal.jl b/src/distributions/lognormal.jl
@@ -12,7 +12,7 @@ parametrized in \\mu and \\sigma^2
 
 * Default link
 """
-function LogNormal end
+LogNormal
 
 function score!(score_til::Matrix{T}, y::T, ::Type{LogNormal}, param::Matrix{T}, t::Int) where T
     score_til[t, 1] = (log(y) - param[t, 1])/param[t, 2]
diff --git a/src/distributions/non_native_dists.jl b/src/distributions/non_native_dists.jl
@@ -1,3 +1,3 @@
 # Define Std as the location scale transformation of the Distributions.jl TDist
-export LocationScaleTDist
-LocationScaleTDist = LocationScale{Float64,TDist{Float64}}
+export TDistLocationScale
+TDistLocationScale = LocationScale{Float64,TDist{Float64}}
diff --git a/src/distributions/normal.jl b/src/distributions/normal.jl
@@ -12,7 +12,7 @@ parametrized in \\mu and \\sigma^2
 
 * Default link
 """
-function Normal end
+Normal
 
 function score!(score_til::Matrix{T}, y::T, ::Type{Normal}, param::Matrix{T}, t::Int) where T
     score_til[t, 1] = (y - param[t, 1])/param[t, 2]
diff --git a/src/distributions/poisson.jl b/src/distributions/poisson.jl
@@ -12,7 +12,7 @@ parametrized in \\lambda
 
 * Default link
 """
-function Poisson end
+Poisson
 
 function score!(score_til::Matrix{T}, y::Int, ::Type{Poisson}, param::Matrix{T}, t::Int) where T
     score_til[t, 1] = y/param[t, 1] - 1
diff --git a/src/distributions/tdist.jl b/src/distributions/tdist.jl
@@ -0,0 +1,61 @@
+"""
+    Student's t
+
+* Parametrization
+parametrized in \\nu
+
+* Score
+
+* Fisher Information
+
+* `time_varying_params` map.
+
+* Default link
+"""
+TDist
+
+function score!(score_til::Matrix{T}, y::T, ::Type{TDist}, param::Matrix{T}, t::Int) where T
+    score_til[t, 1] = 0.5 * (((y^2)*(param[t, 1] + 1)/(param[t, 1] * y^2 + param[t, 1]^2)) -
+                      1/param[t, 1] - 
+                      log((y^2)/(param[t, 1]) + 1) - 
+                      digamma(param[t, 1] / 2) + digamma((param[t, 1] + 1) / 2))
+    return
+end
+
+function fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{TDist}, param::Matrix{T}, t::Int) where T
+    aux.fisher[1, 1] = 0.5 * (0.5 * trigamma(0.5 * param[t, 1]) - 0.5 * trigamma( 0.5 * (param[t, 1] + 1.0) ) -
+                             (param[t, 1] + 5.0) / (param[t, 1] * (param[t, 1] + 3.0) * (param[t, 1] + 1.0)))
+    return
+end
+
+function log_likelihood(::Type{TDist}, y::Vector{T}, param::Matrix{T}, n::Int) where T
+    loglik = zero(T)
+    for t in 1:n
+        loglik -= 0.5 * log(param[t, 1]) + logbeta(0.5, param[t, 1]/2) + ((param[t, 1] + 1)/2) *
+                  log(1 + (y[t]^2)/param[t, 1])
+    end
+    return -loglik
+end
+
+# Links
+function link!(param_tilde::Matrix{T}, ::Type{TDist}, param::Matrix{T}, t::Int) where T 
+    param_tilde[t, 1] = link(LogLink, param[t, 1], zero(T))
+    return
+end
+function unlink!(param::Matrix{T}, ::Type{TDist}, param_tilde::Matrix{T}, t::Int) where T 
+    param[t, 1] = unlink(LogLink, param_tilde[t, 1], zero(T))
+    return
+end
+function jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{TDist}, param::Matrix{T}, t::Int) where T 
+    aux.jac[1] = jacobian_link(LogLink, param[t, 1], zero(T))
+    return
+end
+
+# utils 
+function update_dist(::Type{TDist}, param::Matrix{T}, t::Int) where T
+    return TDist(param[t, 1])
+end 
+
+function num_params(::Type{TDist})
+    return 1
+end
diff --git a/src/distributions/tdistlocationscale.jl b/src/distributions/tdistlocationscale.jl
@@ -1,5 +1,5 @@
 """
-    Location Scale Student's t
+    Student's t location scale
 
 * Parametrization
 parametrized in \\mu, \\sigma^2 and \\nu
@@ -12,9 +12,9 @@ parametrized in \\mu, \\sigma^2 and \\nu
 
 * Default link
 """
-LocationScaleTDist
+TDistLocationScale
 
-function score!(score_til::Matrix{T}, y::T, ::Type{LocationScaleTDist}, param::Matrix{T}, t::Int) where T
+function score!(score_til::Matrix{T}, y::T, ::Type{TDistLocationScale}, param::Matrix{T}, t::Int) where T
     score_til[t, 1] = ((param[t, 3] + 1) * (y - param[t, 1])) / ((y - param[t, 1])^2 + param[t, 2] * param[t, 3])
     score_til[t, 2] = -(param[t, 3] * (param[t, 2] - (y - param[t, 1])^2)) / (2 * param[t, 2] * (param[t, 3] * param[t, 2] + 
                         (y - param[t, 1])^2))
@@ -25,7 +25,7 @@ function score!(score_til::Matrix{T}, y::T, ::Type{LocationScaleTDist}, param::M
     return
 end
 
-function fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{LocationScaleTDist}, param::Matrix{T}, t::Int) where T
+function fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{TDistLocationScale}, param::Matrix{T}, t::Int) where T
     aux.fisher[1, 1] = (param[t, 3] + 1.0)/(param[t, 2] * (param[t, 3] + 3.0))
     aux.fisher[1, 2] = zero(T)
     aux.fisher[2, 1] = zero(T)
@@ -39,7 +39,7 @@ function fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{LocationScaleTDist}
     return
 end
 
-function log_likelihood(::Type{LocationScaleTDist}, y::Vector{T}, param::Matrix{T}, n::Int) where T
+function log_likelihood(::Type{TDistLocationScale}, y::Vector{T}, param::Matrix{T}, n::Int) where T
     loglik = zero(T)
     for t in 1:n
         loglik -= 0.5 * log(param[t, 2]*param[t, 3]) + logbeta(0.5, param[t, 3]/2) + ((param[t, 3] + 1)/2) *
@@ -49,31 +49,31 @@ function log_likelihood(::Type{LocationScaleTDist}, y::Vector{T}, param::Matrix{
 end
 
 # Links
-function link!(param_tilde::Matrix{T}, ::Type{LocationScaleTDist}, param::Matrix{T}, t::Int) where T 
+function link!(param_tilde::Matrix{T}, ::Type{TDistLocationScale}, param::Matrix{T}, t::Int) where T 
     param_tilde[t, 1] = link(IdentityLink, param[t, 1])
     param_tilde[t, 2] = link(LogLink, param[t, 2], zero(T))
     param_tilde[t, 3] = link(LogLink, param[t, 3], zero(T))
     return
 end
-function unlink!(param::Matrix{T}, ::Type{LocationScaleTDist}, param_tilde::Matrix{T}, t::Int) where T 
+function unlink!(param::Matrix{T}, ::Type{TDistLocationScale}, param_tilde::Matrix{T}, t::Int) where T 
     param[t, 1] = unlink(IdentityLink, param_tilde[t, 1])
     param[t, 2] = unlink(LogLink, param_tilde[t, 2], zero(T))
     param[t, 3] = unlink(LogLink, param_tilde[t, 3], zero(T))
     return
 end
-function jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{LocationScaleTDist}, param::Matrix{T}, t::Int) where T 
+function jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{TDistLocationScale}, param::Matrix{T}, t::Int) where T 
     aux.jac[1] = jacobian_link(IdentityLink, param[t, 1])
     aux.jac[2] = jacobian_link(LogLink, param[t, 2], zero(T))
     aux.jac[3] = jacobian_link(LogLink, param[t, 3], zero(T))
     return
 end
 
 # utils 
-function update_dist(::Type{LocationScaleTDist}, param::Matrix{T}, t::Int) where T
+function update_dist(::Type{TDistLocationScale}, param::Matrix{T}, t::Int) where T
     tdist = TDist(param[t, 3])
     return LocationScale(param[t, 1], sqrt(param[t, 2]), tdist)
 end 
 
-function num_params(::Type{LocationScaleTDist})
+function num_params(::Type{TDistLocationScale})
     return 3
 end
diff --git a/src/distributions/weibull.jl b/src/distributions/weibull.jl
@@ -12,7 +12,7 @@ parametrized in \\alpha and \\theta
 
 * Default link
 """
-function Weibull end
+Weibull
 
 function score!(score_til::Matrix{T}, y::T, ::Type{Weibull}, param::Matrix{T}, t::Int) where T
     score_til[t, 1] = (1/param[t, 1]) + log(y/param[t, 2]) * (1 - (y/param[t, 2])^param[t, 1]) 
diff --git a/test/utils.jl b/test/utils.jl
@@ -1,11 +1,12 @@
-const PARAMETER_REQUIRED_DISTS = [Chisq; Chi]
+const PARAMETER_REQUIRED_DISTS = [Chisq; Chi; TDist]
+const FI_NOT_IMPLEMENTED = [Weibull]
 
 struct FakeDist{T<:Real} <: Distributions.ContinuousUnivariateDistribution
     foo::T
 end
 
 function instantiate_dist(D::Type{<:Distribution})
-    if D == LocationScaleTDist
+    if D == TDistLocationScale
         return D(1.0, 1.0, TDist(10))
     elseif D in PARAMETER_REQUIRED_DISTS
         return D(10)
@@ -15,7 +16,7 @@ function instantiate_dist(D::Type{<:Distribution})
 end
 
 function create_sampler(D::Type{<:Distribution}, pars)
-    if D == LocationScaleTDist
+    if D == TDistLocationScale
         return D(pars[1], pars[2], TDist(pars[3]))
     else
         return D(pars...)
@@ -56,10 +57,10 @@ function test_fisher_information(D::Type{<:Distribution}; n::Int = 10^6, seed::I
     aux_lin_alg = AuxiliaryLinAlg{Float64}(n_params)
 
     # Some distributions might not have the fisher information yet available
-    try 
+    if D in FI_NOT_IMPLEMENTED
+        return 
+    else
         ScoreDrivenModels.fisher_information!(aux_lin_alg, D, pars, 1)
-    catch 
-        return
     end
 
     @test var_terms ≈ aux_lin_alg.fisher atol = atol rtol = rtol
