# gnbxml.r                November 2007
 #
 # generate normal random variable via Box-Muller method
 #
  set.seed(5151917)
 #
  N <- 1000         # how many to do
 #
  U <- runif(N)
  V <- runif(N)
  X <- cos(2*pi*U)*sqrt(-2*log(V))
  Y <- sin(2*pi*U)*sqrt(-2*log(V))      # note reverse roles
 # some simple analysis
  XY <- t(rbind(X,Y))          # get right shape
 # one way
  mean(XY)
 # second way
  apply(XY,2,mean)
 # one way
  cov(XY)
 # second way
  apply(XY,2,var)
 # third way
  C <- var(XY)
  C
 # test whether X and Y are correlated
  r <- C[1,2]/sqrt(C[1,1]*C[2,2])
  t <- r*sqrt(N/(1-r*r))
  c(r,t)
 #
 # analyze X's then Y's using Kolmogorov-Smirnov & Anderson-Darling
 #
  sX <- sort(X)
  w <- (2*(1:N)-1)/N
 # first K-S
  f <- pnorm(sX)
  Dnm <- max(f-((1:N)-1)/N)   # Dn-
  Dnp <- max((1:N)/N-f)       # Dn+
  c(Dnp,Dnm)                  # print them
  sqrt(N)*max(Dnp,Dnm)        # standardized statistic
 # now A-D
  w <- (2*(1:N)-1)/N
  lf <- pnorm(sX,log.p=TRUE)
 #              reverse the order
  lomf <- pnorm(rev(sX),lower.tail=F,log.p=TRUE)
  A2X <- -N - sum(w*(lf+lomf))
  A2X
 #                  repeat for Y 
  sY <- sort(Y)
 # first K-S
  f <- pnorm(sY)
  Dnm <- max(f-((1:N)-1)/N)   # Dn-
  Dnp <- max((1:N)/N-f)       # Dn+
  c(Dnp,Dnm)                  # print them
  sqrt(N)*max(Dnp,Dnm)        # standardized statistic
 # now A-D
  lf <- pnorm(sY,log.p=TRUE)
 #              reverse the order
  lomf <- pnorm(rev(sY),lower.tail=F,log.p=TRUE)
  A2Y <- -N - sum(w*(lf+lomf))
  A2Y
 # done
  rm(list=ls())
 q()