# score function of bivariate/multivariate normal distribution

I tried to obtain the score vector (1st derivative of density function w.r.t. parameters) of multivariate normal distribution.

Density function of multivariate normal distribution: \begin{align*} f(x) = \frac{1}{\sqrt{(2\pi)^d |\boldsymbol{\Sigma}|}} \exp\left[-\frac{1}{2}(\boldsymbol{x} -\boldsymbol{\mu})^{\rm{T}} \boldsymbol{\Sigma}^{-1} (\boldsymbol{x} -\boldsymbol{\mu})\right]. \end{align*}

By using matrix algebra, I can obtain this quite easily.

\begin{align*} \boldsymbol{u}(y, \theta) = \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \boldsymbol{\theta}} = \begin{bmatrix} \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \theta_1}\\ \vdots\\ \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \theta_{2d + \frac{d(d-1)}{2}}}\\ \end{bmatrix} = \begin{bmatrix} \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \mu_1}\\ \vdots\\ \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \mu_d}\\ \\ \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \sigma_1^2}\\ \vdots\\ \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \sigma_d^2}\\ \\ \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \sigma_{1,2}}\\ \vdots\\ \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \sigma_{d-1,d}}\\ \end{bmatrix} = \begin{bmatrix} \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \boldsymbol{\mu}}\\ \text{diagonal elements of} \left(\frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \boldsymbol{\Sigma}}\right)\\ \text{upper triangular elements of} \left(\frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \boldsymbol{\Sigma}}\right) \end{bmatrix} \end{align*} where \begin{align*} \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \boldsymbol{\mu}} = \boldsymbol{\Sigma}^{-1}(\boldsymbol{y} - \boldsymbol{\mu}) \quad \mbox{and} \quad \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \boldsymbol{\Sigma}} = -\frac{1}{2} \left(\boldsymbol{\Sigma}^{-1} -\boldsymbol{\Sigma}^{-1}(\boldsymbol{y} - \boldsymbol{\mu})(\boldsymbol{y} - \boldsymbol{\mu})^{\rm{T}}\boldsymbol{\Sigma}^{-1}\right). \end{align*}

As a check up, I tried to derive this again. This time with non-matrix calculus in the case of bivariate normal.

The bivariate normal density: \begin{align*} f(x_1, x_2) &= \frac{1}{2\pi\sigma_1 \sigma_2 \sqrt{1-\rho^2}} \exp\left[-\frac{1}{2(1-\rho^2)}\left(\frac{(x_1 -\mu_1)^2}{\sigma_1^2} -\frac{2\rho(x_1 -\mu_1)(x_2 -\mu_2)}{\sigma_1 \sigma_2} +\frac{(x_2 -\mu_2)^2}{\sigma_2^2} \right)\right]\\ &= \frac{1}{2\pi\sigma_1 \sigma_2 \sqrt{1-\frac{\sigma_{1,2}^2}{\sigma_1^2 \sigma_2^2}}} \exp\left[-\frac{1}{2 \left(1-\frac{\sigma_{1,2}^2}{\sigma_1^2 \sigma_2^2}\right)}\left(\frac{(x_1 -\mu_1)^2}{\sigma_1^2} -\frac{2\sigma_{1,2}(x_1 -\mu_1)(x_2 -\mu_2)}{\sigma_1^2 \sigma_2^2} +\frac{(x_2 -\mu_2)^2}{\sigma_2^2} \right)\right]\\ &= \frac{1}{2\pi \sqrt{\sigma_1^2 \sigma_2^2 -\sigma_{1,2}^2}} \exp\left[-\frac{1}{2 (\sigma_1^2 \sigma_2^2 -\sigma_{1,2}^2)}\left(\sigma_2^2 (x_1 -\mu_1)^2 -2\sigma_{1,2}(x_1 -\mu_1)(x_2 -\mu_2) +\sigma_1^2 (x_2 -\mu_2)^2 \right)\right] \end{align*} Log density: \begin{align*} \log f(x_1, x_2) &= -\log 2\pi -\frac{1}{2} \log \left(\sigma_1^2 \sigma_2^2 -\sigma_{1,2}^2\right) \\ &\quad -\frac{1}{2 (\sigma_1^2 \sigma_2^2 -\sigma_{1,2}^2)}\left(\sigma_2^2 (x_1 -\mu_1)^2 -2\sigma_{1,2}(x_1 -\mu_1)(x_2 -\mu_2) +\sigma_1^2 (x_2 -\mu_2)^2 \right) \end{align*}

