Thread Subject: Conditional Distributions of the Multivariate Normal

I was wondering if someone could help with the following. I
wish to take random draws from the conditional distributions
of the multivariate normal and have produced the following
code, as an example:

clear;

theta = [1 2 3]';
v = ones(3,3)/2;
v(1,1)=3;
v(2,2)=2;
v(3,3)=1;
mu=zeros(3,1);

for i =1:100
    for h = 1:3
    k=setdiff((1:1:3),h);
     mu(h,:)=theta(h,:)+v(h,k)*inv(v(k,k))*(mu(k,:)-theta(k,:));
    sd(h,:)=sqrt(v(h,h)-v(h,k)*inv(v(k,k))*v(h,k)');
    z(i,h:h)=normrnd(mu(h,:),sd(h,:));
  
    end
end
mean(z)
cov(z)

As can be seen, the mean vector is (1,2,3) and the
covariance vector is [3 .5 .5;.5 2 .5;.5 .5 1]. I am using
the algorithm in Robert (1995) to draw from these
conditional distributions. The mean vector seems to come out
ok, but the covariance structure is off, especially the
off-diagonal elements (diagonals are fine).

Can someone help diagnose the problem or propose a better
solution? My ultimate aim is to draw random variates from
the truncated multivariate normal and this is a start.

Donald Lacombe

>I was wondering if someone could help with the following. I
> wish to take random draws from the conditional distributions
> of the multivariate normal and have produced the following
> code, as an example:

Donald, I don't quite understand what you are trying to do. Are you trying
to generate a multivariate sample by generating univariate values
separately? It look suspicious to me that for some passes through the loop,
mu(h,:) may have either the initial 0 value, the value left over from the
previous random draw (previous i value), or the value set during this random
draw (current i value).

Perhaps you could try one of these, and let me know if it does what you
need. Here n is the number of draws you want, and theta and v are the mean
and covariance that you defined.

z = repmat(theta', n, 1) + randn(n,3)*chol(v);
z = mvnrnd(theta,v,n);

-- Tom

Dear Tom,

Thanks for the reply but I'm actually trying to code a Gibbs
sampler for the multivariate normal distribution which
requires all of the conditional distributions, i.e.
(x1|x2,x3), (x2|x1,x3), and (x3|x1,x2). The formula is
fairly well-known (e.g. Robert (1995) Statistics and
Computing) but it doesn't seem that my code is working well,
hence the query.

Again, thanks for the reply and hopefully I've explained the
situation more clearly.

Donald

"Tom Lane" <tlane@mathworks.com> wrote in message
<fcupqn$7cr$1@fred.mathworks.com>...
> >I was wondering if someone could help with the following. I
> > wish to take random draws from the conditional distributions
> > of the multivariate normal and have produced the following
> > code, as an example:
>
> Donald, I don't quite understand what you are trying to
do. Are you trying
> to generate a multivariate sample by generating univariate
values
> separately? It look suspicious to me that for some passes
through the loop,
> mu(h,:) may have either the initial 0 value, the value
left over from the
> previous random draw (previous i value), or the value set
during this random
> draw (current i value).
>
> Perhaps you could try one of these, and let me know if it
does what you
> need. Here n is the number of draws you want, and theta
and v are the mean
> and covariance that you defined.
>
> z = repmat(theta', n, 1) + randn(n,3)*chol(v);
> z = mvnrnd(theta,v,n);
>
> -- Tom
>
>

Dear Tom,

Thanks for the reply but I'm actually trying to code a Gibbs
sampler for the multivariate normal distribution which
requires all of the conditional distributions, i.e.
(x1|x2,x3), (x2|x1,x3), and (x3|x1,x2). The formula is
fairly well-known (e.g. Robert (1995) Statistics and
Computing) but it doesn't seem that my code is working well,
hence the query.

Again, thanks for the reply and hopefully I've explained the
situation more clearly.

Donald

"Tom Lane" <tlane@mathworks.com> wrote in message
<fcupqn$7cr$1@fred.mathworks.com>...
> >I was wondering if someone could help with the following. I
> > wish to take random draws from the conditional distributions
> > of the multivariate normal and have produced the following
> > code, as an example:
>
> Donald, I don't quite understand what you are trying to
do. Are you trying
> to generate a multivariate sample by generating univariate
values
> separately? It look suspicious to me that for some passes
through the loop,
> mu(h,:) may have either the initial 0 value, the value
left over from the
> previous random draw (previous i value), or the value set
during this random
> draw (current i value).
>
> Perhaps you could try one of these, and let me know if it
does what you
> need. Here n is the number of draws you want, and theta
and v are the mean
> and covariance that you defined.
>
> z = repmat(theta', n, 1) + randn(n,3)*chol(v);
> z = mvnrnd(theta,v,n);
>
> -- Tom
>
>

Dear Tom,

Thanks for the reply but I'm actually trying to code a Gibbs
sampler for the multivariate normal distribution which
requires all of the conditional distributions, i.e.
(x1|x2,x3), (x2|x1,x3), and (x3|x1,x2). The formula is
fairly well-known (e.g. Robert (1995) Statistics and
Computing) but it doesn't seem that my code is working well,
hence the query.

Again, thanks for the reply and hopefully I've explained the
situation more clearly.

Donald

"Tom Lane" <tlane@mathworks.com> wrote in message
<fcupqn$7cr$1@fred.mathworks.com>...
> >I was wondering if someone could help with the following. I
> > wish to take random draws from the conditional distributions
> > of the multivariate normal and have produced the following
> > code, as an example:
>
> Donald, I don't quite understand what you are trying to
do. Are you trying
> to generate a multivariate sample by generating univariate
values
> separately? It look suspicious to me that for some passes
through the loop,
> mu(h,:) may have either the initial 0 value, the value
left over from the
> previous random draw (previous i value), or the value set
during this random
> draw (current i value).
>
> Perhaps you could try one of these, and let me know if it
does what you
> need. Here n is the number of draws you want, and theta
and v are the mean
> and covariance that you defined.
>
> z = repmat(theta', n, 1) + randn(n,3)*chol(v);
> z = mvnrnd(theta,v,n);
>
> -- Tom
>
>

Donald wrote:
> I was wondering if someone could help with the following. I
> wish to take random draws from the conditional distributions
> of the multivariate normal and have produced the following
> code, as an example:
>

Donald,
basically your code was ok. I made some minor corrections.

So first compute the mean and standard deviation of the conditionals.
For demonstration let's compute the conditional x_i|x_j=0,x_k=0.

Now mu gives the expectation and sd.^2 the variance of the conditionals.

Then generate independent sample of these distributions.
Since the samples are independent the off-diagonals of cov(z) should be
zero.

And here is the code

theta = [1 2 3]';
v = ones(3,3)/2;
v(1,1)=3;
v(2,2)=2;
v(3,3)=1;

mu=zeros(3,1); sd=zeros(3,1);
conds=mu;

for h = 1:3
     k=setdiff((1:1:3),h);
      mu(h,:)=theta(h,:)+v(h,k)*inv(v(k,k))*(conds(k,:)-theta(k,:));
      sd(h,:)=sqrt(v(h,h)-v(h,k)*inv(v(k,k))*v(h,k)');
end

ns=1000; z=zeros(ns,3);

for i =1:ns
     z(i,:)=normrnd(mu',sd');
end

[mu sd.^2]

mean(z)
cov(z)


Osmo Kaleva

> Thanks for the reply but I'm actually trying to code a Gibbs
> sampler

I see. But in your code it appears that the various z values are
independent. I think you're using mu in a case where you need to use the
previous values of z. See if this makes sense, where I'm using the new
variable zz to hold the most recent values of all three dimensions of z:

theta = [1 2 3]';
v = ones(3,3)/2;
v(1,1)=3;
v(2,2)=2;
v(3,3)=1;
mu=zeros(3,1);
n = 1000;
zz = zeros(3,1); % new
for i =1:n
    for h = 1:3
    k=setdiff((1:1:3),h);
     mu(h,:)=theta(h,:)+v(h,k)*inv(v(k,k))*(zz(k)-theta(k,:)); v % changed
    sd(h,:)=sqrt(v(h,h)-v(h,k)*inv(v(k,k))*v(h,k)');
    z(i,h:h)=normrnd(mu(h,:),sd(h,:));
    zz(h) = z(i,h); % new

    end
end
mean(z)
cov(z)

(I have not verified your formulas -- just correcting what looks to be a
programming bug.)

-- Tom

Dear Tom and Osmo,

Thank you for taking the time to look at this. Osmo's code
is much faster than what I wrote (thank you!) but the off
diagonals in my "toy" example are not 0 so the variables are
not independent. And I believe that Tom is correct in that
the draws need to be used again in the Gibbs sampling algorithm.

Again, thank you for looking at this. I will try the code
and take a look.

Regards,
Don

"Tom Lane" <tlane@mathworks.com> wrote in message
<fd0ia3$ejr$1@fred.mathworks.com>...
> > Thanks for the reply but I'm actually trying to code a Gibbs
> > sampler
>
> I see. But in your code it appears that the various z
values are
> independent. I think you're using mu in a case where you
need to use the
> previous values of z. See if this makes sense, where I'm
using the new
> variable zz to hold the most recent values of all three
dimensions of z:
>
> theta = [1 2 3]';
> v = ones(3,3)/2;
> v(1,1)=3;
> v(2,2)=2;
> v(3,3)=1;
> mu=zeros(3,1);
> n = 1000;
> zz = zeros(3,1); % new
> for i =1:n
> for h = 1:3
> k=setdiff((1:1:3),h);
>
mu(h,:)=theta(h,:)+v(h,k)*inv(v(k,k))*(zz(k)-theta(k,:)); v
 % changed
> sd(h,:)=sqrt(v(h,h)-v(h,k)*inv(v(k,k))*v(h,k)');
> z(i,h:h)=normrnd(mu(h,:),sd(h,:));
> zz(h) = z(i,h); % new
>
> end
> end
> mean(z)
> cov(z)
>
> (I have not verified your formulas -- just correcting what
looks to be a
> programming bug.)
>
> -- Tom
>
>

Dear Tom,

Just a post script to let you know that your addition to the
code worked like a charm. I have reformulated the formulas
to take advantage of a different parameterization of the
covariance matrix ala Zellner. My error was forgetting that
the Gibbs sampling algorithm uses the draws in subsequent
evaluations in the chain.

Once again, thank you for your help. You have made my week!

Regards,
Don


"Tom Lane" <tlane@mathworks.com> wrote in message
<fd0ia3$ejr$1@fred.mathworks.com>...
> > Thanks for the reply but I'm actually trying to code a Gibbs
> > sampler
>
> I see. But in your code it appears that the various z
values are
> independent. I think you're using mu in a case where you
need to use the
> previous values of z. See if this makes sense, where I'm
using the new
> variable zz to hold the most recent values of all three
dimensions of z:
>
> theta = [1 2 3]';
> v = ones(3,3)/2;
> v(1,1)=3;
> v(2,2)=2;
> v(3,3)=1;
> mu=zeros(3,1);
> n = 1000;
> zz = zeros(3,1); % new
> for i =1:n
> for h = 1:3
> k=setdiff((1:1:3),h);
>
mu(h,:)=theta(h,:)+v(h,k)*inv(v(k,k))*(zz(k)-theta(k,:)); v
 % changed
> sd(h,:)=sqrt(v(h,h)-v(h,k)*inv(v(k,k))*v(h,k)');
> z(i,h:h)=normrnd(mu(h,:),sd(h,:));
> zz(h) = z(i,h); % new
>
> end
> end
> mean(z)
> cov(z)
>
> (I have not verified your formulas -- just correcting what
looks to be a
> programming bug.)
>
> -- Tom
>
>

Dear Donald,

I'm working right now on a similar project that requires sampling from the conditional distributions of a multivariate normal using Gibbs sampling. I was wondering what is the Zellner parametrization of the covariance matrix that you mention. I am basically looking for ways to minimize the computational burden of computing the inverse of the covariance matrix, especially when they are large (n*n, n>100).

Thanks,

Jorge

"Donald " <lacombe.nospam@ohio.edu> wrote in message <fd8s28$di$1@fred.mathworks.com>...
> Dear Tom,
>
> Just a post script to let you know that your addition to the
> code worked like a charm. I have reformulated the formulas
> to take advantage of a different parameterization of the
> covariance matrix ala Zellner. My error was forgetting that
> the Gibbs sampling algorithm uses the draws in subsequent
> evaluations in the chain.
>
> Once again, thank you for your help. You have made my week!
>
> Regards,
> Don
>
>
> "Tom Lane" <tlane@mathworks.com> wrote in message
> <fd0ia3$ejr$1@fred.mathworks.com>...
> > > Thanks for the reply but I'm actually trying to code a Gibbs
> > > sampler
> >
> > I see. But in your code it appears that the various z
> values are
> > independent. I think you're using mu in a case where you
> need to use the
> > previous values of z. See if this makes sense, where I'm
> using the new
> > variable zz to hold the most recent values of all three
> dimensions of z:
> >
> > theta = [1 2 3]';
> > v = ones(3,3)/2;
> > v(1,1)=3;
> > v(2,2)=2;
> > v(3,3)=1;
> > mu=zeros(3,1);
> > n = 1000;
> > zz = zeros(3,1); % new
> > for i =1:n
> > for h = 1:3
> > k=setdiff((1:1:3),h);
> >
> mu(h,:)=theta(h,:)+v(h,k)*inv(v(k,k))*(zz(k)-theta(k,:)); v
> % changed
> > sd(h,:)=sqrt(v(h,h)-v(h,k)*inv(v(k,k))*v(h,k)');
> > z(i,h:h)=normrnd(mu(h,:),sd(h,:));
> > zz(h) = z(i,h); % new
> >
> > end
> > end
> > mean(z)
> > cov(z)
> >
> > (I have not verified your formulas -- just correcting what
> looks to be a
> > programming bug.)
> >
> > -- Tom
> >
> >
>