Conditional Normal Random Distribution

I want to generate a random vector (a) from a normal distribution N(mu,sigma) (mu,sigma:known) with a condition that the first n values of vector 'a' are known and fixed (basically its fulfilling boundary conditions).

Is there any way I can use Multivariate normal random numbers function: R = mvnrnd(mu,Sigma) or any other method?

Answer 1

[Paul]

Yes. If you know mu and Sigma of the vector x and the first n values of x are given, then the density of x(n+1:end) is also normal and can be derived from mu, Sigma, and x(1:n). See this link for the math to get the mean and covariance of x(n+1:end) condtioned on x(1:n), then you can use mvnrnd to generate random numbers of x(n+1:end)

[Walter Roberson]

How does knowing x(1:n) help?

It looks to me as if the solution would simply be to take append values drawn from N(mu,sigma) to the end of the known values ?

The required output is a vector and that tells us that only one variable is involved. mvnrnd() does hint that multiple variables are involved, but mvnrnd() generates rows with number of columns according to mu or sigma, with the number of columns implicitly indicating the number of variables involved.

Or is the implication here that there are M variables, where M is the desired length of the output vector, and only one multivariate "drawing" is being done, with the first n outputs known?

[Paul]

My interpretation of the question is as follows:

Let X be a random vector (numel(X) = n + m) where the elements of X have a joint normal density function with mean mu and covariance Sigma. p random vectors of X can be generated with

X = mvnrnd(mu,Sigma,p);

The question implies that the need is to generate random values of X(n+1:n+m) (let's call it Xsub), becaue the first n values of X are given. If we didn't have any information about the known values of X, then the elements of Xsub also have a joint normal density with musub = mu(n+1:n+m) and Sigmasub = Sigma(n+1:n+m,n+1:n+m).

However, the question states that first n values of X are known. With that additional information, the joint denity of the elements of Xsub is still normal, but the mean and covariance of Xsub conditioned on the known values of x(1:n), called mubar and Sigmabar, can be quite different than musub and Sigmasub. The link shows how to compute mubar and Sigmabar based on mu, Sigma, and the given elements of X, which the link calls 'a'.

Note that the link partitions X the opposite of the problem here, i.e., the link shows what to do when the lower part of X is known, but that's easily handled by, for example, reindexing X and changing mu and Sigma appropriately.

Once mubar and Sigmabar are computed, then random vectors of Xsub can be realized by

Xsub = mvnrnd(mubar,Sigmabar,p);

and p realizations of the full vector X would then then be

X = [repmat(a(:).' , p , 1) Xsub]; 

[Vinit Nagda]

@Paul I am getting some error with the size of mubar. Could you please help?

The size of mu is [px1] and for sigma is [pxp]. (p=n+m)

When I calculate mubar (using the formula from the link you shared), I get the size of mubar as [mxn](sigmabar is [mxm]) and when I use this mubar in mvnrnd function, it gives me size error.

[Paul]

mubar should be mx1. If you post your code, maybe someone can help you get it sorted out.

[Vinit Nagda]

@Paul The error is fixed now. The variable with known values was not defined properly. Thank you for your help.

[Paul]

Glad to hear it. Once you get mubar and Sigmabar, you can check your code (for a reasonably small dimension of X) by generating many realizations of X, for example,

X = mvnrnd(mu,Sigma,1e5);

and then only keeping the rows of X where the columns are within a tolerance of the "known values"

keep = abs(X(:,1) - givenvalue1) < tol & abs(X(:,2) - givenvalue2) < tol ; % example with two given values
Xsub = X(keep,:)

and then compare the estimated mean and covariance of Xsub to what they should be

mubarestimate = mean(Xsub);
Sigmabarestimate = cov(Xsub);

[Vinit Nagda]

@Paul Hi, I tried to check and compare your example with my code but the result is not at all similar.

Could you please have a look at my code and see if you find any mistake? Sorry for the trouble.

Thank you very much.

clc
clear all
v=[14.26 7.13 1.91 0.27 0.13 0.34 0.17 0.16 0.04 0.06];
sigma=diag(v);
mu=zeros(size(sigma,1),1);
n=size(sigma,1);
ind=8;%Length of x1
values=[1.6 2.1]'; %known values, x2
s12s22=sigma(1:ind,(ind+1):n)*sigma((ind+1):n,(ind+1):n)^(-1);
mubar=mu((1:ind),1)+s12s22*(values-mu(((ind+1):n),1));
sigmabar=sigma(1:ind,1:ind)-s12s22*sigma((ind+1):n,1:ind);

[Paul]

Hi @Vinit Nagda

Offhand I don't see anything wrong with your code, but there may be a few things you can do to improve it for clarity.

Basically, it would be clearer to define variables in accordance with the equations you're going to implement, instead of using the explicit indexing in those equations. That is, define the variables mu1, mu2, Sigma11, Sigma12, Sigma21 and Sigma22 based on mu, sigma, and ind, and then use the mu* and Sigma* variables as needed. Sacrifice a very small amount of efficiency (at leadt for this example) for a lot of clarity.

Perhaps rename the variable s12s22 to s12s22inv.

Instead of hardcoding ind to a value, compute it based on numel(v) and numel(values).

Finally the use of ^(-1) is disouraged for a matrix inverse. In this instance, use the / operator, i.e.,

s12s22inv = Sigma12/Sigma22;

Having said all of that, this example data might not be useful for illustrating the concept of conditional probability density. The reason is that sigma (which I would call Sigma) is diagonal, which means that all 10 of the random variables in X are jointly independent. Joint independence means that knowledge of the values of any subset of the random variables does not give any additional information about the others that will change their probability densities. Hence, the probability density of the unknown variables is the same as what you started with, which is exactly what your code shows, i.e. mubar = mu(1:8) and Sigmabar = Sigma11).

What different result were you expecting?

[Vinit Nagda]

@Paul Thank you very much for your response.

Actually, I didn't know what result to expect and I just compared your example with my code and thought something was wrong. But, I now understand the joint independence concept and thanks a lot for sharing that.

Also, thanks for your suggestions, I will try to make my code look more code-like with simplifications and clarity.

Just, one last final thing (I will then close this thread and not disturb you more). As in the link you shared, lower part of x(x2) is known and in my case x_1 is known so I have to do reindexing.

The order of elements in random normal vector and mu vector can be easily reversed, I just want to confirm the reindexing of Sigma matrix. Flipping the elements of matrix across diagonal and reversing the diagonal elements should work right?

Thank you. I very much appreciate your time and help.

[Paul]

@Vinit Nagda

I didn't provide an example, so I'm still not sure what you've compared to your code.

Instead of reindexing your problem to fit the formula, it's probably easier to modfiy the formula to fit your problem. Just reverse the 1's and 2's.

mubar    = mu2 + Sigma21/Sigma11*(a - mu1);
Sigmabar = Sigma22 - Sigma21/Sigma11*Sigma12;

Feel free to come back if you have any more questions, particularly if you want to post your code with an example. However, if you do, don't use a 10-element vector in the example. Maybe use a 3-element vector X and show the mubar and Sigmabar of X(3) given known values of X(1) and X(2).