# Aaron Segment 39

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

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

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