Aaron Segment 39

From Computational Statistics Course Wiki
Jump to navigation Jump to search
To Calculate
1. Suppose the domain of a model are the five integers x = \{1,2,3,4,5\}, and that your proposal distribution is: "When x_1 = 2,3,4, choose with equal probability x_2 = x_1 \pm 1. For x_1=1 always choose x_2 =2. For x_1=5 always choose x_2 =4. What is the ratio of q's that goes into the acceptance probability \alpha(x_1,x_2) for all the possible values of x_1 and x_2?
A: Below is my code and graphs
clear all
clc

%% Running a chain with these rules
x=[];
x(1)=randi(5);
A=[1 2 3 4 5];
inter=[2 3 4];
for i=1:1000
    if ismember(x(i),inter)==1
        p=rand(1);
        if p<=.5;
            x(i+1)=x(i)+1;
        else
            x(i+1)=x(i)-1;
        end
    elseif x(i)==5
            x(i+1)=4;
    else
        x(i+1)=2;
    end
    
end
figure(1)
hist(x,100)
p=[sum(x==1) sum(x==2) sum(x==3) sum(x==4) sum(x==5)];
psep=p./sum(p)

%% Pulling random samples and determining acceptance probs

y=[];
y(1)=randi(5);
den=0;
for j=1:1000
    s=randi(5);
    if ismember(y(j),inter)==1
        test=[y(j)+1 y(j)-1];
        if ismember(s,test)==1
            y(j+1)=s;
        else
            y(j+1)=y(j);
        end
    elseif y(j)==5;
        if s==4
            y(j+1)=s;
        else
            y(j+1)=y(j);
        end
    else
        if s==2;
            y(j+1)=s;
        else
            y(j+1)=y(j);
        end
    end
    if y(j+1)==y(j)
        den=den+1;
    end
    
end
figure(2)
hist(y,100)
allacceptperc=den/j

% Testing answer to first question with x1=3, x2=4
three=0;
l=0;
five=0;
for k=2:length(y)
    if y(k)==4
        l=l+1;
        if y(k-1)==3
            three=three+1;
        elseif y(k-1)==5;
            five=five+1;
        end
    end
end
threep=three/k
fivep=five/k
%% Comments for calculate problems
% If x2=3 then the ratio of qs is 1 if x1=2,4 0 otherwise
%
% If x2=2,4 then the ratio of qs with x1=3 is 1. 0 otherwise
%
% IFf x2=2,4 with x=1,5 then the ratio is 1, 0 otherwise
%
% IF x2=1,5 then if x1=2,4 respectively, then the q ratio is 1, 0 otherwise
%



psep =

    0.1249    0.2438    0.2438    0.2567    0.1309


allacceptperc =

    0.6740


threep =

    0.0340


fivep =

    0.0350

AaronSeg39Hist1.jpg AaronSeg39hist2.jpg

2. Suppose the domain of a model is -inf < x < inf and your proposal distribution is (perversely), q(x_2|x_1) = \begin{cases}\tfrac{7}{2}\exp[-7(x_2-x_1)],\quad & x_2 \ge x_1 \\ \tfrac{5}{2}\exp[-5(x_1-x_2)],\quad & x_2 < x_1 \end{cases}
clear all
clc

%% I'm attempting to pull samples from N(5,9) on the interval 0 to 5
% This is something I can test very easily using invCDF or a Matlab rand
% q will be normal dist centered at xi with std 1.
% this is symetric so alpha is just pi(x)/pi(x1+)
alpha=[];
x(1)=3;
for i=1:10000
    test=normrnd(x(i),1);
    stat=normpdf(test,5,9)/normpdf(x(i),5,9);
    alpha(i)=min(1,stat);
    alphcomp=rand(1);
    if alpha(i)>alphcomp
        x(i+1)=test;
    else
        x(i+1)=x(i);
    end
end
figure(1)
hist(x,100)
%xlim([0,8])
mean(x)
std(x)
y=x(x>=0);
comp=y(y<=7);
figure(2)
hist(comp,100)
normsv=normrnd(5,9,[10000 1]);
figure(3)
hist(normsv,1000)
xlim([0 7])
Here is an additional attempt at Gibbs sampling with something I could compare it to. The results were almost identical to pulling samples using the inverse cdf method (thats exciting)!
%% Attempting Gibbs sampling with f(x,y)=nchoosek(16 x)*y^(x-3)*(1-y)^(19-x)
% therefore f(x|y)=binomial(16,y)
% f(y|x)=Beta(x+2,20-x)
clear all 
clc

y(1)=1;
for i=1:1000
    x(i)=binornd(16,y(i));
    y(i+1)=betarnd(x(i)+2,20-x(i));
end
figure(1)
hist(x,100)

AaronGibbsSeg39.jpg

Sketch this distribution as a function of x_2-x_1. Then, write down an expression for the ratio of q's that goes into the acceptance probability \alpha(x_1,x_2).

To Think About
1. Suppose an urn contains 7 large orange balls, 3 medium purple balls, and 5 small green balls. When balls are drawn randomly, the larger ones are more likely to be drawn, in the proportions large:medium:small = 6:4:3. You want to draw exactly 6 balls, one at a time without replacement. How would you use Gibbs sampling to learn: (a) How often do you get 4 orange plus 2 of the same (non-orange) color? (b) What is the expectation (mean) of the product of the number of purple and number of green balls drawn?
2. How would you do the same problem computationally but without Gibbs sampling?
3. How would you do the same problem non-stochastically (e.g., obtain answers to 12 significant figures)? (Hint: This is known as the Wallenius non-central hypergeometric distribution.)
[Answers: 0.155342 and 1.34699]