Elements of the score vector: \begin{align*} \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \theta_1} &= \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \mu_1} = \frac{1}{(1-\rho^2)}\left(\frac{x_1 -\mu_1}{\sigma_1^2} -\frac{\rho(x_2 -\mu_2)}{\sigma_1 \sigma_2}\right)\\ \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \theta_2} &= \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \mu_2} = \frac{1}{(1-\rho^2)}\left(\frac{x_2 -\mu_2}{\sigma_2^2} -\frac{\rho(x_2 -\mu_1)}{\sigma_1 \sigma_2}\right)\\ \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \theta_3} &= \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \sigma_1^2} = \frac{1}{2\left(\sigma_1^2 \sigma_2^2 -\sigma_{1,2}^2 \right)} \left(-\sigma_2^2 +\frac{\sigma_2^2}{\sigma_1^2 \sigma_2^2 -\sigma_{1,2}^2}z -(x_2-\mu_2)^2 \right)\\ \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \theta_4} &= \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \sigma_2^2} = \frac{1}{2\left(\sigma_1^2 \sigma_2^2 -\sigma_{1,2}^2 \right)} \left(-\sigma_1^2 +\frac{\sigma_1^2}{\sigma_1^2 \sigma_2^2 -\sigma_{1,2}^2}z -(x_1-\mu_1)^2 \right)\\ \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \theta_5} &= \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \sigma_{1,2}} = \frac{1}{\sigma_1^2 \sigma_2^2 -\sigma_{1,2}^2} \left(\sigma_{1,2} +(x_1 -\mu_1)(x_2 -\mu_2) -\frac{\sigma_{1,2}}{\sigma_1^2 \sigma_2^2 -\sigma_{1,2}^2}z\right) \end{align*} where \begin{align*} z = \left(\sigma_2^2 (x_1 -\mu_1)^2 -2\sigma_{1,2}(x_1 -\mu_1)(x_2 -\mu_2) +\sigma_1^2 (x_2 -\mu_2)^2 \right). \end{align*}

Then, I implemented both versions in R.

Implementation of the matrix version:

score.mvn.func = function(y, mu, Sigma) {
# This function analytically computes the score vector of the multivariate normal distribution.

# Measure the dimension of the data
d = length(y)

# Safety check
{
# Check whether the input has a proper dimension.
if (
(length(mu) != d) | (length(y) != d) | (dim(Sigma) != d) | (dim(Sigma) != d)
) {
stop("The dimension of the input is not proper.")
}
}

# Frame to write down the result
score.vec = rep(NA, 2*d + choose(d,2))

# Differentiate with respect to 'mu'.
score.vec[1:d] = solve(Sigma)%*%(y - mu)
# Differentiate with respect to 'Sigma'.
logf.by.Sigma = -(1/2)*(solve(Sigma) - solve(Sigma)%*%(y - mu)%*%t(y - mu)%*%solve(Sigma))
# Differentiated by variance
score.vec[(d+1):(2*d)] = diag(logf.by.Sigma)
# Differentiated by covariance
score.vec[(2*d+1):length(score.vec)] = logf.by.Sigma[upper.tri(logf.by.Sigma)]

return(score.vec)
}


Implementation of the non-matrix version:

