# ghutmm.r                   November 2007, 2008
 # compute bounds for ratio of uniforms method
 #  for generating from multivariate Huber distribution
 #
 #  boundary functions and derivatives
 #
 # first one is u(z)
  u <- function(z) {
     if( z >  breakz ) ul <- m*log(1+z)-m*(1+z)+k*k/2
      else ul <- m*log(1+z)-(m*m/(k*k))*(1+z)*(1+z)/2
     u <- exp((ul+m)/2) }               # sqrt(f(z))
 #
 # second is up(z)=u'(z), the derivative of u(z)
  up <- function(z) {
     if( z > breakz ) up <- m/(1+z) - m
      else up <- m/(1+z) - (m*m/(k*k))*(1+z) }
 #
 # next is v(z)
  v <- function(z) {
     if( z >  breakz ) ul <- m*log(1+z)-m*(1+z)+k*k/2
      else ul <- m*log(1+z)-(m*m/(k*k))*(1+z)*(1+z)/2
     v <- z*exp((ul+m)/2) }             # z*sqrt(f(z))
 #
 # last is vp(z)=v'(z), the derivative of v(z)
  vp <- function(z) {
     if( z > breakz ) vp <- m/(1+z) - m
      else vp <- m/(1+z) - (m*m/(k*k))*(1+z) 
     vp <- 1/z + vp/2 }
 #
 # end of functions
 #
 # set parameters
 #
   m <- 4                     # dimension - 1 so d = 5
   k <- 1.5                   # Huber cutoff
   breakz <- k*k/m - 1        # cutoff for transformed variable
   print("m,k,breakz")
   print(c(m,k,breakz))
 #
 # find boundary box
 #
 #                         limits of search
   ustar <- uniroot(f=up, lower=-.95, upper=(1+sqrt(1+2*m)) )
   ust <- u(ustar$root)
   vstarp <- uniroot(f=vp, lower=.02, upper=2*(1+sqrt(1+2*m)) )
   vstp <- v(vstarp$root)
   vstarm <- uniroot(f=vp, lower=-.95, upper=-.02 )
   vstm <- v(vstarm$root)
   print("ust,vstp,vstm")
   print(c(ust,vstp,vstm))
 #
 # got the box, now generate some
 #
   set.seed(5151917)
   N <- 10000
 # R of U on transformed variable
   u <- ust*runif(N)                    # uniform( 0, ust )
   v <- ( (vstp-vstm)*runif(N) + vstm ) # uniform( vstm, vstp )
   UV <- t(matrix(c(u,v),N,2))
 #
 # define function for testing
 #
   rofu <- function(uv) {
   U <- uv[1]
   V <- uv[2]
   z <- V/U
   if( z > -1 ) w <- m*log(1+z)
    else w <- -1e100                    # instead of log(0)
   y <- m*(1+z)/k
   if( y < k ) w <- w - y*y/2
    else w <- w - k*y + k*k/2
   if( U > exp((w+m)/2) ) rofu <- 0     # reject -- returns a zero
    else  rofu <- m*(1+z)/k }           # accept -- returns positive
 #
 # now apply function
 #
   y <- apply(UV,2,rofu)      # this is what I want
 #
 # only so many are good
 #
   ygood <- y[(y>0)]
   n <- length(ygood)
   n                          # number accepted
   mean(ygood)                # true is 3.348267
   var(ygood)
   sqrt(var(ygood)/(n*(n-1))) # standard error -- use correct n
 #  done
  rm(list=ls())
  q()