abcnorm <- function(x,m=0,c=1,sigma=1,iters=1000,errprop=0.01){
#
# USe of the ABC method accepting a fixed proportion of the error values and with a sufficient statistic.
#  
  n <- length(x)
  xbar <- mean(x)
  niters <- ceiling(iters/errprop)
  mu <- rnorm(niters,m,sigma/sqrt(c))
  xbarsim <- rnorm(niters,mu,sigma/sqrt(n))
  errsim <- abs(xbarsim-xbar)
  e <- sort(errsim,index.return=TRUE)
  err <- e$x
  e <- e$ix
  mu <- mu[e]
  err <- err[iters]
  mu <- mu[1:iters]
  return(list('mu'=mu,'err'=err))
}

abcnormmedian <- function(x,m=0,c=1,sigma=1,iters=1000,errprop=0.01){
#
# USe of the ABC method accepting a fixed proportion of the error values and with an insufficient statistic.
#  
  n <- length(x)
  xmed <- median(x)
  niters <- ceiling(iters/errprop)
  mu <- rnorm(niters,m,sigma/sqrt(c))
  xmedsim <- rep(0,niters)
  for (i in 1:niters){
    xmedsim[i] <- median(rnorm(n,mu[i],sigma))
  }
  errsim <- abs(xmedsim-xmed)
  e <- sort(errsim,index.return=TRUE)
  err <- e$x
  e <- e$ix
  mu <- mu[e]
  err <- err[iters]
  mu <- mu[1:iters]
  return(list('mu'=mu,'err'=err))
}

abcnorms2 <- function(x,m=0,c=1,sigma=1,iters=1000,errprop=0.01){
  #
  # USe of the ABC method accepting a fixed proportion of the error values and with an insufficient statistic.
  #  
  n <- length(x)
  xs2 <- var(x)
  niters <- ceiling(iters/errprop)
  mu <- rnorm(niters,m,sigma/sqrt(c))
  xs2sim <- rep(0,niters)
  for (i in 1:niters){
    xs2sim[i] <- var(rnorm(n,mu[i],sigma))
  }
  errsim <- abs(xs2sim-xs2)
  e <- sort(errsim,index.return=TRUE)
  err <- e$x
  e <- e$ix
  mu <- mu[e]
  err <- err[iters]
  mu <- mu[1:iters]
  return(list('mu'=mu,'err'=err))
}

#
# Sample code for the original method ...
#
truemu <- 1
sigma <- 1
n <- 50
x <- rnorm(n,truemu,sigma)
#
m <- 0
c <- 1
musd <- sigma/sqrt(c+n)
mumean <- (c*m+sum(x))/(c+n)
mugrid <-mumean-3*musd+6*musd*c(0:1000)/1000
fmu <- dnorm(mugrid,mumean,musd)
#
iters <- 1000
errprop <- 0.005
e <- abcnorm(x,m,c,sigma,iters,errprop)
#
mu <- e$mu
err <- e$err
d <- density(mu)
hist(mu,probability=TRUE,xlim=c(min(mugrid[1],d$x[1]),max(mugrid[1001],d$x)),ylim=c(0,max(fmu,d$y)),nclass=20,xlab='mu',ylab='f',main='')
lines(d,type='l',lwd=2,col='green')
lines(mugrid,fmu,type='l',lwd=2,col='red')
#
# ... and for the second method
#
e <- abcnormmedian(x,m,c,sigma,iters,errprop)
mu <- e$mu
err <- e$err
d <- density(mu)
hist(mu,probability=TRUE,xlim=c(min(mugrid[1],d$x[1]),max(mugrid[1001],d$x)),ylim=c(0,max(fmu,d$y)),nclass=20,xlab='mu',ylab='f',main='')
lines(d,type='l',lwd=2,col='green')
lines(mugrid,fmu,type='l',lwd=2,col='red')
#
# ... and for the third
#
e <- abcnorms2(x,m,c,sigma,iters,errprop)
mu <- e$mu
err <- e$err
d <- density(mu)
mugridpri <- -4+c(0:1000)*8/1000
fmupri <- dnorm(mugridpri,m,sqrt(1/c))
hist(mu,probability=TRUE,xlim=c(min(mugridpri[1],mugrid[1],d$x[1]),max(mugridpri[1001],mugrid[1001],d$x)),ylim=c(0,max(fmu,d$y)),nclass=20,xlab='mu',ylab='f',main='')
lines(d,type='l',lwd=2,col='green')
lines(mugrid,fmu,type='l',lwd=2,col='red')
lines(mugridpri,fmupri,type='l',lwd=2,col='blue')