surfacePlots

Some of my students have been asking if there is a way to make interactive, 3D plots in R. Thanks to recent advances, the answer is yes! In this project, we will build a bivariate probability density function and [TO DO] explore its probabilities. Be sure to click (or mouseover) the graphs to get them to appear.

This document uses the following algebra function (exercise 9.12 in A Modern Introduction to Probability and Statistics by F.M. Dekking) \[f(x,y) = K(3x^2 + 8xy), \quad 0 \leq x \leq 1, \quad 0 \leq y \leq 2\]

Normalization

First we need to find that normalization constant K so that the volume formed represents one cubic unit (for later probability questions). Fortunately, double integrals can be done rather efficiently in Hans Werner Borchers’ pracma package (“Practical Numerical Math Functions”)

#install.packges("pracma") #if needed
library("pracma")
K <- 1

f <- function(x,y){K*(3*x^2 + 8*x*y)}
xa <- 0
xb <- 1
ya <- 0
yb <- 2

quad2d(f, xa, xb, ya, yb)

## [1] 10

We can see that without the normalization constant, the volume would be 10 cubic units, so we will need to normalize the function by that amount:

K <- 1 / quad2d(f, xa, xb, ya, yb)
f <- function(x,y){K*(3*x^2 + 8*x*y)}

Surface Plot

The plotly package in R (last updated November 12, 2016) creates interactive plots in R, including 3D surface plots!

Be sure to try out the code in your own RStudio IDE to feel the interactive elements.

#install.packages("plotly") #if needed
library("plotly")

x <- seq(xa, xb, 0.01)
y <- seq(ya, yb, 0.01)
z <- outer(x,y,f)
plot_ly(z = z, type = "surface") %>%
  layout(title = "f(x,y) = 0.1(3x^2 + 8xy)")

Displaying Probability

Here we will try to display \(P(2X \leq Y)\) with a plotly surface. This strategy entails making two separate surfaces and displaying them simultaneously.

# some space allocation
z_yes <- outer(x,y)
z_no <- outer(x,y)

for(i in x){
  for(j in y){
    z_yes[i,j] <- ifelse(2*i <= j, z[i,j], 0) #z = 0 outside desired region
  }
}
for(i in x){
  for(j in y){
    z_no[i,j] <- ifelse(2*i > j, z[i,j], 0) #z = 0 inside desired region
  }
}

plot_ly(showscale = FALSE) %>%
  add_surface(z = ~z_yes, opacity = 1.00) %>%
  add_surface(z = ~z_no, opacity = 0.32)

Conclusions

At the moment, the pracma package can compute double integrals easily over rectangular bounds (i.e. not variable bounds), and the plotly surface integrals act similarly. Later work should be able to garner volume slicing.