==== ワイブル確率紙プロット ====
=== weibull.probability.plot() ===
ワイブル確率紙にデータセットを、プロットして最小二乗法で直線を当てはめる関数。
複数のデータセットを同時にプロットすることができる。
weibull.probability.plot <- function(data=data.df,
log.buff=c(0.2, 0.2),
points.pch=NULL,
points.color=NULL,
lines.color=NULL,
color="grey",
trimming=NULL)
{
x <- data[,1]
x.condition <- data[,2]
log.x.min <- min(log10(x))
log.x.max <- max(log10(x))
data.condition <- unique(x.condition)
weibull.prob.transform <- function(p) { # ワイブル分布の累積分布関数を直線に変更するよう確率の尺度を変更
return(log10(log10(1/(1-p))))
}
log.axis <- function(z) # 対数軸を描く関数
{
z <- floor(log10(z)) # 対数にしたときの整数部
log.min <- min(z) # 最小値
z2 <- 1:10*10^log.min # 値の範囲をカバーするように
n <- max(z)-log.min # 10 倍しながら順次，右の位置に目盛りを描く
z2 <- rep(z2, n+1)*10^rep(0:n, each=10) # 対数目盛り位置の数値
log.z2 <- log10(z2) # 目盛りを描く位置
axis(1, at=log.z2, labels=z2) # log.z2 の位置に，z2 という数値を描く
abline(v=log.z2, col=color) # 垂直格子線を描く
}
max.n <- 0
for( i in data.condition ) {
log.x <- log10(sort(x[x.condition==i]))
n <- length(log.x)
max.n <- max(max.n, n)
}
log.x <- log10(sort(x))
log.x <- append(min(log.x)-log.buff[1], log.x)
log.x <- append(log.x, max(log.x)+log.buff[2])
# print(max(log.x))
# print(min(log.x))
# print(max.n)
probs <- c(10^(-10:0), 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99, 99.9, 99.99, 99.999)/100
probs <- probs[probs > 0.1/max.n]
probs.transformed <- weibull.prob.transform(probs)
plot(c(log.x[1], log.x[length(log.x)]),
c(probs.transformed[1], probs.transformed[length(probs.transformed)]),
type="n", xaxt="n", yaxt="n",
xlab="Observed Value",
ylab="Cumulative Percent",
main="Weibull Probability Paper")
abline(h=probs.transformed, col="grey") # 水平の格子線
log.axis(x) # 横軸を描く
axis(2, at=probs.transformed, labels=probs) # 縦軸を描く
for( i in data.condition) {
log.x <- log10(sort(x[x.condition==i]))
n <- length(log.x)
y <- 1:n/(n+1)
k <- match(i, data.condition)
if(!is.null(points.pch)&&!is.null(lines.color)) {
points(log.x, weibull.prob.transform(y), pch=points.pch[k], col=points.color[k])
} else {
points(log.x, weibull.prob.transform(y), pch=20)
}
if( (!is.null(trimming)) & (length(trimming)==2) ) {
i.min <- round(n*trimming[1])
i.max <- round(n*trimming[2])
} else {
i.min <- 1
i.max <- n
}
y.to.fit <- y[i.min:i.max]
log.x.to.fit <- log.x[i.min:i.max]
lm.fit <- lm(weibull.prob.transform(y.to.fit)~log.x.to.fit)
print(summary(lm.fit)$coefficients)
print(summary(lm.fit)$r.squared)
if(!is.null(lines.color)) {
abline(lm.fit, col=lines.color[k])
} else {
abline(lm.fit)
}
}
}
使い方は
wb.data <- data.frame(time=c(1,2,3,4,5,6,7,8,9,10), condition=c("A","A","B","A","B","B","A","B","A","B"))
weibull.probability.plot(data=wb.data,
log.buff=c(0.2, 0.2),
points.pch=NULL,
points.color=NULL,
lines.color=NULL,
color="grey",
trimming=NULL)
など。
=== weibull.probability.plot.2() ===
ワイブル確率紙に複数のデータセットを、同時にプロットして最小二乗法で直線を当てはめる関数。
複数のデータセットを同時にプロットすることができる上、打ち切りデータにも対応している。
weibull.probability.plot.2 <- function(data=data.df,
log.buff=c(0.2, 0.2),
points.pch=NULL,
points.color=NULL,
lines.color=NULL,
color="grey",
trimming=NULL,
debugging=FALSE,
xlim=NULL)
{
x <- data[,1]
x.condition <- data[,2]
x.status <- data[,3]
# x <- data$time
# x.condition <- data$condition
# x.status <- data$status
log.x.min <- min(log10(x))
log.x.max <- max(log10(x))
data.condition <- unique(sort(x.condition))
weibull.prob.transform <- function(p) { # ワイブル分布の累積分布関数を直線に変更するよう確率の尺度を変更
return(log10(log10(1/(1-p))))
# return(log10(p))
# return(p)
}
weibull.prob.transform2 <- function(p) { # ワイブル分布の累積分布関数を直線に変更するよう確率の尺度を変更
# return(log10(log10(1/(1-p))))
return(log10(p))
}
log.axis <- function(z) # 対数軸を描く関数
{
z <- floor(log10(z)) # 対数にしたときの整数部
log.min <- min(z) # 最小値
z2 <- 1:10*10^log.min # 値の範囲をカバーするように
n <- max(z)-log.min # 10 倍しながら順次，右の位置に目盛りを描く
z2 <- rep(z2, n+1)*10^rep(0:n, each=10) # 対数目盛り位置の数値
log.z2 <- log10(z2) # 目盛りを描く位置
axis(1, at=log.z2, labels=z2) # log.z2 の位置に，z2 という数値を描く
abline(v=log.z2, col=color) # 垂直格子線を描く
}
max.n <- 0
for( i in data.condition ) {
log.x <- log10(sort(x[x.condition==i]))
n <- length(log.x)
max.n <- max(max.n, n)
}
log.x <- log10(sort(x))
log.x <- append(min(log.x)-log.buff[1], log.x)
log.x <- append(log.x, max(log.x)+log.buff[2])
# print(max(log.x))
# print(min(log.x))
# print(max.n)
probs <- c(10^(-10:0), 5, 10, 20, 30, 40, 50, 60, 70, 80, 90, 95, 99, 99.9, 99.99, 99.999)
# probs <- c(10^(-10:0),2:10,20,30,40,50,100,200,500)
probs <- probs[probs > 10/max.n]
probs.transformed <- weibull.prob.transform(probs/100)
# probs.transformed <- weibull.prob.transform2(probs/100)
if(is.null(xlim)){
plot(c(log.x[1], log.x[length(log.x)]),
c(probs.transformed[1], probs.transformed[length(probs.transformed)]),
type="n",
xaxt="n",
yaxt="n",
xlab="Observed Value",
ylab="F(t) %",
main="Weibull Probability Paper")
# main="Cumulative Hazard Paper")
}else{
plot(c(log.x[1], log.x[length(log.x)]),
c(probs.transformed[1], probs.transformed[length(probs.transformed)]),
type="n",
xaxt="n",
yaxt="n",
xlab="Observed Value",
ylab="F(t) %",
main="Weibull Probability Paper",
xlim=xlim)
# main="Cumulative Hazard Paper")
}
abline(h=probs.transformed, col="grey") # 水平の格子線
log.axis(x) # 横軸を描く
axis(2, at=probs.transformed, labels=probs) # 縦軸を描く
# i <- data.condition[1]
for( i in data.condition) {
if(debugging==TRUE) print(i);
log.x <- log10(x[x.condition==i])
st <- x.status[x.condition==i]
x.sort <- sort.list(log.x)
log.x <- log.x[x.sort]
st <- st[x.sort]
n <- length(log.x)
hazard <- (st/c(length(st):1))
# Hazard <- cumsum(hazard) * 100
Hazard <- cumsum(hazard)
# y <- { 1-exp(-Hazard/100) }
y <- 1-exp(-Hazard)
# y <- Hazard
y <- y[st==1] # statusが1を抽出
log.x <- log.x[st==1] # statusが1を抽出
k <- match(i, data.condition)
if(!is.null(points.pch)&&!is.null(lines.color)) {
points(log.x, weibull.prob.transform(y), pch=points.pch[k], col=points.color[k])
} else {
points(log.x, weibull.prob.transform(y), pch=20,col=weibull.color(i))
# points(log.x, weibull.prob.transform2(y), pch=20)
}
if( (!is.null(trimming)) & (length(trimming)==2) ) {
i.min <- round(n*trimming[1])
i.max <- round(n*trimming[2])
} else {
i.min <- 1
i.max <- length(log.x)
}
y.to.fit <- y[i.min:i.max]
log.x.to.fit <- log.x[i.min:i.max]
lm.fit <- lm(weibull.prob.transform(y.to.fit)~log.x.to.fit) # for plot
# lm.fit <- lm(weibull.prob.transform2(y.to.fit)~log.x.to.fit) # for plot
# lm.fit <- lm(y.to.fit~log.x.to.fit) # for plot
print(lm.fit)
print(summary(lm.fit)$coefficients)
print(summary(lm.fit)$r.squared)
if(!is.null(lines.color)) {
abline(lm.fit, col=lines.color[k])
} else {
abline(lm.fit, col=weibull.color(i))
}
# for m,eta
y <- Hazard
y <- y[st==1] # statusが1を抽出
y.to.fit <- y[i.min:i.max]
lm.fit <- lm(weibull.prob.transform2(y.to.fit)~log.x.to.fit)
m <- summary(lm.fit)$coefficients[2]
# eta <- 10^((log10(100)-summary(lm.fit)$coefficients[1])/m)
eta<- 10^(-summary(lm.fit)$coefficients[1]/m)
print(paste("m",m,"eta",sprintf("%7.3e",eta)))
mu <- eta * gamma(1 + 1/m)
sigma <- eta * ( gamma(1 + 2/m) - (gamma(1 + 1/m))^2 )^(1/2)
print(paste("mu",sprintf("%7.3e", mu),"sigma",sprintf("%7.3e",sigma)))
cat("B10:"); cat( sprintf("%7.3e", eta*(-log(1-0.1))^(1/m) ));
cat(" B1:"); cat( sprintf("%7.3e", eta*(-log(1-0.01))^(1/m) ));
cat(" B0.1:"); cat( sprintf("%7.3e", eta*(-log(1-0.001))^(1/m) ));
cat(" B0.01:"); cat( sprintf("%7.3e", eta*(-log(1-0.0001))^(1/m) )); cat("\n")
}
}
使い方は次の通り。
wb.data.2 <- data.frame(time=c(1,3,4,4,5,5,5,2,2,4),
condition=c(1:10)/c(1:10), status=c(1,1,1,1,1,1,1,1,1,1))
weibull.probability.plot.2(data=wb.data.2,
log.buff=c(0.2, 0.2),
points.pch=NULL,
points.color=NULL,
lines.color=NULL,
color="grey",
trimming=NULL)
weibull.color()
weibull.color <- function(i=i)
{
if(i == 1){
return("blue")
}
if(i == 2){
return("red")
}
if(i == 3){
return("black")
}
if(i == 4){
return("darkturquoise")
}
if(i == 5){
return("darkorange")
}
if(i == 6){
return("dimgrey")
}
}