# Generating normal random samples given covariance matrix and observations on some components?

Suppose we have n by n covariance matrix of n normal random variables. Is there program where we plug in sample values for k of the n variables and the program, by using the covariance matrix and these k values, generates samples for the remaining n-k variables ?

• we plug in sample values for k of the n variables Please tell more clearly what it is. Do you want to say that you already have some variables with the wanted covariances and now you want to generate more variables following with each other and with the old ones acoording to a full, extended cov/ matrix? – ttnphns Nov 17 '13 at 20:39

Let $Z \sim N_n(μ_Z, Σ_Z)$ and $Z= \left(\begin{array}\\Y\\X\end{array}\right)$, where $X$ has been observed. What you want is to generate from the distribution of $Y|X$.

If $μ_Z$ and $Σ_Z$ are partitioned as follows

$$\mu_Z = \begin{bmatrix} \mu_1 \\ \mu_2 \end{bmatrix}$$

corresponding to the means of the $Y$ and $X$ components respectively, where the lengths of $\mu_1$ and $\mu_2$ are $n-k$ and $k$ respectively, and

$$\Sigma_Z = \begin{bmatrix} \Sigma_{11} &\Sigma_{12} \\ \Sigma_{21} & \Sigma_{22} \end{bmatrix}$$

with $\Sigma_{11}$ being $(n-k) \times (n-k)$, $\Sigma_{12}=\Sigma_{21}^T$ being $(n-k) \times k$ and $\Sigma_{22}$ being $k \times k$, then this is a standard result:

The distribution of $Y$ conditional on $X = \mathbf{x}$ is multivariate normal $(Y|X = \mathbf{x}) \sim N(μ, Σ)$ where

$$\mu = \mu_1 + \Sigma_{12} \Sigma_{22}^{-1} \left( \mathbf{x} - \mu_2 \right)$$

and covariance matrix

$$\Sigma = \Sigma_{11} - \Sigma_{12} \Sigma_{22}^{-1} \Sigma_{21}$$

A derivation is outlined here.

Once you have $\mu$ and $\Sigma$, you can generate as for any multivariate normal, which if you already know them might, for example, be done via Cholesky decomposition.

However, if you don't have that conditional mean and variance already available, the above computation of $\mu$ and $\Sigma$ is not the most efficient or stable way to approach the problem, though; a good way would probably involve working with decompositions more directly, for example avoiding explicit computation of $\Sigma_{22}^{-1}$ at all.

One approach (most likely not the best approach):

If $L$ is the lower triangular Cholesky factor of $\Sigma_{22}$ (that is $LL^T=\Sigma_{22}$ for lower triangular $L$), then you can avoid computation of the inverse by computing $\Sigma_{12} \Sigma_{22}^{-1} \Sigma_{21}$ and $\Sigma_{12} \Sigma_{22}^{-1} \left( \mathbf{x} - \mu_2 \right)$ by solving linear systems involving $L$ or $R=L^T$. Note that

$\Sigma_{12} \Sigma_{22}^{-1} \Sigma_{21} = \Sigma_{12} (LL^T)^{-1} \Sigma_{21}= \Sigma_{12} (L^T)^{-1}L^{-1} \Sigma_{21} = A^TA$, say.

So you need to solve the system $LA = \Sigma_{21}$ to find $A$ (or equivalently solve $BR=\Sigma_{12}$ for $B=A^T$ where $R=L^T$ is the upper, or right Cholesky factor).

With that solved, you could compute $\Sigma = \Sigma_{11}-A^TA$, and if you also solve $L^Tb = \left( \mathbf{x} - \mu_2 \right)$ for $b$, then $\mu=A^Tb$.

[However, there's probably a way to compute the decomposition of $\Sigma$ more directly; I'll try to think about how to that.]

These calculations are readily done in a number of packages.

In R, for example, assuming we already have variables muz, Sigz and x, n and k, being the mean and variance-covariance matrix of the original normal $n$-vector, and the observed $k$-vector as well as the values values of those dimensions, you might do this:

nk <- n-k
ry <- 1:nk
rx <- (nk+1):n

mu1 <- muz[ry]
mu2 <- muz[rx]
S11 <- Sigz[ry,ry]
S12 <- Sigz[ry,rx]   # not used, I just get it for completeness
S21 <- t(S12)
S22 <- Sigz[rx,rx]   # so far just setting up all the partitions

R <- chol(S22)
L <- t(R)

A <- solve(L,S21)
b <- solve(R,x-mu2)
mu <- mu1 + crossprod(A,b)
Sig <- S11 - crossprod(A)


The samples could then be generated by computing the Cholesky decomposition of Sig:

R2 <- chol(Sig)
zeta <- matrix(rnorm(nk*nsim),nr=nk)
dev <- crossprod(R2,zeta)
y <- drop(mu)+dev


provided one has previously defined nsim.

Alternatively, the generation part could be done by installing the mvtnorm package and calling rmvnorm with mu and Sig as the mean and variance arguments.

The calculations could as easily be done in say Matlab, or any number of other packages with similar ease.

Example:

Sigz <- matrix(c(5,3,3,4),nr=2)
muz <- c(10,8)
x=9; n=2; k=1


followed by running the first section of above code to produce mu and Sig:

> mu
[,1]
[1,] 10.75
> Sig
[,1]
[1,] 2.75


and then

nsim <- 3


followed by the second section of code to produce the random values:

> y
[,1]     [,2]     [,3]
[1,] 11.50796 10.76874 13.42943


If $n-k$ was greater than 1 each column would be a vector of values (of length $n-k$), and there would be nsim columns.