Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

How to write Matlab code for a multivariate Beta distribution

  1. Dec 20, 2011 #1
    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 : )
     
  2. jcsd
  3. Dec 21, 2011 #2
    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!
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: How to write Matlab code for a multivariate Beta distribution
Loading...