Home > Software design >  box muller equations returning NA values
box muller equations returning NA values

Time:09-28

  #generating 100 uniformly distributed numbers
  u1 <- runif(100,0,1)
  u2 <- runif(100,0,1)
  

  x1 <- function(x, y) {
    return(sqrt(-2 * log(x) * cos(2 * pi * y)))
  }
  x2 <- function(x, y) {
    return(sqrt(-2 * log(x) * sin(2 * pi * y)))
  }
  #applying x1
  x1_vals <- mapply(x1, u1, u2)
  #applying x2
  x2_vals <- mapply(x2, u1, u2)

Hi I want to write a box muller function, this is part of my attempt above, I'm trying to avoid using loops as much as possible(for loops especially) However, I keep getting NA values for my x1/x2 functions. I can't figure out whats wrong.

CodePudding user response:

Those functions are already vectorized. Don't need to wrap mapply around them. Just give the two vectors to x1 and x2

x1(u1,u2)
  [1]       NaN       NaN       NaN       NaN 0.3088330
  [6]       NaN 0.7866889       NaN       NaN 1.7102801
 [11] 2.1886770       NaN 1.5473627 1.0644378 0.8499208
snipped remaining 100 values

Vectorization is a feature of the R language. If the expressions in the function body can all take vectors and return vectors of equal length then no loop wrapper is needed.

You are getting NA's because the domain of the arguments to sin and cos are causing both positive and negative values to be given to sqrt.

CodePudding user response:

Here is the Box-Muller formula corrected.

x1 <- function(x, y) sqrt(-2 * log(x)) * cos(2 * pi * y)
x2 <- function(x, y) sqrt(-2 * log(x)) * sin(2 * pi * y)

u1 <- runif(1000,0,1)
u2 <- runif(1000,0,1)

# applying x1
x1_vals <- x1(u1, u2)
all(is.finite(x1_vals))
#> [1] TRUE

# applying x2
x2_vals <- x2(u1, u2)
all(is.finite(x2_vals))
#> [1] TRUE

old_par <- par(mfrow = c(1, 2))
hist(x1_vals, freq = FALSE)
curve(dnorm, from = -4, to = 4, add = TRUE)
hist(x2_vals, freq = FALSE)
curve(dnorm, from = -4, to = 4, add = TRUE)

par(old_par)

Created on 2022-09-28 with reprex v2.0.2

  •  Tags:  
  • r
  • Related