Я использую связку, чтобы посмотреть на вероятность возникновения событий, основываясь на продолжительности и величине событий.Я могу создать контуры для интервалов повторения с наблюдаемыми и смоделированными данными в графике Base R, но я не могу понять, как воспроизвести в ggplot2.Почему бы просто не создавать графики в базовой графике и двигаться дальше, может быть, вам интересно?Потому что я включаю графики в краткий сводный отчет и хочу согласовать их с множеством других графиков.Ниже приведен пример кода.Я знаю, что использование местоположения, масштаба и формы для распределения GEV для создания случайных отклонений для получения того же распределения не является идеальным, но это лучший способ, который я мог бы придумать, чтобы создать несколько воспроизводимый пример, несмотря на плохую корреляциюв конце.В базе R контуры генерируются из матрицы моделируемых данных.Возможно ли это в ggplot2?
library(evd)
library(copula)
dur <- rgev(500, 2.854659, 2.170122, -0.007829)
mag <- rgev(500, 0.02482, 0.01996, 0.04603)
fDurGev <- fgev(dur)
fMagGev <- fgev(mag)
durVec <- dgev(dur, fDurGev[[1]][1], fDurGev[[1]][2], fDurGev[[1]][3])
magVec <- dgev(mag, fMagGev[[1]][1], fMagGev[[1]][2], fMagGev[[1]][3])
durMagMat <- as.matrix(cbind(duration = durVec, magnitude = magVec))
theta <- coef(fitCopula(claytonCopula(dim = 2), durMagMat, method = "itau"))
clayCop <- claytonCopula(theta, dim = 2)
fCopDurMag <- pCopula(durMagMat, clayCop)
copPts <- data.frame(duration = dur, magnitude = mag, copNEP = fCopDurMag,
copEP = (1 - fCopDurMag), copRI = (1 / fCopDurMag))
fSim <- seq(0.05, 0.99998, length.out = 1000)
quaDur <- qgev(fSim, fDurGev[[1]][1], fDurGev[[1]][2], fDurGev[[1]][3])
quaMag <- qgev(fSim, fMagGev[[1]][1], fMagGev[[1]][2], fMagGev[[1]][3])
expDurMagMat <- cbind(expand.grid(fSim, fSim)$Var1, expand.grid(fSim,
fSim)$Var2)
simPred <- pCopula(expDurMagMat, clayCop)
simPredMat <- matrix(simPred, 1000, 1000)
simDF <- data.frame(simDur = quaDur, simMag = quaMag, simPredMat)
rndPred <- data.frame(rCopula(5000, clayCop))
rndPred$rndDur <- qgev(rndPred[,1], fDurGev[[1]][1], fDurGev[[1]][2],
fDurGev[[1]][3])
rndPred$rndMag <- qgev(rndPred[,2], fMagGev[[1]][1], fMagGev[[1]][2],
fMagGev[[1]][3])
RI <- c(1.25, 2 ,5, 10, 20, 50, 100, 200, 500)
NEP <- 1 - (1 / RI)
plot(rndPred$rndDur, rndPred$rndMag, col = "light grey", cex = 0.5, xlab =
"Duration (time)", ylab = "Magnitude (x)")
points(copPts[,1], copPts[,2], col = "red", cex = 0.5)
contour(simDF$simDur, simDF$simMag, simPredMat, levels = NEP, labels = RI,
xaxs = 'i', yaxs = 'i', labcex = 0.6, lwd = 1, col = "black", add =
TRUE, method = "flattest", vfont = c("sans serif", "plain"))
А теперь для моей попытки воссоздать в ggplot2 (который не рисует контуры).
library(dplyr)
simDF <- data.frame(dur = expDurMagMat[, 1], mag = expDurMagMat[, 2], NEP = simPred)
simDF <- simDF %>%
dplyr::mutate(quaDur = qgev(NEP, fDurGev[[1]][1], fDurGev[[1]][2], fDurGev[[1]][3])) %>%
dplyr::mutate(quaMag = qgev(NEP, fMagGev[[1]][1], fMagGev[[1]][2], fMagGev[[1]][3]))
library(ggplot2)
ggplot(data = rndPred, aes(x = rndDur, y = rndMag)) +
geom_point(color = "light grey", alpha = 0.5) +
labs(x = "Duration (time)", y = "Magnitude (x)") +
geom_point(data = copPts, aes(x = duration, y = magnitude),
color = "red") +
geom_contour(data = simDF, aes(x = quaDur, y = quaMag, z = NEP),
inherit.aes = FALSE, breaks = NEP) +
theme_classic()
Спасибо всем, кто может помочь.