Search code examples
rfunction3dstatisticssample

Function that produces a sample of random vectors (X,Y,Z) in R

We have a spherical globe of radius 1, centre 1,1. There is spot randomly located on the globes surface. We are generating independent Unif(-1,1) random variables X,Y,Z which will be the coordinates of the random point. Divide (X,Y,Z) by √{X2 +Y2 +Z2} to get a point 1m from the centre of the globe. **

Write a function sample3d that produces a sample of random vectors (X, Y, Z), each of which is a point from a uniform distribution on the globe’s surface. Calling this function by the command sample3d(n) should produce an n × 3 array, where each row is a vector (X, Y, Z).

I have managed it up to ** but cannot create the sample3d, any help would be appreciated!

Solution

  • I had fun trying to come up with something helpful, here is what I got:

    First, I defined a function norm.2 which calculates the two norm of a vector:

    norm.2 <- function(x, na.rm){
  if(length(dim(x)) != 0){
    v1.logical <- ifelse(missing(na.rm), FALSE, TRUE)
    return(sqrt(colSums(x^2, na.rm = v1.logical)))
  }
  if(length(dim(x)) == 0){
    v1.logical <- ifelse(missing(na.rm), FALSE, TRUE)
    return(sqrt(sum(x^2, na.rm = v1.logical)))
  }
}

    Then, I proceeded to define a function sample3d which will give n points in the 3-dimensional space, all of which with norm 1, i.e. on the 2-sphere:

    sample3d <- function(n, Boundary){

  M <- ifelse(missing(Boundary), sqrt(2e+300), Boundary)

  x1 <- runif(n, min = -M, max = M)
  x2 <- runif(n, min = -M, max = M)
  x3 <- runif(n, min = -M, max = M)
  x <- t(cbind(x1,x2,x3))

  p <- t(t(x)/norm.2(x))

  df <- data.frame(t(p), stringsAsFactors = FALSE)

  return(df)  
}

    Here is a sample of the result:

    > head(sample3d(10000))
            x1         x2         x3
1  0.321159709 -0.5014622 -0.8033630
2  0.488181408  0.5547003 -0.6737852
3 -0.661576495 -0.4592729  0.5927773
4 -0.333447393  0.9331249 -0.1345016
5 -0.009070263  0.4267690  0.9043152
6 -0.375122328 -0.2393661 -0.8955373

    Now, using the plotly package, we can have some fun visualizing it:

    library(plotly)
dat <- sample3d(100000)
p <- plot_ly(dat, x = ~x1, y = ~x2, z = ~x3, color = norm.2(t(dat)), colors = c('#BF382A', '#0C4B8E')) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'X'),
                  yaxis = list(title = 'Y'),
                  zaxis = list(title = 'Z')))

    Yielding: