Applied Statistics and Data Mining

文書の過去の版を表示しています。

draw.reliability.data.diagram <- function(x,sort=TRUE,...) {
  if(is.numeric(x)==TRUE) {
    n <- length(x)
    if(sort==TRUE) {
      x.time <- x[sort.list(x,decreasing=TRUE)]
      x.status <- rep(1,n)
      id <- c(1:n)
    } else {
      x.time <- x
      x.status <- rep(1,n)
      id <- c(1:n)
    }
  } else if (is.data.frame(x)==TRUE) {
    n <- dim(x)[1]
    if(sort==TRUE) {
      x.time <- x[sort.list(x[,1],decreasing=TRUE),1]
      x.status <- x[sort.list(x[,1],decreasing=TRUE),2]
      id <- row.names(x)[sort.list(x[,1],decreasing=TRUE)]
    } else {
      x.time <- x[,1]
      x.status <- x[,2]
      id <- row.names(x)
    }
  } else if (is.matrix(x)==TRUE) {
    n <- dim(x)[1]
    if(sort==TRUE) {
      x.time <- x[sort.list(x[,1],decreasing=TRUE),1]
      x.status <- x[sort.list(x[,1],decreasing=TRUE),2]
      id <- rawnames(x)[sort.list(x[,1],decreasing=TRUE)]
    } else {
      x.time <- x[,1]
      x.status <- x[,2]
      id <- rawnames(x)
    }
  }
  plot(x.time,c(1:n),
       xlim=c(0,max(x.time)),
       type="n",
       ylab="#",
       xlab="Lifetime data",yaxt="n")
  points(x.time,c(1:n),pch=c(1,4)[x.status+1],...)
  for( i in c(1:n) ) {
    lines(c(0,x.time[i]),c(i,i),...)
  }
}

compare.histgram.and.densities <- function(X,hist=TRUE) {
  require(MASS)
  est.weibull <- fitdistr(X$time,densfun="weibull")
  est.lnorm <- fitdistr(X$time,densfun="lognormal")
  est.gamma <- fitdistr(X$time,densfun="gamma")
  if(hist==TRUE) {
    hist(X$time,freq=FALSE,ylim=c(0,0.3),
         xlab="Lifetime")
    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),
       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))
}