0

I am trying to use the Metropolis Hastings algorithm with a random walk sampler to simulate samples from a function $$ in matlab, but something is wrong with my code. The proposal density is the uniform PDF on the ellipse 2s^2 + 3t^2 ≤ 1/4. Can I use the acceptance rejection method to sample from the proposal density?

N=5000;
alpha = @(x1,x2,y1,y2) (min(1,f(y1,y2)/f(x1,x2)));
X = zeros(2,N);
accept = false;
n = 0;
while n < 5000
    accept = false;
    while ~accept
        s = 1-rand*(2);
        t = 1-rand*(2);
        val = 2*s^2 + 3*t^2;
        % check acceptance
        accept = val <= 1/4;
    end
    % and then draw uniformly distributed points checking that u< alpha?
    u = rand();
    c = u < alpha(X(1,i-1),X(2,i-1),X(1,i-1)+s,X(2,i-1)+t);
    X(1,i) = c*s + X(1,i-1);
    X(2,i) = c*t + X(2,i-1);
    n = n+1;
end
figure;
plot(X(1,:), X(2,:), 'r+');
Natalie_94
  • 79
  • 3
  • 12

2 Answers2

0

You may just want to use the native implementation of matlab mhsample.

Regarding your code, there are a few things missing: - function alpha, - loop variable i (it might be just n but it is not suited for indexing since it starts at zero). And you should always allocate memory in matlab if you want to fill it dynamically, i.e. X in your case.

max
  • 3,915
  • 2
  • 9
  • 25
0

To expand on the suggestions by @max, the code appears to work if you change the i indices to n and replace

n = 0;

with

n = 2; X(:,1) = [.1,.1];

It would probably be better to assign X(:,1) to random values within your accept region (using the same code you use later), and/or include a burn-in period.

Depending upon what you are going to do with this, it may also make things cleaner to evaluate the argument to sin in the f function to keep it within 0 to 2 pi (likely by shifting the value by 2 pi if it exceeds those bounds)

Jim Quirk
  • 606
  • 5
  • 18