Как получить данные класса Uvtop продукта анализа POT-пакета в R

Question

Я использую пакет POT для выполнения определенных вычислений в R. Результаты анализа хранятся в объекте класса uvtop. Теперь я хотел бы экспортировать результат анализа, а не просто отобразить его в окне R.

Вот пример использования данных из этого пакета.

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)

С этим я получаю изображение ниже:

Хорошо, но я хочу экспортировать необходимые данные для воспроизведения графика уровня возврата (нижняя правая панель) в другом месте, то есть значения, представленные кружками, сплошной и обеими пунктирными линиями.

Answer 1

Чтобы получить данные, которые отображаются для уровня возврата, нам нужно углубиться в функцию retlev. По сути, я сделал все возможное, чтобы вырезать все графики и построить data.frame из необходимых точек.

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")

Наложение извлеченных данных против retlev графика. CI немного отличаются из-за выборки.

Answer 2

Если вы не хотите читать файл с другими приложениями, кроме R, просто сохраните его с:

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

или

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

Чтобы прочитать его обратно, загрузите библиотеку, которая сгенерировала данные, а затем используйте

load("myfilename.Rdata")

чтобы загрузить объект обратно в рабочее пространство или

mle <- readRDS("myfilename.Rds")

save сохраняет больше сред вместе с объектом, чем saveRDS, в зависимости от того, как реализована библиотека, это может иметь значение. Я бы предложил использовать save, если наборы данных не становятся слишком большими.