Can anyone convert this code from matlab to R? Here is what the code is doing:

Simulate X from beta-binomial distribution, by using Gibbs sampler. (step 1): X ? Bin(N,r),

(step 2): r ? Beta(X + ?, N ? X + ?).

Here is the MatLab code:

N = 100;

M = 1000;

a = 5; b=2;

r_alku = rand(1); %Starting value for r

X = zeros(M,1); r = zeros(M,1);

X(1) = binornd(N,r_alku);

r(1) = betarnd(X(1)+a,N-X(1)+b);

for i = 2:M

X(i) = binornd(N,r(i-1));

r(i) = betarnd(X(i)+a,N-X(i)+b);

end

I want to do this in r, where x~bin(n,y) and y~beta(X + ?, N ? X + ?)