Metropolis–Hastings algorithm で乱数生成するためのスクリプト

mh.m

function x1 = mh(x0,func,xgen)
  ## input:
  ##   x0   - input distribution
  ##   func - target distribution function
  ##   xgen - random generator for next step 
  ## output:
  ##   x1   - resultant distribution
  xp = xgen(x0);
  y0 = func(x0); y1 = func(xp);
  a  = y1./(y0+(y0==0)) + (y0==0);
  b  = rand(size(y0));

  x1 = x0; x1(a>b) = xp(a>b);
endfunction