R: How to Visualize an Integral

Question

I am working with the R programming language.

Below, here is some R code (adapted from here: R: Making a Transparent Plot) for making a 3D Bell Curve (Gaussian/Normal Distribution):

library(plotly)
library(mvtnorm)  


mu <- c(0, 0)  
Sigma <- matrix(c(1, 0, 0, 1), nrow=2)  

x <- seq(-5, 5, length.out = 50)
y <- seq(-5, 5, length.out = 50)
grid <- expand.grid(x=x, y=y)

grid$z <- apply(grid, 1, function(x) dmvnorm(x, mean=mu, sigma=Sigma))

fig <- plot_ly(data = grid, x = ~x, y = ~y, z = ~z, type = 'scatter3d', mode = 'lines',
               line = list(color = '#0066FF', width = 2)) %>%
    layout(title = "3D Normal Distribution Wireframe Plot",
           scene = list(xaxis = list(title = "X",
                                     gridcolor = 'rgb(255, 255, 255)',
                                     zerolinecolor = 'rgb(255, 255, 255)',
                                     showbackground = TRUE,
                                     backgroundcolor = 'rgb(230, 230,230)'),
                        yaxis = list(title = "Y",
                                     gridcolor = 'rgb(255, 255, 255)',
                                     zerolinecolor = 'rgb(255, 255, 255)',
                                     showbackground = TRUE,
                                     backgroundcolor = 'rgb(230, 230,230)'),
                        zaxis = list(title = "Z",
                                     gridcolor = 'rgb(255, 255, 255)',
                                     zerolinecolor = 'rgb(255, 255, 255)',
                                     showbackground = TRUE,
                                     backgroundcolor = 'rgb(230, 230,230)')),
           showlegend = FALSE)

fig

My Question: I want to "shade in" the volume corresponding to the integral of this function (i.e. the bell curve) between the regions X = -1 to 1 and Y = -1 to 1.

Conceptually, I think this means that I need to select every point in this diagram and check if this point lies within the boundary (-1,1) and (-1,1) : If yes, then shade else no.

Here is my attempt to do this:

library(plotly)
library(mvtnorm)  

mu <- c(0, 0)  
Sigma <- matrix(c(1, 0, 0, 1), nrow=2)  

x <- seq(-5, 5, length.out = 50)
y <- seq(-5, 5, length.out = 50)
grid <- expand.grid(x=x, y=y)

grid$z <- apply(grid, 1, function(x) dmvnorm(x, mean=mu, sigma=Sigma))

grid$color <- '#0066FF'

grid$color[grid$x >= -1 & grid$x <= 1 & grid$y >= -1 & grid$y <= 1] <- 'red'

fig <- plot_ly(data = grid, x = ~x, y = ~y, z = ~z, type = 'scatter3d', mode = 'lines',
               line = list(color = ~color, width = 2)) %>%
    layout(title = "3D Normal Distribution Wireframe Plot",
           scene = list(xaxis = list(title = "X",
                                     gridcolor = 'rgb(255, 255, 255)',
                                     zerolinecolor = 'rgb(255, 255, 255)',
                                     showbackground = TRUE,
                                     backgroundcolor = 'rgb(230, 230,230)'),
                        yaxis = list(title = "Y",
                                     gridcolor = 'rgb(255, 255, 255)',
                                     zerolinecolor = 'rgb(255, 255, 255)',
                                     showbackground = TRUE,
                                     backgroundcolor = 'rgb(230, 230,230)'),
                        zaxis = list(title = "Z",
                                     gridcolor = 'rgb(255, 255, 255)',
                                     zerolinecolor = 'rgb(255, 255, 255)',
                                     showbackground = TRUE,
                                     backgroundcolor = 'rgb(230, 230,230)')),
           showlegend = FALSE)


fig

But I am not sure if I have done this correctly. In general, is there a way to visualize the results of an integral within R?

Thanks!

References:

• https://community.plotly.com/t/fill-volume-under-the-surface/64944/2

Answer 1

It seems that you are looking for a shaded volume to be placed under the surface, representing its integral. To do this, we need a 3D grid where the voxels are turned on or off according to whether they are inside or outside the volume of interest. This is the same basic principle as outlined in this answer

For your example, we can set up the data as follows:

library(mvtnorm)
library(plotly)

mu <- c(0, 0)  
Sigma <- matrix(c(1, 0, 0, 1), nrow=2)  
grid <- expand.grid(x = seq(-5, 5, len = 50), y = seq(-5, 5, len = 50))

grid$z <- dmvnorm(grid, mu, Sigma)

grid3d <- expand.grid(x = seq(-1, 1, len = 100), 
                      y = seq(-1, 1, len = 100),
                      z = seq(0, 0.16, len = 100))

grid3d$value <- ifelse(grid3d[[3]] < dmvnorm(grid3d[1:2], mu, Sigma), 1, 0)

The plot would then be:

plot_ly() %>%
  add_trace(data = grid, x = ~x, y = ~y, z = ~z, type = 'scatter3d', 
            mode = 'lines', line = list(color = '#0066FF', width = 2),
            showlegend = FALSE) %>%
  add_trace(type = "volume", data = grid3d, x = ~x, y = ~y, z = ~z,
            value = ~value, opacity = 0.5, isomin = 0.1, showscale = FALSE,
            colorscale = list(c(0, 1), c("red", "red"))) %>%
  layout(title = "3D Normal Distribution Wireframe Plot",
         scene = list(xaxis = list(title = "X",
                                   gridcolor = 'rgb(255, 255, 255)',
                                   zerolinecolor = 'rgb(255, 255, 255)',
                                   showbackground = TRUE,
                                   backgroundcolor = 'rgb(230, 230,230)'),
                      yaxis = list(title = "Y",
                                   gridcolor = 'rgb(255, 255, 255)',
                                   zerolinecolor = 'rgb(255, 255, 255)',
                                   showbackground = TRUE,
                                   backgroundcolor = 'rgb(230, 230,230)'),
                      zaxis = list(title = "Z",
                                   gridcolor = 'rgb(255, 255, 255)',
                                   zerolinecolor = 'rgb(255, 255, 255)',
                                   showbackground = TRUE,
                                   backgroundcolor = 'rgb(230, 230,230)')))