Due May 1 @ 11:59PM
versioninfo()
Consider a linear mixed effects model $$ \mathbf{Y}_i = \mathbf{X}_i \boldsymbol{\beta} + \mathbf{Z}_i \boldsymbol{\gamma} + \boldsymbol{\epsilon}_i, \quad i=1,\ldots,n, $$ where
Write down the log-likelihood of the $i$-th datum $(\mathbf{y}_i, \mathbf{X}_i, \mathbf{Z}_i)$ given parameters $(\boldsymbol{\beta}, \boldsymbol{\Sigma}, \sigma^2)$.
Use the following template to define a type LmmObs
that holds an LMM datum $(\mathbf{y}_i, \mathbf{X}_i, \mathbf{Z}_i)$.
# define a type that holds LMM datum
struct LmmObs{T <: AbstractFloat}
# data
y :: Vector{T}
X :: Matrix{T}
Z :: Matrix{T}
# working arrays
# whatever intermediate arrays you may want to pre-allocate
res :: Vector{T}
storage_q :: Vector{T}
ztz :: Matrix{T}
storage_qq :: Matrix{T}
end
# constructor
function LmmObs(
y::Vector{T},
X::Matrix{T},
Z::Matrix{T}) where T <: AbstractFloat
res = similar(y)
storage_q = Vector{T}(undef, size(Z, 2))
ztz = transpose(Z) * Z
storage_qq = similar(ztz)
LmmObs(y, X, Z, res, storage_q, ztz, storage_qq)
end
Write a function, with interface
logl!(obs, β, L, σ²)
that evaluates the log-likelihood of the $i$-th datum. Here L
is the lower triangular Cholesky factor from the Cholesky decomposition Σ=LL'
. Make your code efficient in the $n_i \gg q$ case. Think the intensive longitudinal measurement setting.
function logl!(
obs :: LmmObs{T},
β :: Vector{T},
L :: Matrix{T},
σ² :: T) where T <: AbstractFloat
n, p, q = size(obs.X, 1), size(obs.X, 2), size(obs.Z, 2)
# TODO: compute and return the log-likelihood
sleep(1e-3) # wait 1 ms as if your code takes 1ms
return 0.0
end
Hint: This function shouldn't be very long. Mine, obeying 80-character rule, is 25 lines. If you find yourself writing very long code, you're on the wrong track. Think about algorithm first then use BLAS functions to reduce memory allocations.
Compare your result (both accuracy and timing) to the Distributions.jl package using following data.
using BenchmarkTools, Distributions, LinearAlgebra, Random
Random.seed!(257)
# dimension
n, p, q = 2000, 5, 3
# predictors
X = [ones(n) randn(n, p - 1)]
Z = [ones(n) randn(n, q - 1)]
# parameter values
β = [2.0; -1.0; rand(p - 2)]
σ² = 1.5
Σ = fill(0.1, q, q) + 0.9I
# generate y
y = X * β + Z * rand(MvNormal(Σ)) + sqrt(σ²) * randn(n)
# form an LmmObs object
obs = LmmObs(y, X, Z)
This is the standard way to evaluate log-density of a multivariate normal, using the Distributions.jl package. Let's evaluate the log-likelihood of this datum.
μ = X * β
Ω = Z * Σ * transpose(Z) + σ² * I
mvn = MvNormal(μ, Symmetric(Ω)) # MVN(μ, Σ)
logpdf(mvn, y)
Check that your answer matches that from Distributions.jl
L = Matrix(cholesky(Σ).L)
logl!(obs, β, L, σ²)
You will lose all 15 + 30 + 30 = 75 points if the following statement throws AssertError
.
@assert logl!(obs, β, Matrix(cholesky(Σ).L), σ²) ≈ logpdf(mvn, y)
Benchmarking your code and compare to the Distributions.jl function logpdf
.
# benchmark the `logpdf` function in Distribution.jl
bm1 = @benchmark logpdf($mvn, $y)
# benchmark your implementation
L = Matrix(cholesky(Σ).L)
bm2 = @benchmark logl!($obs, $β, $L, $σ²)
The points you will get is $$ \frac{x}{1000} \times 30, $$ where $x$ is the speedup of your program against the standard method.
# this is the points you'll get
clamp(median(bm1).time / median(bm2).time / 1000 * 30, 0, 30)
Hint: Apparently I am using 1000 as denominator because I expect your code to be at least $1000 \times$ faster than the standard method.
You want to avoid memory allocation in the "hot" function logl!
. You will lose 1 point for each 1 KiB = 1024 bytes
memory allocation. In other words, the points you get for this question is
clamp(30 - median(bm2).memory / 1024, 0, 30)
Hint: I am able to reduce the memory allocation to 0 bytes.
Coding style, Git workflow, etc. For reproducibity, make sure we (TA and myself) can run your Jupyter Notebook. That is how we grade Q4 and Q5. If we cannot run it, you will get zero points.