# Aaron Segment 39

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
% 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


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)


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).