0

R の ggplot2 パッケージを使用して、信頼性分析用のワイブル確率グリッドを構築しようとしています。2 つのパラメトリック ワイブルの累積分布関数が直線に適合するように、時間と故障確率などのデータを変換しました。グリッドにライン。

y 軸については、パーセンテージを表示する二重対数スケールを変換する関数を作成しました。

ここで、適切なブレークとラベルを x 軸に合わせたいと思います。単位は時間ですが、軸は対数です。短い失敗時間の場合、軸にラベルを付けて整数で分割する必要があります。たとえば、c(100, 1000, 2000, 5000) のようなベクトルです。

失敗回数が多い場合 (時間 > 5000)、x >= 4 の対数スケール 10 ^(x) に切り替えたいと考えています。

x 軸上の 2 つのものを組み合わせるのに問題があり、誰かがこの問題を解決するのを手伝ってくれることを願っています。

以下に、私のコードの「興味深い」部分があります。対数 (10 ^(x)) の場合のためだけに、x 軸を分割してラベルを付けたことがわかります。しかし、私が言ったように、私は1つの軸に組み合わせを持ちたい.
これは私が今までコーディングしたものです。

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