CODEWITHSUNDEEP
rplot
Pythonist
Pythonist

Reputation: 2679

How to get data of class uvtop product of a POT package analysis in R

I'm using the POT package to carry out certain calculations in R. The output of the analysis is stored in a object of class uvtop. Now, I'd like to export the result of the analysis, rather than just plot it within an R window.

Here it goes an example, using sample data from this package.

data(ardieres)
events1 <- clust(ardieres, u = 6, tim.cond = 8/365, clust.max = TRUE)
npy <- length(events1[,"obs"]) / (diff(range(ardieres[,"time"], na.rm
= TRUE)) - diff(ardieres[c(20945,20947),"time"]))
mle <- fitgpd(events1[,"obs"], thresh = 6, est = "mle")
par(mfrow=c(2,2))
plot(mle, npy = npy)

With this, I get the image below: test output

OK, but what I want is to export the necessary data to reproduce the Return Level Plot (bottom right panel) somewhere else, i.e. the values represented by circles, the solid and both dashed lines.

Upvotes: 1

Views: 143

Answers (2)

emilliman5
emilliman5

Reputation: 5966

To get the data that is plotted for Return Level, we have to dig into the retlev function. Basically, I did my best to strip all of the plotting out and construct a data.frame of the required points.

getRetLevData <- function (fitted, npy) {
  data  <- fitted$exceed
  loc   <- fitted$threshold[1]
  scale <- fitted$param["scale"]
  shape <- fitted$param["shape"]
  n     <- fitted$nat

  pot.fun <- function(T) {
    p <- rp2prob(T, npy)[, "prob"]
    return(qgpd(p, loc, scale, shape))
  }

  eps <- 10^(-3)

  if (!is.null(fitted$noy)){ 
    npy <- n/fitted$noy
    } else if (missing(npy)) {
    warning("Argument ``npy'' is missing. Setting it to 1.")
    npy <- 1
  }

  xlimsup    <- prob2rp((n - 0.35)/n, npy)[, "retper"]
  fittedLine <- pot.fun(seq(1/npy + eps, xlimsup, length.out = n))
  p_emp      <- (1:n - 0.35)/n
  xPoints    <- 1/(npy * (1 - p_emp))
  empPoints  <- sort(data)
  samp       <- rgpd(1000 * n, loc, scale, shape)
  samp       <- matrix(samp, n, 1000)
  samp       <- apply(samp, 2, sort)
  samp       <- apply(samp, 1, sort)
  ci_inf     <- samp[25, ]
  ci_sup     <- samp[975, ]

  rst <- data.frame(xPoints, fittedLine, empPoints, 
                    ci_inf, ci_sup)
}

x <- getRetLevData(mle, npy)
head(x)
#    fittedX fittedLine  xPoints empPoints   ci_inf   ci_sup
#1  1.001000   6.003716 1.011535      6.09 6.001557 6.239971
#2  3.891288  11.678503 1.029810      6.09 6.014536 6.363070
#3  6.781577  14.402517 1.048758      6.09 6.042065 6.470195
#4  9.671865  16.282306 1.068416      6.19 6.074348 6.583290
#5 12.562153  17.740710 1.088825      6.44 6.114193 6.684118
#6 15.452441  18.942354 1.110029      6.45 6.146098 6.812058
write.csv(x, "my_pot_results.csv")

An overlay of the extracted data vs retlev plot. The CI's are a bit different because of sampling.

enter image description here

Upvotes: 1

snaut
snaut

Reputation: 2535

If you don't want to read the file with other applications than R, just save it with:

save(mle, file="myfilename.Rdata")

or

saveRDS(mle, file="myfilename.Rds")

To read it back in, load the library that generated the data and then use 

load("myfilename.Rdata")

to load the object back to your workspace or 

mle <- readRDS("myfilename.Rds")

save saves more environments together with the object than saveRDS depending on how the library is implemented this might make a difference. I'd suggest to use save unless the datasets are getting too large.

Upvotes: 0

Related Questions