R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.1  (2005-06-20), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

>  # 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")
[1] "m,k,breakz"
>    print(c(m,k,breakz))
[1]  4.0000  1.5000 -0.4375
>  #
>  # 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")
[1] "ust,vstp,vstm"
>    print(c(ust,vstp,vstm))
[1]  1.7550547  0.9500833 -0.5931599
>  #
>  # 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
[1] 7237
>    mean(ygood)                # true is 3.348267
[1] 3.333835
>    var(ygood)
[1] 2.151553
>    sqrt(var(ygood)/(n*(n-1))) # standard error -- use correct n
[1] 0.0002026971
>  #  done
>   rm(list=ls())
>   q()