score.bvn.func = function(y, mu, Sigma) {
# This function analytically computes the score vector of the bivariate normal distribution.

# Measure the dimension of the data
d = 2

# Safety check
{
# Check whether the input has a proper dimension.
if (
(length(mu) != d) | (length(y) != d) | (dim(Sigma) != d) | (dim(Sigma) != d)
) {
stop("The dimension of the input is not proper.")
}
}

# Transform covariate to correlation
rho = Sigma[1,2]/sqrt(Sigma[1,1]*Sigma[2,2])

# Frame to write down the result
score.vec = rep(NA, 2*d + choose(d,2))

# Differentiate with respect to 'mu'.
score.vec = 1/(1 - rho^2)*
(
(y - mu)/(Sigma[1,1]) - (rho*(y - mu))/sqrt(Sigma[1,1]*Sigma[2,2])
)

score.vec = 1/(1 - rho^2)*
(
(y - mu)/(Sigma[2,2]) - (rho*(y - mu))/sqrt(Sigma[1,1]*Sigma[2,2])
)

score.vec = 1/(2*(Sigma[1,1]*Sigma[2,2] - Sigma[1,2]^2))*(
-Sigma[2,2] + Sigma[2,2]/(Sigma[1,1]*Sigma[2,2] - Sigma[1,2]^2)*(
Sigma[2,2]*(y - mu)^2 - 2*Sigma[1,2]*(y - mu)*(y - mu) + Sigma[1,1]*(y - mu)^2
) -
(y - mu)^2
)

score.vec = 1/(2*(Sigma[1,1]*Sigma[2,2] - Sigma[1,2]^2))*(
-Sigma[1,1] + Sigma[1,1]/(Sigma[1,1]*Sigma[2,2] - Sigma[1,2]^2)*(
Sigma[2,2]*(y - mu)^2 - 2*Sigma[1,2]*(y - mu)*(y - mu) + Sigma[1,1]*(y - mu)^2
) -
(y - mu)^2
)

score.vec = (1/(Sigma[1,1]*Sigma[2,2] - Sigma[1,2]^2))*(
Sigma[1,2] + (y - mu)*(y - mu) - (Sigma[1,2]/(Sigma[1,1]*Sigma[2,2] - Sigma[1,2]^2))*(
Sigma[2,2]*(y - mu)^2 - 2*Sigma[1,2]*(y - mu)*(y - mu) + Sigma[1,1]*(y - mu)^2
)
)

return(score.vec)
}


Then, I compared the 2 methods:

> # Simulate data
> set.seed(1)
> X = cbind(rnorm(100,3,1),rnorm(100,2,2))
> mu = colMeans(X)
> Sigma = cov(X)
>
> # Compute score vector with 2 different methods
> score.vec.method.1 = score.mvn.func(y = X[1,], mu = mu, Sigma = Sigma)
> score.vec.method.2 = score.bvn.func(y = X[1,], mu = mu, Sigma = Sigma)
>
> # Display the result
> score.vec.method.1
 -0.91214631 -0.31788468 -0.20375656 -0.08570979  0.14468974
> score.vec.method.2
 -0.91214631 -0.31788468 -0.20375656 -0.08570979  0.28937949


As expected, they give the same result, except for $u_5 = \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \sigma_{1,2}}$

I tried it with different data points and it seems like that $u_5$ from the non-matrix method is exactly 2 times of $u_5$ from the matrix method.

So, \begin{align*} 2\left(\frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \boldsymbol{\Sigma}}\right)_{1,2} = \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \sigma_{1,2}} \end{align*}

instead of \begin{align*} \left(\frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \boldsymbol{\Sigma}}\right)_{1,2} = \frac{\partial \log f(y, \boldsymbol{\theta})}{\partial \sigma_{1,2}} \end{align*}

What did I do wrong? I triple checked both my analytical derivation and the code, but I couldn't find anything that can cause this problem...

• I think that the problem is in the matrix differentiation. Remember that the covariance matrix is symmetric. For details I recomend you to review "A Matrix Handbook for Statisticians" by G. A. F. Seber (2007). Apr 4 '18 at 5:02

$\frac{\partial l}{\partial \Sigma}=-\Sigma^{-1}+\frac{diag(\Sigma^{-1})}{2}+\Sigma^{-1}(x-\mu)(x-\mu)'\Sigma^{-1}-\frac{diag(\Sigma^{-1}(x-\mu)(x-\mu)'\Sigma^{-1})}{2}$
where $l$ is the log-likelihood function.