逆Gamma函数在R语言里面怎么表示

Python032

逆Gamma函数在R语言里面怎么表示,第1张

逆Gamma函数R语言的自带包里没有,

一个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(α,β).