fn = function(x,a=10,b=-0.2){
a*exp(b*x)
}

fn = function(x,a=0,b=1){
a+b*x
}

samplea = function(x,y,a,b,s,da){
astar = runif(1,a-da,a+da)
logalpha = sum(
dnorm(y,mean=fn(x,astar,b),sd=s,log=T) -
dnorm(y,mean=fn(x,a,b),sd=s,log=T))
logu = log(runif(1,0,1))
acc = (logu < logalpha)
a = acc*astar + (1-acc)*a
list(a=a,acc=acc)
}

sampleb = function(x,y,a,b,s,db){
bstar = runif(1,b-db,b+db)
logalpha = sum(
dnorm(y,mean=fn(x,a,bstar),sd=s,log=T) -
dnorm(y,mean=fn(x,a,b),sd=s,log=T))
logu = log(runif(1,0,1))
acc = (logu < logalpha)
b = acc*bstar + (1-acc)*b
list(b=b,acc=acc)
}

samples = function(x,y,a,b,s,ds){
sstar = abs(runif(1,s-ds,s+ds))
logalpha = sum(
dnorm(y,mean=fn(x,a,b),sd=sstar,log=T) -
dnorm(y,mean=fn(x,a,b),sd=s,log=T)) -
log(sstar/s)
logu = log(runif(1,0,1))
acc = (logu < logalpha)
s = acc*sstar + (1-acc)*s
list(s=s,acc=acc)
}

sampleX = function(x,y,a,b,s,X,dX){
N = length(X)
Xstar = X + runif(N,-dX,dX)
logalpha = (
dnorm(y,mean=fn(Xstar,a,b),sd=s,log=T) +
dnorm(x,mean=Xstar,sd=s,log=T) -
dnorm(y,mean=fn(X,a,b),sd=s,log=T) -
dnorm(x,mean=X,sd=s,log=T) )
logu = log(runif(N,0,1))
acc = (logu < logalpha)
X = acc*Xstar + (1-acc)*X
# X = ifelse(acc,Xstar,X)
list(X=X,acc=acc)
}

samplen = function(n=1000,x,y,
a=0,b=1,s=2,
da=2,db=0.5,ds=2,dX=2)
{
N = length(x)
acca = accb = accs = 0
accX = rep(0,N)
A = B = S = rep(NaN,n)
XX = matrix(NaN,nrow=N,ncol=n)
X = x # Initialize latent variables
for(i in 1:n){
z = samplea(X,y,a,b,s,da)
A[i] = a = z$a
acca = acca + z$acc
z = sampleb(X,y,a,b,s,db)
B[i] = b = z$b
accb = accb + z$acc
z = samples(X,y,a,b,s,ds)
S[i] = s = z$s
accs = accs + z$acc
z = sampleX(x,y,a,b,s,X,dX)
XX[,i] = X = z$X
accX = accX + z$acc
}
invisible(list(
a=A, b=B, s=S, X=XX,
acca=acca/n, accb=accb/n,
accs=accs/n, accX=accX/n))
}