一个R包:SuppDists里有
dinvGauss,pinvGauss,qinvGauss,rinvGauss
详情见帮助文件。
如有疑问,请追问!
R里的伪随机数怎么取的不得而知,但逆变换法应该是在分布函数已知的情况下最方便的做法吧。
我们从最简单的指数分布来测试吧。方法1用逆变换,方法2用伪随机也就是R里的built-in.
最后比较每种方法和各自,还有和对方的最大绝对值差值的分布。
方法1:
invF <- function(x){ log(1/(1-x))}D <- vector()
for (i in 1:100){
a <- runif(1e4) # 两组1万个0-1均匀分布随机数
b <- runif(1e4)
ran1 <- sapply(a, invF)
ran2 <- sapply(b, invF)
D <- c(D, max(abs(ran1 - ran2)))
}
summary(D)
Min. 1st Qu. Median Mean 3rd Qu. Max.
7.651 9.100 9.762 9.875 10.540 13.940
方法2:
D <- vector()for (i in 1:100){
a <- rexp(1e4, rate = 1) # 两组1万个参数为1的指数分布随机数
b <- rexp(1e4, rate = 1)
D <- c(D, max(abs(a - b)))
}
summary(D)
Min. 1st Qu. Median Mean 3rd Qu. Max.
7.765 8.977 9.574 9.727 10.590 12.540
两种方法混合:a是逆变换,b是伪随机
invF <- function(x){ log(1/(1-x))}D <- vector()
for (i in 1:100){
a <- runif(1e4) # 1万个0-1均匀分布随机数
b <- rexp(1e4, rate = 1) # 1万个参数为1的指数分布随机数
ran1 <- sapply(a, invF)
D <- c(D, max(abs(ran1 - b)))
}
summary(D)
Min. 1st Qu. Median Mean 3rd Qu. Max.
7.679 8.701 9.385 9.536 10.150 13.610
有没有发现。。根本没什么不同。
WinBugs中没有直接生成逆伽马分布的函数,需要利用伽马分布和逆伽马分布的关系,若X~Gamma(alpha, beta),则1/X ~ Inverse Gamma(alpha, beta)。
t ~ IG(0.001, 1000)转换成代码:
#
tt ~ dgamma(1.0E-3,1.0E3)
t <- 1 / tt
#
伽玛分布(Gamma Distribution)是统计学的一种连续概率函数。Gamma分布中的参数α称为形状参
数(shape parameter),β称为尺度参数(scale parameter)。
实验定义:
假设随机变量X为 等到第α件事发生所需之等候时间
当两随机变量服从Gamma分布,互相独立,且单位时间内频率相同时,Gamma分布具有加成性
数学表达式
若随机变量X具有概率密度
其中α>0,β>0,则称随机变量X服从参数α,β的伽马分布,记作G(α,β).