adrec<-function(n,h,hd,yL,yU,...){
# Funktio adaptiiviseen hylkäyspoimintaan
# Muutettu viimeksi 3.3.2010
# n: generoitavien satunnaislukujen määrä
# h: tiheysfunktion logaritmi
# hd: tiheysfunktion logaritmin derivaatta
# yL: arvoalueen alaraja
# yU: arvoalueen yläraja 
# ...: parametrejä funkioille h ja hd


eps <- 1e-6
x0 <- uniroot(hd,lower=yL,upper=yU,...)$root  #Maximum point of h
hd2 <- (hd(x0+eps,...)-hd(x0-eps,...))/(2*eps) #Second derivative of h at maximimum

y <- c(x0-1/sqrt(-hd2),x0+1/sqrt(-hd2))

r <- numeric(n)
v <- 0
for(i in 1:n)
{
  reject<-TRUE
  while(reject)   
  {
     if(v==0){
       k <- length(y)
       a <- b <- p <- z<- numeric(k)          
       z[1]<-yL
       z<-c(z,yU)
       for(j in 1:(k-1)) 
          z[j+1]<-y[j]+(h(y[j],...)-h(y[j+1],...)+(y[j+1]-y[j])*hd(y[j+1],...))/(hd(y[j+1],...)-hd(y[j],...))
       for(j in 1:k){
          a[j]<-h(y[j],...)-y[j]*hd(y[j],...)
          b[j]<-hd(y[j],...)   
          p[j] <- exp(a[j])*(exp(b[j]*z[j+1])-exp(b[j]*z[j]))/b[j] 
       }     
     } 
      
     j <- sample(k,size=1,prob=p)             
     u <- 1/b[j]*log(exp(b[j]*z[j])+b[j]*p[j]/exp(a[j])*runif(1))
     v <- rbinom(1,size=1,prob=exp(h(u,...)-a[j]-b[j]*u))
     if(v==1){
       reject <- FALSE
       r[i]<-u
     }  
     else
       y<-sort(c(y,u)) 
   }
}
r
} 

# Kokeilu: generoidaan 1000 havaintoa normaalijakaumasta N(5,7)

hkoe<-function(x,mu,sigma2){
-0.5*log(sigma2)-0.5*(x-mu)^2/sigma2
}
hdkoe<-function(x,mu,sigma2){
-(x-mu)/sigma2
}

y<-adrec(1000,hkoe,hdkoe,-100,100,mu=5,sigma2=7)

qqnorm(y)
qqline(y)