0

我正在尝试借助 R 中的 ggplot2 包构建用于可靠性分析的 Weibull 概率网格。我已经转换了数据，即时间和故障概率，使得两个参数 weibull 的累积分布函数符合直线网格中的线。

对于 y 轴，我编写了一个函数来转换显示百分比的双对数刻度。

现在我想在 x 轴上放置适当的中断和标签。单位只是时间，但轴是对数。对于较小的故障时间，我希望用整数标记和分解轴，例如像 c(100, 1000, 2000, 5000) 这样的向量。

对于更高的故障时间（时间> 5000），我想切换到对数刻度 10 ^(x)，x >= 4。

我在将 x 轴上的两件事结合起来时遇到了麻烦，希望有人可以帮助我解决这个问题。

下面是我的代码中“有趣”的部分。你可以看到我已经为对数 (10 ^(x)) 情况打破并标记了 x 轴。但正如我所说，我希望在一个轴上进行组合。
这是我到目前为止编码的内容。

library(ggplot2)
library(scales)
library(gridExtra)

set.seed(101)

# length of data vector has to be plugged in the following, i.e. how many obs are given  
data <- numeric(length = 5000)

for (sim in 1:length(data)) {
    data[sim] <- rweibull(n = 1, shape = 2.2, scale = 35000)
}

data <- data.frame(time = round(data, digits = 0), event = 1)


# Estimating failure probability for complete data with the inverse of the incomplete beta function 
# Referring to "Monte Carlo Pivotal Confidence Bounds for Weibull Analysis, with Implementations in R, p. 3 ff., Jurgen Symnyck Filip De Bal"
betaPlotPos <- function(r, n) {
    prob <- qbeta(.5, r, n - r + 1)
}

RankingData <- function(data) {
    n <- nrow(data)
    data$rank <- rank(data$time, ties.method = "first")
    data <- data[order(data$rank),]
    data <- cbind(data, beta.prob = betaPlotPos(data$rank, n))
}
data <- RankingData(data)

FTrans <- function(p) {
    log( - log(1 - p))
}

parameterOLS <- function(x, y) {
    b <- cov(log(x), FTrans(y)) / var(FTrans(y))
    a <- mean(log(x)) - b * mean(FTrans(y))
    beta.hat <- 1 / b
    eta.hat <- exp(a)
    rbind(beta.hat, eta.hat)
}
estimation <- data.frame((parameterOLS(data$time, data$beta.prob)))
colnames(estimation) <- "estimates"

rg <- function(t) {
    prob <- estimation[1, 1] * log(t) - estimation[1, 1] * log(estimation[2, 1])
}
data <- cbind(data, prob.points = rg(data$time))

scaleConverter <- function(x) {
    round((1 - (exp( - exp(x)))) * 100, digits = 3)
}

wb.grid <- ggplot(data = data, aes(x = time, y = FTrans(beta.prob))) +
             geom_point(shape = 1, color = "blue", size = 1.5) +
             geom_line(aes(x = time, y = prob.points), color = "red", size = 1.2) +
             scale_x_log10(breaks = trans_breaks("log10", function(x) 10 ^ as.integer(x)), labels = trans_format("log10", math_format(10 ^ .x))) +
             scale_y_continuous(breaks = if (data$beta.prob[length(data$beta.prob) / 2] <= .5) {
               c(-4.600149, -2.970195, -2.250367, -0.3665129, -.000328, 0.8340324, 1.097189, 1.52718, 2.220327)
               } else {
               c(-4.600149, -2.970195, -2.250367, -0.3665129)
               },
               labels = scaleConverter) +
             annotation_logticks(base = 10, sides = "b") +
             theme_bw() +
             theme(plot.title = element_text(size = 24, color = "black", face = "bold"), axis.text = element_text(size = rel(1)),
               axis.title = element_text(size = 20, color = "black"),
               panel.background = element_rect(color = "black", size = 2.5), panel.grid.major = element_line(color = "gray50", size = .5)) +
               xlab("Lifetime in kilometer [km]") +
               ylab("Unreliability in Percent [%]") +
               ggtitle("Two Parametric Weibull Distribution")

希望我的语言足以解释我的问题。

最好的，蒂姆

4

0 回答 0