差分

このページの2つのバージョン間の差分を表示します。

この比較画面へのリンク

両方とも前のリビジョン前のリビジョン
次のリビジョン		前のリビジョン
r:survival:code [2018/12/09 22:01] – watalur:survival:code [2020/10/20 16:16] (現在) – jin
行 1: 行 1:
 <code> <code>
-draw.reliability.data.diagram <- function(x,sort=TRUE,...) {+plot.reliability.data.diagram <- function(x,sort=TRUE,...) {
   if(is.numeric(x)==TRUE) {   if(is.numeric(x)==TRUE) {
     n <- length(x)     n <- length(x)
行 39: 行 39:
        type="n",        type="n",
        ylab="#",        ylab="#",
-       xlab="Lifetime data",yaxt="n")+       xlab="Lifetime data", 
 +       yaxt="n") 
 +  axis(2,1:n,labels=c(n:1))
   points(x.time,c(1:n),pch=c(1,4)[x.status+1],...)   points(x.time,c(1:n),pch=c(1,4)[x.status+1],...)
   for( i in c(1:n) ) {   for( i in c(1:n) ) {
行 46: 行 48:
 } }
  
-compare.histgram.and.densities <- function(X,hist=TRUE) {+compare.histgram.and.densities <- function(X,hist=TRUE,plot=TRUE) {
   require(MASS)   require(MASS)
   est.weibull <- fitdistr(X$time,densfun="weibull")   est.weibull <- fitdistr(X$time,densfun="weibull")
   est.lnorm <- fitdistr(X$time,densfun="lognormal")   est.lnorm <- fitdistr(X$time,densfun="lognormal")
   est.gamma <- fitdistr(X$time,densfun="gamma")   est.gamma <- fitdistr(X$time,densfun="gamma")
-  if(hist==TRUE) { +  if(plot==TRUE) { 
-    hist(X$time,freq=FALSE,ylim=c(0,0.3), +    if(hist==TRUE) { 
-         xlab="Lifetime") +      hist(X$time,freq=FALSE,ylim=c(0,0.3), 
-    plot(function(x){return(dlnorm(x,meanlog=est.lnorm$estimate[1], +           xlab="Lifetime") 
-                                   sdlog=est.lnorm$estimate[2]))},+      plot(function(x){return(dlnorm(x,meanlog=est.lnorm$estimate[1], 
 +                                     sdlog=est.lnorm$estimate[2]))}, 
 +           xlim=c(0,12), 
 +           col="blue", 
 +           add=TRUE) 
 +    } else { 
 +      plot(function(x){return(dlnorm(x,meanlog=est.lnorm$estimate[1], 
 +                                     sdlog=est.lnorm$estimate[2]))}, 
 +           xlim=c(0,12), 
 +           col="blue",xlab="Lifetime",ylab="Density") 
 +    } 
 +    plot(function(x){return(dweibull(x,shape=est.weibull$estimate[1], 
 +                                     scale=est.weibull$estimate[2]))},
          xlim=c(0,12),          xlim=c(0,12),
-         col="blue",+         col="red",
          add=TRUE)          add=TRUE)
-  } else { +    plot(function(x){return(dgamma(x,shape=est.gamma$estimate[1], 
-    plot(function(x){return(dlnorm(x,meanlog=est.lnorm$estimate[1], +                                   rate=est.gamma$estimate[2]))},
-                                   sdlog=est.lnorm$estimate[2]))},+
          xlim=c(0,12),          xlim=c(0,12),
-         col="blue",xlab="Lifetime",ylab="Density")+         col="green", 
 +         add=TRUE) 
 +    legend("topright", 
 +           col=c("blue","red","green"), 
 +           lty=c(1,1,1), 
 +           legend=c("lognormal","weibull","gamma"))
   }   }
-  plot(function(x){return(dweibull(x,shape=est.weibull$estimate[1], 
-                                   scale=est.weibull$estimate[2]))}, 
-       xlim=c(0,12), 
-       col="red", 
-       add=TRUE) 
-  plot(function(x){return(dgamma(x,shape=est.gamma$estimate[1], 
-                                 rate=est.gamma$estimate[2]))}, 
-       xlim=c(0,12), 
-       col="green", 
-       add=TRUE) 
-  legend("topright", 
-         col=c("blue","red","green"), 
-         lty=c(1,1,1), 
-         legend=c("lognormal","weibull","gamma")) 
   return(list(weibull=est.weibull,lnorm=est.lnorm,gamma=est.gamma))   return(list(weibull=est.weibull,lnorm=est.lnorm,gamma=est.gamma))
 } }