Skip to content

Commit d2ba869

Browse files
mmkim1210mmkim1210
committed
mo test
1 parent cdd76ab commit d2ba869

3 files changed

Lines changed: 41 additions & 64 deletions

File tree

‎test/eigen_test.jl‎

Lines changed: 15 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -4,10 +4,10 @@ using MultiResponseVarianceComponentModels
44
using BenchmarkTools, LinearAlgebra, Profile, Random, StatsBase, Test
55

66
const MRVCModels = MultiResponseVarianceComponentModels
7-
rng = MersenneTwister(456)
7+
Random.seed!(456)
88

99
n, p, d, m = 855, 3, 4, 2
10-
X = [ones(n) randn(rng, n, p - 1)] # design matrix including intercept
10+
X = [ones(n) randn(n, p - 1)] # design matrix including intercept
1111
# V[1] has entries i * (n - j + 1) for j ≥ i, then scaled to be a correlation matrix
1212
# V[2] is identity
1313
V = Vector{Matrix{Float64}}(undef, m)
@@ -25,7 +25,7 @@ B_true = 2 * rand(p, d) # uniform on [0, 2]
2525
for k in 1:m
2626
Ω_true .+= kron(Σ_true[k], V[k])
2727
end
28-
y = vec(X * B_true) + cholesky(Symmetric(Ω_true)).L * randn(rng, n * d)
28+
y = vec(X * B_true) + cholesky(Symmetric(Ω_true)).L * randn(n * d)
2929
Y = reshape(y, n, d)
3030

3131
@testset "constructor two component" begin
@@ -55,12 +55,12 @@ model = MRVCModel(Y, X, V)
5555
for k in 1:m
5656
println("||Σ_true[$k] - Σ̂[$k]|| = $(norm(Σ_true[k] - model.Σ[k]))")
5757
end
58-
# @test norm(model2.B - model.B) < 2e-4 # 6.986221757875547e-5
59-
# @test norm(model2.Σ[1] - model.Σ[1]) < 2e-4 # 0.0001772074443345503
60-
# @test norm(model2.Σ[2] - model.Σ[2]) < 2e-4 # 1.6274369068536495e-5
61-
# @test abs2(model2.logl[1] - model.logl[1]) < 2e-4 # 4.7903903529347e-8
62-
# @test norm(model2.Bcov - model.Bcov) < 2e-4 # 4.423117207077623e-5
63-
# @test norm(model2.Σcov - model.Σcov) < 2e-4 # 9.284043626500212e-6
58+
@test norm(model2.B - model.B) ≈ 0.00014716893555402874
59+
@test norm(model2.Σ[1] - model.Σ[1]) ≈ 0.00016413437480828665
60+
@test norm(model2.Σ[2] - model.Σ[2]) ≈ 1.959162421647704e-5
61+
@test abs2(model2.logl[1] - model.logl[1]) ≈ 1.1554151166406821e-7
62+
@test norm(model2.Bcov - model.Bcov) ≈ 4.381181547309963e-5
63+
@test norm(model2.Σcov - model.Σcov) ≈ 7.364691385280541e-6
6464
end
6565

6666
model2 = MRTVCModel(Y, X, V, reml = true)
@@ -80,12 +80,12 @@ model = MRVCModel(Y, X, V, reml = true)
8080
for k in 1:m
8181
println("||Σ_true[$k] - Σ̂[$k]|| = $(norm(Σ_true[k] - model.Σ[k]))")
8282
end
83-
# @test norm(model2.B_reml - model.B_reml) < 9e-14 # 8.886934507200367e-14
84-
# @test norm(model2.Σ[1] - model.Σ[1]) < 9e-14 # 1.9800578221407427e-14
85-
# @test norm(model2.Σ[2] - model.Σ[2]) < 9e-14 # 2.9834222365790734e-15
86-
# @test abs2(model2.logl[1] - model.logl[1]) < 9e-14 # 1.987301421658649e-22
87-
# @test norm(model2.Bcov_reml - model.Bcov_reml) < 9e-14 # 5.04991691887946e-15
88-
# @test norm(model2.Σcov - model.Σcov) < 9e-14 # 1.0311308785032728e-15
83+
@test norm(model2.B_reml - model.B_reml) ≈ 1.0372715743775087e-14
84+
@test norm(model2.Σ[1] - model.Σ[1]) ≈ 1.5009812336408224e-14
85+
@test norm(model2.Σ[2] - model.Σ[2]) ≈ 1.914067895379648e-15
86+
@test abs2(model2.logl[1] - model.logl[1]) ≈ 5.169878828456423e-24
87+
@test norm(model2.Bcov_reml - model.Bcov_reml) ≈ 3.3256512190621786e-15
88+
@test norm(model2.Σcov - model.Σcov) ≈ 8.216977255740531e-16
8989
end
9090

9191
end

‎test/fit_test.jl‎

Lines changed: 25 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -2,14 +2,12 @@ module FitTest
22

33
using MultiResponseVarianceComponentModels
44
using BenchmarkTools, LinearAlgebra, Profile, Random, StatsBase, Test
5-
import LinearAlgebra: copytri!
65

76
const MRVCModels = MultiResponseVarianceComponentModels
8-
9-
rng = MersenneTwister(123)
7+
Random.seed!(123)
108

