fix scaling (#33)

guilhermebodin · raphaelsaavedra · web-flow · commit 4df2542072d5 · 2019-11-27T16:53:22.000-03:00
* fix scaling

* Update src/distributions/README.md

Co-Authored-By: Raphael Saavedra &lt;raphael.saavedra93@gmail.com&gt;
diff --git a/Project.toml b/Project.toml
@@ -9,16 +9,16 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
-[extras]
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
+[compat]
+Distributions = "~0.21"
+Optim = "~0.19"
+SpecialFunctions = "~0.8"
+julia = "~1"
 
 [targets]
 test = ["Test", "Random", "HypothesisTests"]
 
-[compat]
-julia = "~1"
-Distributions = "~0.21"
-Optim = "~0.19"
-SpecialFunctions = "~0.8"
+[extras]
+HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/src/distributions/README.md b/src/distributions/README.md
@@ -9,6 +9,6 @@ Each distribution must have the following methods:
 * link interface
     * param_to_param_tilde
     * param_tilde_to_param
-    * jacobian_param_tilde
+    * jacobian_link
 * update_dist
 * num_params
diff --git a/src/distributions/beta.jl b/src/distributions/beta.jl
@@ -46,15 +46,16 @@ function param_tilde_to_param(::Type{Beta}, param_tilde::Vector{T}) where T
         param_tilde_to_param(ExponentialLink, param_tilde[2], zero(T))
     ]
 end
-function jacobian_param_tilde(::Type{Beta}, param_tilde::Vector{T}) where T 
+function jacobian_link(::Type{Beta}, 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))
+        jacobian_link(ExponentialLink, param_tilde[1], zero(T));
+        jacobian_link(ExponentialLink, param_tilde[2], zero(T))
     ])
 end
 
 # utils
 function update_dist(::Type{Beta}, param::Vector{T}) where T
+    small_threshold!(param, T(1e-8))
     return Beta(param[1], param[2])
 end 
 
diff --git a/src/distributions/common_interface.jl b/src/distributions/common_interface.jl
@@ -23,8 +23,8 @@ end
 function param_tilde_to_param(D::Type{<:Distribution}, param_tilde::Vector{T}) where T 
     return error("param_tilde_to_param not implemented for $D distribution")
 end
-function jacobian_param_tilde(D::Type{<:Distribution}, param_tilde::Vector{T}) where T 
-    return error("jacobian_param_tilde not implemented for $D distribution")
+function jacobian_link(D::Type{<:Distribution}, param_tilde::Vector{T}) where T 
+    return error("jacobian_link not implemented for $D distribution")
 end
 function update_dist(D::Type{<:Distribution}, param::Vector{T}) where T
     return error("update_dist not implemented for $D distribution")
diff --git a/src/distributions/gamma.jl b/src/distributions/gamma.jl
@@ -43,10 +43,10 @@ function param_tilde_to_param(::Type{Gamma}, param_tilde::Vector{T}) where T
         param_tilde_to_param(ExponentialLink, param_tilde[2], zero(T))
     ]
 end
-function jacobian_param_tilde(::Type{Gamma}, param_tilde::Vector{T}) where T 
+function jacobian_link(::Type{Gamma}, 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))
+        jacobian_link(ExponentialLink, param_tilde[1], zero(T));
+        jacobian_link(ExponentialLink, param_tilde[2], zero(T))
     ])
 end
 
diff --git a/src/distributions/log_normal.jl b/src/distributions/log_normal.jl
@@ -40,10 +40,10 @@ function param_tilde_to_param(::Type{LogNormal}, param_tilde::Vector{T}) where T
         param_tilde_to_param(ExponentialLink, param_tilde[2], zero(T))
     ]
 end
