Question: MCMC Metropolis-Hastings algorithm - how to improve calculation time?

Dears

Can you give me advice on how to improve the calculation time for Binomial-Beta mcmc Metropolis-Hastings algorithm?

The calculation time for 50 000 samples for my laptop (Dell XPS 9710) takes approx 22 sec.
At the same time, in MathCAD Prime 8.0, the same amount of samples is calculated (using the same procedure) in approx 0.3 s.

======

mcmcBinomialBeta:=proc(a,b,n,k,init,N,burn,thin);

local curr, prop, r, u, fbeta, fbinom, i, p;

curr := <seq(0.0, N)>;
prop := <seq(0.0, N)>;
r := <seq(0.0, N)>;
u := <seq(0.0, N)>;
curr[1] := init;
prop[1] := 0.0;

fbeta := (p, a, b) -> p^(a - 1)*(1 - p)^(-1 + b)/Beta(a, b); 
fbinom := (k, n, p) -> binomial(n, k)*p^k*(1 - p)^(n - k);

for i from 2 to N do 
    prop[i] := Quantile(Normal(curr[i - 1], thin), evalf(rand()/10^12), numeric=true); 
        if 0 < prop[i] and prop[i] < 1 then 
                r[i] := fbinom(k, n, prop[i])*fbeta(prop[i], a, b)/(fbinom(k, n, curr[i - 1])*fbeta(curr[i - 1], a, b)); 
            else 
                r[i] := 0; 
        end if; 
    u[i] := evalf(rand()/10^12); 
        if u[i] < r[i] then 
            curr[i] := prop[i]; 
        else 
            curr[i] := curr[i - 1]; 
        end if; 
end do;
return curr[burn..];

end proc:

======

I will be appreciated any help

wzelik