119
n, p, d, m = 855, 3, 4, 3
12-
X = [ones(n) randn(rng, n, p - 1)] # design matrix including intercept
10+
X = [ones(n) randn(n, p - 1)] # design matrix including intercept
1311
# V[1] is an AR1(ρ) matrix, with entries ρ^|i-j|
1412
# V[2] has entries i * (n - j + 1) for j ≥ i, then scaled to be a correlation matrix
1513
# V[3] is identity
@@ -30,21 +28,36 @@ B_true = 2 * rand(p, d) # uniform on [0, 2]
3028
for k in 1:m
3129
Ω_true .+= kron(Σ_true[k], V[k])
3230
end
33-
y = vec(X * B_true) + cholesky(Symmetric(Ω_true)).L * randn(rng, n * d)
31+
y = vec(X * B_true) + cholesky(Symmetric(Ω_true)).L * randn(n * d)
3432
Y = reshape(y, n, d)
3533

34+
@testset "constructor" begin
35+
model = MRVCModel(Y, X[:, 1], V)
36+
model = MRVCModel(Y[:, 1], X, V)
37+
model = MRVCModel(Y[:, 1], X[:, 1], V)
38+
model = MRVCModel(Y, V)
39+
model = MRVCModel(Y[:, 1], V)
40+
model = MRVCModel(Y, X, V[end])
41+
model = MRVCModel(Y, V[end])
42+
model = MRVCModel(Y, X, V, se = false)
43+
model = MRVCModel(Y, X, V)
44+
end
45+
3646
model = MRVCModel(Y, X, V)
3747

3848
@testset "fit! by MLE with MM" begin
39-
@time MRVCModels.fit!(model, algo = :MM, verbose = false, maxiter = 100)
49+
@time MRVCModels.fit!(model, algo = :MM, maxiter = 150)
4050
println("||B_true - B̂|| = $(norm(B_true - model.B))")
4151
for k in 1:m
4252
println("||Σ_true[$k] - Σ̂[$k]|| = $(norm(Σ_true[k] - model.Σ[k]))")
4353
end
54+
println("logl = $(model.logl[1])")
55+
display(model.Bcov)
56+
display(model.Σcov)
4457
end
4558

4659
@testset "fit! by MLE with EM" begin
47-
@time MRVCModels.fit!(model, algo = :EM, verbose = false, maxiter = 100)
60+
@time MRVCModels.fit!(model, algo = :EM, maxiter = 150)
4861
println("||B_true - B̂|| = $(norm(B_true - model.B))")
4962
for k in 1:m
5063
println("||Σ_true[$k] - Σ̂[$k]|| = $(norm(Σ_true[k] - model.Σ[k]))")
@@ -54,54 +67,18 @@ end
5467
model = MRVCModel(Y, X, V; reml = true)
5568

5669
@testset "fit! by REML with MM" begin
57-
@time MRVCModels.fit!(model, algo = :MM, verbose = false, maxiter = 100)
70+
@time MRVCModels.fit!(model, algo = :MM, maxiter = 150)
5871
println("||B_true - B̂|| = $(norm(B_true - model.B_reml))")
5972
for k in 1:m
6073
println("||Σ_true[$k] - Σ̂[$k]|| = $(norm(Σ_true[k] - model.Σ[k]))")
6174
end
75+
println("logl = $(model.logl[1])")
76+
display(model.Bcov_reml)
77+
display(model.Σcov)
6278
end
6379

6480
@testset "fit! by REML with EM" begin
65-
@time MRVCModels.fit!(model, algo = :EM, verbose = false, maxiter = 100)
66-
println("||B_true - B̂|| = $(norm(B_true - model.B_reml))")
67-
for k in 1:m
68-
println("||Σ_true[$k] - Σ̂[$k]|| = $(norm(Σ_true[k] - model.Σ[k]))")
69-
end
70-
end
71-
72-
m = 2
73-
# V[1] has entries i * (n - j + 1) for j ≥ i, then scaled to be a correlation matrix
74-
# V[2] is identity
75-
V = Vector{Matrix{Float64}}(undef, m)
76-
V[1] = [j ≥ i ? i * (n - j + 1) : j * (n - i + 1) for i in 1:n, j in 1:n]
77-
StatsBase.cov2cor!(V[1], [sqrt(V[1][i, i]) for i in 1:n])
78-
V[2] = Matrix(UniformScaling(1.0), n, n)
79-
# true parameter values
80-
Σ_true = [
81-
Matrix(UniformScaling(0.2), d, d),
82-
Matrix(UniformScaling(0.6), d, d)
83-
]
84-
Ω_true = zeros(n * d, n * d)
85-
for k in 1:m
86-
Ω_true .+= kron(Σ_true[k], V[k])
87-
end
88-
y = vec(X * B_true) + cholesky(Symmetric(Ω_true)).L * randn(rng, n * d)
89-
Y = reshape(y, n, d)
90-
91-
model = MRTVCModel(Y, X, V)
92-
93-
@testset "fit! two component by MLE with MM" begin
94-
@time MRVCModels.fit!(model, algo = :MM, verbose = false, maxiter = 100)
95-
println("||B_true - B̂|| = $(norm(B_true - model.B))")
96-
for k in 1:m
97-
println("||Σ_true[$k] - Σ̂[$k]|| = $(norm(Σ_true[k] - model.Σ[k]))")
98-
end
99-
end
100-
101-
model = MRTVCModel(Y, X, V; reml = true)
102-
103-
@testset "fit! two component by REML with MM" begin
104-
@time MRVCModels.fit!(model, algo = :MM, verbose = false, maxiter = 100)
81+
@time MRVCModels.fit!(model, algo = :EM, maxiter = 150)
10582
println("||B_true - B̂|| = $(norm(B_true - model.B_reml))")
10683
for k in 1:m
10784
println("||Σ_true[$k] - Σ̂[$k]|| = $(norm(Σ_true[k] - model.Σ[k]))")

‎test/runtests.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,3 @@
1+
include("fit_test.jl")
12
include("eigen_test.jl")
2-
# include("fit_test.jl")
33
# include("mvcalculus_test.jl")

0 commit comments

Comments
 (0)