-function jacobian_param_tilde(::Type{LogNormal}, param_tilde::Vector{T}) where T 
+function jacobian_link(::Type{LogNormal}, param_tilde::Vector{T}) where T 
     return Diagonal([
-        jacobian_param_tilde(IdentityLink, param_tilde[1]);
-        jacobian_param_tilde(ExponentialLink, param_tilde[2], zero(T))
+        jacobian_link(IdentityLink, param_tilde[1]);
+        jacobian_link(ExponentialLink, param_tilde[2], zero(T))
     ])
 end
 
diff --git a/src/distributions/normal.jl b/src/distributions/normal.jl
@@ -43,10 +43,10 @@ function param_tilde_to_param(::Type{Normal}, param_tilde::Vector{T}) where T
         param_tilde_to_param(ExponentialLink, param_tilde[2], zero(T))
     ]
 end
-function jacobian_param_tilde(::Type{Normal}, param_tilde::Vector{T}) where T 
+function jacobian_link(::Type{Normal}, param_tilde::Vector{T}) where T 
     return Diagonal([
-        jacobian_param_tilde(IdentityLink, param_tilde[1]);
-        jacobian_param_tilde(ExponentialLink, param_tilde[2], zero(T))
+        jacobian_link(IdentityLink, param_tilde[1]);
+        jacobian_link(ExponentialLink, param_tilde[2], zero(T))
     ])
 end
 
diff --git a/src/distributions/poisson.jl b/src/distributions/poisson.jl
@@ -26,7 +26,7 @@ end
 # Links
 param_to_param_tilde(::Type{Poisson}, param::Vector{T}) where T = param_to_param_tilde.(ExponentialLink, param, zero(T))
 param_tilde_to_param(::Type{Poisson}, param_tilde::Vector{T}) where T = param_tilde_to_param.(ExponentialLink, param_tilde, zero(T))
-jacobian_param_tilde(::Type{Poisson}, param_tilde::Vector{T}) where T = Diagonal(jacobian_param_tilde.(ExponentialLink, param_tilde, zero(T)))
+jacobian_link(::Type{Poisson}, param_tilde::Vector{T}) where T = Diagonal(jacobian_link.(ExponentialLink, param_tilde, zero(T)))
 
 # utils
 function update_dist(::Type{Poisson}, param::Vector{T}) where T
diff --git a/src/distributions/weibull.jl b/src/distributions/weibull.jl
@@ -40,10 +40,10 @@ function param_tilde_to_param(::Type{Weibull}, param_tilde::Vector{T}) where T
         param_tilde_to_param(ExponentialLink, param_tilde[2], zero(T))
     ]
 end
-function jacobian_param_tilde(::Type{Weibull}, param_tilde::Vector{T}) where T 
+function jacobian_link(::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))
+        jacobian_link(ExponentialLink, param_tilde[1], zero(T));
+        jacobian_link(ExponentialLink, param_tilde[2], zero(T))
     ])
 end
 
diff --git a/src/link_functions.jl b/src/link_functions.jl
@@ -4,13 +4,13 @@ struct IdentityLink <: Link end
 
 param_to_param_tilde(::Type{IdentityLink}, param::T) where T = param
 param_tilde_to_param(::Type{IdentityLink}, param_tilde::T) where T = param_tilde
-jacobian_param_tilde(::Type{IdentityLink}, param_tilde::T) where T = one(T)
+jacobian_link(::Type{IdentityLink}, param_tilde::T) where T = one(T)
 
 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)
+jacobian_link(::Type{ExponentialLink}, param_tilde::T, lower_bound::T) where T = 1/(param_tilde - lower_bound)
 
 struct LogitLink <: Link end
 
diff --git a/src/score.jl b/src/score.jl
@@ -1,13 +1,12 @@
 function score_tilde(y::T, D::Type{<:Distribution}, param::Vector{T}, param_tilde::Vector{T}, scaling::T) where T
     @assert scaling in (0.0, 1/2, 1.0)
