How to write Matlab code for a multivariate Beta distribution

[glacier302]

Hello,

I am working on a problem in which I first sample two unknown probabilities, p and q, from Beta distributions, and then I want to sample both of them at the same time from a multivariate Beta distribution.

This is the code that I have for sampling p and q individually from Beta distributions (I have already defined a, b, c, d, x, n, and s):


p=betarnd(a+x,b+n-x);
q=betarnd(c+s,d+x+2*(n-x)-s)

Now I'm trying to figure out how to sample p and q jointly from a multivariate Beta distribution.

I know that the Dirichlet distribution is a generalization of the multivariate Beta distribution, and after searching online, I think that I can sample p and q individually from the Dirichlet distribution by doing the following:

alpha_p = [a+x b+n-x];
alpha_q=[c+s d+x+2*(n-x)-s];

d_p=gamrnd(alpha_p, 1);
theta_p=d_p./sum(d_p);

d_q=gamrnd(alpha_q, 1);
theta_q=d_q./sum(d_q);

This results in 2x1 vectors theta_p and theta_q, with the entries in each vector adding up to 1. I'm letting p equal the first entry of theta_p, and q equal the first entry of theta_q.

I'm not sure if this is the right method for the multivariate Beta distribution, but even if it is,
I still don't know how to sample p and q from the Dirichlet distribution at the same time.

Any help would be much appreciated : )

 

[glacier302]

Maybe it would help if I provided more information about the problem:

I'm working on a Gibbs sampler problem. For the regular Gibbs sampler, I used beta distributions to sample the probabilities p and q. I have a vector I that I'm updating at each iteration, which depends on the current values of p and q.

Now I'm trying to write the code for the blocked Gibbs sampler, so I want to sample p and q at the same time instead of individually (right?) I figured that since I used the beta distribution to sample p and q individually, I should use the multivariate beta distribution (dirichlet) to sample them jointly. But I'm not sure how to do this.

This is my code for the regular Gibbs sampler:

a=1;
b=10;
c=10;
d=10;
Data=dlmread('Data.txt');
J=5000;

I=(Data~=1);
n=length(I);

k=1;

for j=1:J

x=sum(I);
s=sum(Data.*(2*ones(n,1)-I)./2);
p=betarnd(a+x,b+n-x);
q=betarnd(c+s,d+x+2*(n-x)-s);

I = -ones(n,1);

for i=1:n

if Data(i)==0
r=rand(1,1);
if r<p*(1-q)/(p*(1-q)+(1-p)*(1-q)^2)
I(i)=1;
else
I(i)=0;
end
end

if Data(i)==1
I(i)=0;
end

if Data(i)==2
r=rand(1,1);
if r<p*q/(p*q+(1-p)*q^2)
I(i)=1;
else
I(i)=0;
end
end

end

if any(I==-1)
warning('sample failed');
end

k=k+1;

end



Any help would be much appreciated!