7.2.10 威布尔分布(Weibull)
威布尔分布是指数分布的一个推广,此时“失效率”(正确的说法是“危险率”——在生存分析中会经常遇到)不是一个常量。用它来模型化某个机器生存时间(没有失败的时间)。(这里“失效率”就是随着机器使用年代的增加,出现问题的概率。)
密度函数的形式为:exp(-a t^b).
%G
curve(dexp(x), xlim=c(0,3), ylim=c(0,2))
curve(dweibull(x,1), lty=3, lwd=3, add=T)
curve(dweibull(x,2), col='red', add=T)
curve(dweibull(x,.8), col='blue', add=T)
title(main="Weibull Probability Distribution Function")
legend(par('usr')[2], par('usr')[4], xjust=1,
c('Exponential', 'Weibull, shape=1',
'Weibull, shape=2', 'Weibull, shape=.8'),
lwd=c(1,3,1,1),
lty=c(1,3,1,1),
col=c(par("fg"), par("fg"), 'red', 'blue'))
%--
7.2.11伽玛分布( Gamma )
它是独立的指数分布随机变量的和,也是指数分布的推广,通常用来描述生存时间(例如,一个可靠机器(的生存时间),如果遇到三个连续的问题就被程序化停止,那么每一个问题(出现的时间)都可以用指数分布律来描述)。
泊松过程中到达时间的分布就是伽玛分布。
%G
curve( dgamma(x,1,1), xlim=c(0,5) )
curve( dgamma(x,2,1), add=T, col='red' )
curve( dgamma(x,3,1), add=T, col='green' )
curve( dgamma(x,4,1), add=T, col='blue' )
curve( dgamma(x,5,1), add=T, col='orange' )
title(main="Gamma probability distribution function")
legend(par('usr')[2], par('usr')[4], xjust=1,
c('k=1 (Exponential distribution)', 'k=2', 'k=3', 'k=4', 'k=5'),
lwd=1, lty=1,
col=c(par('fg'), 'red', 'green', 'blue', 'orange') )
%--
%G
n <- 500
x1 <- rexp(n,17)
x2 <- rexp(n,17)
x3 <- rexp(n,17)
x <- x1 + x2 + x3
# Simpler, but less readable:
# k <- 3
# x <- drop(apply( matrix( rexp(n*k,17), nr=n, nc=k ), 1, sum))
y <- qgamma(ppoints(n),3,17)
plot( sort(x) ~ sort(y), log='xy' )
abline(0,1, col='red')
title("Comparision: gamma distribution and sum of exponential r.v.")
%--
可以参考:
http://www.math.uah.edu/statold/special/special3.html
7.2.12 Beta 分布
在许多参考书上都可以看到这个定义(不太具有启发性,稍后我将给出更直观的定义):
如果X和T 是独立的随机变量,分别服从参数为(a,r)和(b,r)的Gamma分布, 则X/(X+Y)服从参数为(a,b)的Beta分布。
%G
curve( dbeta(x,1,1), xlim=c(0,1), ylim=c(0,4) )
curve( dbeta(x,2,1), add=T, col='red' )
curve( dbeta(x,3,1), add=T, col='green' )
curve( dbeta(x,4,1), add=T, col='blue' )
curve( dbeta(x,2,2), add=T, lty=2, lwd=2, col='red' )
curve( dbeta(x,3,2), add=T, lty=2, lwd=2, col='green' )
curve( dbeta(x,4,2), add=T, lty=2, lwd=2, col='blue' )
curve( dbeta(x,2,3), add=T, lty=3, lwd=3, col='red' )
curve( dbeta(x,3,3), add=T, lty=3, lwd=3, col='green' )
curve( dbeta(x,4,3), add=T, lty=3, lwd=3, col='blue' )
title(main="Beta distribution")
legend(par('usr')[1], par('usr')[4], xjust=0,
c('(1,1)', '(2,1)', '(3,1)', '(4,1)',
'(2,2)', '(3,2)', '(4,2)',
'(2,3)', '(3,3)', '(4,3)' ),
lwd=1, #c(1,1,1,1, 2,2,2, 3,3,3),
lty=c(1,1,1,1, 2,2,2, 3,3,3),
col=c(par('fg'), 'red', 'green', 'blue',
'red', 'green', 'blue',
'red', 'green', 'blue' ))
%--
如果X1,X2,...,Xn是独立的随机变量,服从均匀分布U(0,1),则 max(X1,X2,...,Xn)服从参数为(n,1)的beta分布。
%G
N <- 500
n <- 5
y <- drop(apply( matrix( runif(n*N), nr=N, nc=n), 1, max ))
x <- qbeta(ppoints(N), n, 1)
plot( sort(y) ~ x )
abline(0,1, col='red')
title("Order statistic and Beta distribution")
%--
其它的次序统计量(X1,X2,...,Xn中第k大元素)也服从beta分布(参数不同而已)。
%G
N <- 500
n <- 5
k <- 3
y <- drop(apply( matrix( runif(n*N), nr=n, nc=N), 2, sort )[n-k,])
x <- qbeta(ppoints(N), n-k, k+1) # Exercice: Where do those
# coefficients come from?
plot( sort(y) ~ x )
abline(0,1, col='red')
title("Order statistics and Beta distribution")
%--