Derivatives for original model

.gitignore

@@ -3,3 +3,4 @@
 *.jl.cov
 *.jl.mem
 /Manifest.toml
+test_code.jl

Manifest.toml

@@ -2,7 +2,7 @@
 
 julia_version = "1.10.2"
 manifest_format = "2.0"
-project_hash = "74736c4a072a7d041f11add2f815261972da8a49"
+project_hash = "fc4f4328fbdcbb9066a8f2f59f300b0c66de4908"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]

Project.toml

@@ -7,7 +7,9 @@ version = "1.0.0-DEV"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]

src/UnobservedCountEstimation.jl

@@ -3,8 +3,9 @@ module UnobservedCountEstimation
 # Write your package code here.
 
 # Imports
-using Optim, Statistics, GLM, Distributions, DataFrames
+using Optim, Statistics, GLM, Distributions, DataFrames, SpecialFunctions, LinearAlgebra
 
+include("zhang_likelihood.jl")
 include("original_model.jl")
 include("interface.jl")
 

src/original_model.jl

@@ -8,12 +8,23 @@ function zhang_model(m, N, n; start = "glm")
     if start == "glm"
         start = coef(glm(@formula(y ~ x1 + x2 + 0), df, Poisson(), LogLink()))
     else
-        start = coef(m(@formula(log(y) ~ x1 + x2 + 0), df))
+        start = coef(lm(@formula(log(y) ~ x1 + x2 + 0), df))
     end # end if
 
     #= mm  = glm(@formula(y ~ x1 + x2 + 0), df, Poisson(), LogLink())
     ols =  lm(@formula(log(y) ~ x1 + x2 + 0), df) =#
+    X = ones(length(n))
+    X = X[:, :]
+
+    Z = ones(length(n))
+    Z = Z[:, :]
+
+    log_l = log_lik_original_model(start[1], start[2], 1, m, N, n, X, Z)
+
+    grad_l = grad_log_lik_original_model(start[1], start[2], 1, m, N, n, X, Z)
+
+    hes_l = hess_log_lik_original_model(start[1], start[2], 1, m, N, n, X, Z)
 
     #[coef(ols), coef(mm)]
-    start
+    [start, log_l, grad_l, hes_l]
 end # end function

src/zhang_likelihood.jl

@@ -0,0 +1,40 @@
+# likelihood and derivatives go here
+
+function log_lik_original_model(α, β, ϕ, m, N, n, X, Z)
+    μ = (N .^ (X * α)) .* ((n ./ N) .^ (Z * β))
+    sum(m .* log.(μ) .- (m .+ ϕ) .* log.(μ .+ ϕ) .+ ϕ .* log.(ϕ) .- loggamma.(ϕ) .- loggamma.(m .+ 1) .+ loggamma.(m .+ ϕ))
+end # end function
+
+function grad_log_lik_original_model(α, β, ϕ, m, N, n, X, Z)
+    # SpecialFunctions.jl has digamma and trigamma function
+    μ = (N .^ (X * α)) .* ((n ./ N) .^ (Z * β))
+    dμ = m ./ μ .- (ϕ .+ m) ./ (μ .+ ϕ)
+
+    dα = dμ .* μ .* log.(N)
+    dβ = dμ .* μ .* log.(n ./ N)
+    dϕ = sum(SpecialFunctions.digamma.(ϕ .+ m) .- log.(μ .+ ϕ) .+ log.(ϕ) .- SpecialFunctions.digamma.(ϕ) .+ 1 .- (ϕ .+ m) ./ (ϕ .+ μ))
+
+    vcat(X' * dα, Z' * dβ, dϕ)
+end # end function
+
+function hess_log_lik_original_model(α, β, ϕ, m, N, n, X, Z)
+    μ = (N .^ (X * α)) .* ((n ./ N) .^ (Z * β))
+    dμ = m ./ μ .- (ϕ .+ m) ./ (μ .+ ϕ)
+    dμ_2 = (ϕ .+ m) ./ (μ .+ ϕ) .^ 2 .- m ./ μ .^ 2
+    dμdϕ = - (μ .- m) ./ (μ .+ ϕ) .^ 2
+
+    #dϕ = sum(SpecialFunctions.digamma.(ϕ .+ m) .- log.(μ .+ ϕ) .+ log.(ϕ) .- SpecialFunctions.digamma.(ϕ) .+ 1 .- (ϕ .+ m) ./ (ϕ .+ μ))
+    dϕ_2 = sum(SpecialFunctions.trigamma.(ϕ .+ m) .- SpecialFunctions.trigamma.(ϕ) .- 2 ./ (ϕ .+ μ) .+ (ϕ .+ m) ./ (ϕ .+ μ) .^ 2 .+ ϕ ^ -1)
+
+    dα_2 = dμ_2 .* (μ .* log.(N)) .^ 2 .+ dμ .* μ .* (log.(N) .^ 2)
+    dβ_2 = dμ_2 .* (μ .* log.(n ./ N)) .^ 2 .+ dμ .* μ .* (log.(n ./ N) .^ 2)
+    dαdβ = dμ_2 .* μ .* log.(N) .* μ .* log.(n ./ N) .+ dμ .* μ .* log.(N) .* log.(n ./ N)
+    dαdϕ = dμdϕ .* μ .* log.(N)
+    dβdϕ = dμdϕ .* μ .* log.(n ./ N)
+
+    vcat(
+        hcat(X' * (dα_2 .* X), X' * (dαdβ .* Z), X' * dαdϕ),
+        hcat((X' * (dαdβ .* Z))', Z' * (dβ_2 .* Z), Z' * dβdϕ),
+        hcat((X' * dαdϕ)', (Z' * dβdϕ)', dϕ_2)
+    )
+end # end function

test/runtests.jl

@@ -5,5 +5,5 @@ using Test
 
 @testset "UnobservedCountEstimation.jl" begin
     # TODO  
-    @test UnobservedCountEstimation.placeholder(65) == 65
+    @test placeholder(65) == 65
 end