function y=monte_carlo(beta,m,k,kappa,a0,b0,varagin)

%MONTE_CARLO    Compute \pi(\beta|\gamma) using Monte Carlo integration.
%               The random draws are from a Gamma-Inverse Gamma
%               distribution.

%Last Modification: 19 March 2008

y=zeros(m,1);


if (nargin==6)

    r=gamrnd(a0,1/b0,1,m);
    r1=find(r);
    lambda=1./r(r1);
    c=gamrnd(kappa*ones(1,length(lambda)),lambda);
    i1=find(c<inf);
    l=length(i1);
    
    y=(c(i1)).^(-k/2).*exp(-0.5./c(i1)*(beta'*beta));

    y=sum(y)/l;

else

    lambda=1./randraw('gamma',[0 1/b0 a0],1,m);
    l=0;

    for j=1:m

        c=randraw('gamma',[0 lambda(j) kappa]);
        
        if (c<inf)
            l=l+1;
            y=y+c^(-k/2)*exp(-0.5/c*(beta'*beta));

        end

    end

    y=y/l;


end
