R/isingibbs.R

In bayess: Bayesian Essentials with R

isingibbs=function(niter,n,m=n,beta){
  # initialization
  x=sample(c(0,1),n*m,prob=c(0.5,0.5),replace=TRUE)
  x=matrix(x,n,m)

  # Gibbs runs
  for (i in 1:niter){
    sampl1=sample(1:n)
    sampl2=sample(1:m)
    for (k in 1:n){
    for (l in 1:m){
     n0=xneig4(x,sampl1[k],sampl2[l],0)
     n1=xneig4(x,sampl1[k],sampl2[l],1)
     x[sampl1[k],sampl2[l]]=sample(c(0,1),1,prob=exp(beta*c(n0,n1)))
     }}}
  x
}