-    # Evaluate Jacobian and score tilde (\tilde \nabla)
-    jac = jacobian_param_tilde(D, param_tilde)
-    score_til = jac*score(y, D, param)
-    # if scaling == 0 do nothing
-    if scaling == 1/2
-        score_til = scaling_invsqrt(jac, D, param)*score_til
+
+    if scaling == 0
+        score_til = scaling_identity(y, D, param)
+    elseif scaling == 1/2
+        score_til = scaling_invsqrt(y, D, param)
     elseif scaling == 1.0
-        score_til = scaling_inv(jac, D, param)*score_til
+        score_til = scaling_inv(y, D, param)
     end
 
     NaN2zero!(score_til)
@@ -17,15 +16,19 @@ function score_tilde(y::T, D::Type{<:Distribution}, param::Vector{T}, param_tild
 end
 
 # Scalings
-function scaling_invsqrt(jac::AbstractMatrix{T}, D::Type{<:Distribution}, param::Vector{T}) where T
-    #TODO improve performace and check if this is right
-    # chol = cholesky(fisher_information(dist, param))
-    # cholmat = Matrix(chol)
-    error("check how to do this properly")
-    return jac'*inv(cholmat)*jac
+function scaling_identity(y::T, D::Type{<:Distribution}, param::Vector{T}) where T
+    jac = jacobian_link(D, param)
+    return inv(jac)*score(y, D, param)
+end
+
+function scaling_invsqrt(y::T, D::Type{<:Distribution}, param::Vector{T}) where T
+    jac = jacobian_link(D, param)
+    inv_fisher = inv(fisher_information(D, param))
+    J_tilde = cholesky(jac*inv_fisher*jac').L
+    return J_tilde * inv(jac) * score(y, D, param)
 end
 
-function scaling_inv(jac::AbstractMatrix{T}, D::Type{<:Distribution}, param::Vector{T}) where T
-    #TODO improve performace
-    return jac'*inv(fisher_information(D, param))*jac
+function scaling_inv(y::T, D::Type{<:Distribution}, param::Vector{T}) where T
+    jac = jacobian_link(D, param)
+    return jac*inv(fisher_information(D, param))*score(y, D, param)
 end
diff --git a/src/utils.jl b/src/utils.jl
@@ -6,34 +6,34 @@ function check_model_estimated(len::Int)
     return false
 end
 
-function NaN2zero!(score_til::Vector{T}) where T
-    for i in eachindex(score_til)
-        if isnan(score_til[i])
-            score_til[i] = zero(T) 
+function NaN2zero!(v::Vector{T}) where T
+    for i in eachindex(v)
+        if isnan(v[i])
+            v[i] = zero(T) 
         end
     end
     return 
 end
 
-function big_threshold!(score_til::Vector{T}, threshold::T) where T
-    for i in eachindex(score_til)
-        if score_til[i] >= threshold
-            score_til[i] = threshold 
+function big_threshold!(v::Vector{T}, threshold::T) where T
+    for i in eachindex(v)
+        if v[i] >= threshold
+            v[i] = threshold 
         end
-        if score_til[i] <= -threshold
-            score_til[i] = -threshold 
+        if v[i] <= -threshold
+            v[i] = -threshold 
         end
     end
     return 
 end
 
-function small_threshold!(score_til::Vector{T}, threshold::T) where T
-    for i in eachindex(score_til)
-        if score_til[i] <= threshold && score_til[i] >= 0
-            score_til[i] = threshold 
+function small_threshold!(v::Vector{T}, threshold::T) where T
+    for i in eachindex(v)
+        if v[i] <= threshold && v[i] >= 0
+            v[i] = threshold 
         end
-        if score_til[i] >= -threshold && score_til[i] <= 0
-            score_til[i] = -threshold 
+        if v[i] >= -threshold && v[i] <= 0
+            v[i] = -threshold 
         end
     end
     return 
diff --git a/test/test_distributions.jl b/test/test_distributions.jl
@@ -1,4 +1,4 @@
-function test_score_mean(D::Type{<:Distribution{Univariate,Continuous}}; n::Int = 10^7, seed::Int = 10,
+function test_score_mean(D::Type{<:Distribution}; n::Int = 10^7, seed::Int = 10,
                             atol::Float64 = 1e-3, rtol::Float64 = 1e-3)
     Random.seed!(seed)
     dist = D()
@@ -10,10 +10,6 @@ function test_score_mean(D::Type{<:Distribution{Univariate,Continuous}}; n::Int
     avg ./= n
     @test avg ≈ zeros(SDM.num_params(D)) atol = atol rtol = rtol
 end
-function test_score_mean(D::Type{<:Distribution}; atol::Float64 = 1e-3, rtol::Float64 = 1e-3)
-    @warn("Score of $D is not currently tested because it has countably infinite support.")
-    return 
-end
 
 function test_loglik(D::Type{<:Distribution}; atol::Float64 = 1e-3, rtol::Float64 = 1e-3,
                      seed::Int = 13, n::Int = 100)
diff --git a/test/test_estimate.jl b/test/test_estimate.jl
@@ -6,7 +6,7 @@ function simulate_GAS_1_1(D::Type{<:Distribution}, scaling::Float64, ω::Vector{
     gas.ω = ω
     gas.A[1] = A
     gas.B[1] = B
-    series, param = simulate(gas, 1000)
+    series, param = simulate(gas, 5000)
 
     return series
 end
@@ -57,7 +57,7 @@ end
         ω = [0.1, 0.1]
         A = [0.5 0; 0 0.5]
         B = [0.5 0; 0 0.5]
-        simulation = simulate_GAS_1_1(Beta, 0.0, ω, A, B, 13)
+        simulation = simulate_GAS_1_1(Beta, 0.0, ω, A, B, 1)
         @testset "Estimation by passing number of seeds" begin
             gas = GAS(1, 1, Beta, 0.0)
             estimate!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
@@ -81,20 +81,24 @@ end
             estimate!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
             test_coefficients_GAS_1_1(gas, ω, A, B)
         end
-        # @testset "Scaling = 0.5" begin
-        #     gas = GAS(1, 1, LogNormal, 0.5)
-        #     estimate!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
-        #     test_coefficients_GAS_1_1(gas)
-        # end
-        # @testset "Scaling = 1.0" begin
-        #     ω = [0.5, 0.5]
-        #     A = [0.1 0; 0 0.1]
-        #     B = [0.1 0; 0 0.1]
-        #     simulation = simulate_GAS_1_1(LogNormal, 1.0, ω, A, B, 123)
-        #     gas = GAS(1, 1, LogNormal, 1.0)
-        #     estimate!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
-        #     test_coefficients_GAS_1_1(gas)
-        # end
+        @testset "Scaling = 0.5" begin
+            ω = [0.1, 0.1]
+            A = [0.5 0; 0 0.5]
+            B = [0.5 0; 0 0.5]
+            simulation = simulate_GAS_1_1(LogNormal, 0.5, ω, A, B, 13)
+            gas = GAS(1, 1, LogNormal, 0.5)
+            estimate!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
+            test_coefficients_GAS_1_1(gas, ω, A, B)
+        end
+        @testset "Scaling = 1.0" begin
+            ω = [0.5, 0.5]
+            A = [0.1 0; 0 0.1]
+            B = [0.1 0; 0 0.1]
+            simulation = simulate_GAS_1_1(LogNormal, 1.0, ω, A, B, 3)
+            gas = GAS(1, 1, LogNormal, 1.0)
+            estimate!(gas, simulation; verbose = 1, opt_method = ScoreDrivenModels.LBFGS(gas, 3))
+            test_coefficients_GAS_1_1(gas, ω, A, B)
+        end
         @testset "GAS([1, 12], [1, 12])" begin
             simulation = simulate_GAS_1_12(LogNormal, 0.0)
             gas = GAS([1, 12], [1, 12], LogNormal, 0.0)
