差分
このページの2つのバージョン間の差分を表示します。
両方とも前のリビジョン前のリビジョン次のリビジョン | 前のリビジョン | ||
r:survival:code [2018/12/09 22:00] – watalu | r:survival:code [2020/10/20 16:16] (現在) – jin | ||
---|---|---|---|
行 1: | 行 1: | ||
< | < | ||
- | draw.reliability.data.diagram <- function(x, | + | plot.reliability.data.diagram <- function(x, |
if(is.numeric(x)==TRUE) { | if(is.numeric(x)==TRUE) { | ||
n <- length(x) | n <- length(x) | ||
行 39: | 行 39: | ||
| | ||
| | ||
- | | + | |
+ | yaxt=" | ||
+ | axis(2, | ||
points(x.time, | points(x.time, | ||
for( i in c(1:n) ) { | for( i in c(1:n) ) { | ||
行 46: | 行 48: | ||
} | } | ||
- | compare.histgram.and.densities <- function(X, | + | compare.histgram.and.densities <- function(X, |
require(MASS) | require(MASS) | ||
est.weibull <- fitdistr(X$time, | est.weibull <- fitdistr(X$time, | ||
est.lnorm <- fitdistr(X$time, | est.lnorm <- fitdistr(X$time, | ||
est.gamma <- fitdistr(X$time, | est.gamma <- fitdistr(X$time, | ||
- | if(hist==TRUE) { | + | if(plot==TRUE) { |
- | hist(X$time, | + | |
- | | + | |
- | plot(function(x){return(dlnorm(x, | + | |
- | | + | plot(function(x){return(dlnorm(x, |
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | } | ||
+ | plot(function(x){return(dweibull(x, | ||
+ | | ||
| | ||
- | | + | |
| | ||
- | } else { | + | |
- | | + | rate=est.gamma$estimate[2]))}, |
- | sdlog=est.lnorm$estimate[2]))}, | + | |
| | ||
- | | + | |
+ | add=TRUE) | ||
+ | legend("topright", | ||
+ | col=c(" | ||
+ | | ||
+ | | ||
} | } | ||
- | plot(function(x){return(dweibull(x, | ||
- | | ||
- | | ||
- | | ||
- | | ||
- | plot(function(x){return(dgamma(x, | ||
- | | ||
- | | ||
- | | ||
- | | ||
- | legend(" | ||
- | | ||
- | | ||
- | | ||
return(list(weibull=est.weibull, | return(list(weibull=est.weibull, | ||
} | } | ||
- | </ | ||
+ | compare.KM.and.distributions <- function(fits, | ||
+ | require(survival) | ||
+ | km <- survfit(Surv(time, | ||
+ | plot(km, | ||
+ | | ||
+ | | ||
+ | plot(function(x) pweibull(x, | ||
+ | lower.tail=FALSE), | ||
+ | | ||
+ | plot(function(x) plnorm(x, | ||
+ | lower.tail=FALSE), | ||
+ | | ||
+ | plot(function(x) pgamma(x, | ||
+ | lower.tail=FALSE), | ||
+ | | ||
+ | legend(" | ||
+ | | ||
+ | paste(" | ||
+ | paste(" | ||
+ | paste(" | ||
+ | | ||
+ | } | ||
+ | |||
+ | compare.KM.and.distributions.2 <- function(data) { | ||
+ | require(survival) | ||
+ | km <- survfit(Surv(time, | ||
+ | fits <- compare.histgram.and.densities(data, | ||
+ | survreg.weibull <- survreg(Surv(time, | ||
+ | fits$weibull$estimate[1] <- 1/ | ||
+ | fits$weibull$estimate[2] <- exp(survreg.weibull$coefficients[1]) | ||
+ | survreg.lognormal <- survreg(Surv(time, | ||
+ | fits$lnorm$estimate[1] <- survreg.lognormal$coefficients[1] | ||
+ | fits$lnorm$estimate[2] <- survreg.lognormal$scale | ||
+ | plot(km, | ||
+ | | ||
+ | | ||
+ | plot(function(x) pweibull(x, | ||
+ | lower.tail=FALSE), | ||
+ | | ||
+ | plot(function(x) plnorm(x, | ||
+ | lower.tail=FALSE), | ||
+ | | ||
+ | legend(" | ||
+ | | ||
+ | paste(" | ||
+ | paste(" | ||
+ | | ||
+ | } | ||
+ | </ |