CODEWITHSUNDEEP
rplot3d
LaTeXFan
LaTeXFan

Reputation: 1231

Plot a generic surface and contour in R

I have the following data

var.asym <- function(alpha1, alpha2, xi, beta, n){
  term11 <- alpha1*(1-alpha1)^(2*xi-1)
  term12 <- alpha1*(1-alpha1)^(xi-1)*(1-alpha2)^xi
  term22 <- alpha2*(1-alpha2)^(2*xi-1)
  Sigma <- matrix(c(term11, term12, term12, term22), nrow=2, byrow=TRUE)
  Sigma*beta^2/n
}

mop.jacob.inv <- function(alpha1, alpha2, xi, beta){
  term11 <- -qgpd(alpha1, xi, 0, beta)/xi - beta*(1-alpha1)^xi*log(1-alpha1)/xi
  term12 <- qgpd(alpha1, xi, 0, beta)/beta
  term21 <- -qgpd(alpha2, xi, 0, beta)/xi - beta*(1-alpha2)^xi*log(1-alpha2)/xi
  term22 <- qgpd(alpha2, xi, 0, beta)/beta
  jacob <- matrix(c(term11, term12, term21, term22), nrow=2, byrow=TRUE)
  jacob.inv <- solve(jacob)
  jacob.inv
}

var.asym2 <- function(alpha1, alpha2) var.asym(alpha1, alpha2, 0.2, 1, 1000)
mop.jacob.inv2 <- function(alpha1, alpha2) mop.jacob.inv(alpha1, alpha2, 0.2, 1)
object <- function(alpha1, alpha2){
  term1 <- mop.jacob.inv2(alpha1, alpha2)%*%var.asym2(alpha1, alpha2)%*%t(mop.jacob.inv2(alpha1, alpha2))
  sum(diag(term1))
}

x <- seq(0.01, 0.98, by=0.01)
y <- seq(x[1]+0.01, 0.99, by=0.01)
xy <- cbind(rep(x[1], length(x)), y)

for(i in 2:length(x)){
  y <- seq(x[i]+0.01, 0.99, by=0.01)
  xy <- rbind(xy, cbind(rep(x[i], length(x)-i+1), y))
}

object.xy <- rep(0, 4851)

for(i in 1:4851){
  object.xy[i] <- object(xy[i, 1], xy[i, 2])
}

Now I want to plot a surface of (xy[, 1], xy[, 2], object.xy). Is there a way to do so in R? I tried persp and contour function but it did not seem to be appropriate for this case since they both require increasing sequences x and y. I guess a more general question would be how to make contour plot when we are given a sequence of triplets (x, y, z).

Upvotes: 0

Views: 211

Answers (1)

bramtayl
bramtayl

Reputation: 4024

library(dplyr)
library(tidyr)
library(magrittr)

long_data = 
  data.frame(
    x = xy[,1] %>% round(2),
    y = xy[,2] %>% round(2),
    z = object.xy)

wide_data = 
  long_data %>%
  spread(x, z)

y = wide_data$y
wide_data %<>% select(-y)
x = names(wide_data) %>% as.numeric
z = wide_data %>% as.matrix
persp(x, y, z)
contour(x, y, z)

Dunno why the round helps, but it does. The reshape was necessary to build a matrix from x, y, z data. Note that the contour lines coalesce into a black dot because of the huge narrow peak in the data.

Upvotes: 1

Related Questions