模拟 SAS 与 R 中伽玛分布的概率

Question

我有一些

SAS

代码，我正在尝试将其翻译为

R

。在某一时刻，

SAS

代码基本上通过以下方式模拟了概率

bhrsim

。我不确定如何将其翻译为

R

，但请在下面展示我的尝试。我认为我以前从未尝试过以这种方式模拟概率。首先是

SAS

代码：

data aa;

seed1    = 11111111 ;
seed2    = 33333333 ;
bhralpha = 2 ;
bhrbeta  = 2 ;

%let rep=50;

do i=1 to &rep;

call rangam(seed1,bhralpha,xgam1);
call rangam(seed2,bhrbeta,xgam2);

bhrsim=xgam1/(xgam1+xgam2);

output;
end;

proc print;
    var seed1 seed2 bhralpha bhrbeta xgam1 xgam2 bhrsim ;
run;

我有点惊讶

call rangam

函数没有同时使用

alpha

和

beta

参数。

这是上面

SAS

代码的输出，格式为

R

:

SAS.bhrsim <- c(0.43771, 0.28190, 0.47057, 0.30099, 0.25343,
                0.67331, 0.66465, 0.85924, 0.80344, 0.07354,
                0.16228, 0.92554, 0.52711, 0.56944, 0.10370,
                0.05858, 0.89869, 0.88515, 0.74498, 0.17853,
                0.71263, 0.49174, 0.21546, 0.42798, 0.26264,
                0.52674, 0.41922, 0.71119, 0.92044, 0.69456,
                0.33825, 0.55214, 0.42025, 0.93093, 0.20075,
                0.50655, 0.53586, 0.41479, 0.91006, 0.61604,
                0.71669, 0.72271, 0.91053, 0.73377, 0.50403,
                0.28722, 0.54455, 0.26749, 0.58494, 0.17943)

mean(SAS.bhrsim)
#[1] 0.5226472

这是我尝试将其翻译为

R

。我希望有人可以检查我的

R

代码是否正确并指出任何错误。特别是我没想到必须使用这一行：

bhrsim <- xgam1 / (xgam1 + xgam2)

in

R

，但这可能只是因为我以前从未这样做过。

set.seed(1234)

rep <- 50

bhralpha <- 2
bhrbeta  <- 2

xgam1 <- rgamma(rep, shape = bhralpha, scale = bhrbeta)
xgam2 <- rgamma(rep, shape = bhralpha, scale = bhrbeta)

bhrsim <- xgam1 / (xgam1 + xgam2)

mean(xgam1)
#[1] 3.698261

mean(xgam2)
#[1] 4.427575

mean(bhrsim)
#[1] 0.4477643

R.bhrsim <- c(0.34655843, 0.71945276, 0.20861699, 0.30611308, 0.55605485,
              0.55064680, 0.69470936, 0.45576462, 0.30185399, 0.10438540,
              0.48265790, 0.75655704, 0.42791822, 0.26072882, 0.45523966,
              0.45975115, 0.21874035, 0.22447877, 0.34634855, 0.38155777,
              0.11296304, 0.48800715, 0.37407339, 0.50943619, 0.71306467,
              0.85602900, 0.88723711, 0.05755638, 0.63186913, 0.13553985,
              0.67087093, 0.26015097, 0.26414524, 0.14128396, 0.71854283,
              0.77900243, 0.53769594, 0.46026569, 0.09196446, 0.67574945,
              0.24611917, 0.54923784, 0.60613130, 0.67733881, 0.22317935,
              0.64392641, 0.43528504, 0.44700128, 0.55522003, 0.38119564)

Answer 1

根据调用 rangam() 的 SAS 文档，call rangam(seed, a, x)

 使用随机种子 

a 生成一个 Gamma 随机偏差，其中 shape

 等于 

seed

，并将结果分配给 

x

（并更新种子 

internally — 使用相同的 call rangam(...)

 参数第二次调用 

seed

 使用 

RNG 流中的新值（这个习惯用法对我来说似乎很奇怪，但我知道什么？）。所以 R 等效项是set.seed(seed)

（一次，在循环开始时）；

x <- rgamma(1, a)

（R 允许您通过指定其他参数来设置速率或比例！= 1；第一个参数是请求的偏差数量。）

所以我相信 R 等价物是

set.seed(111111)
xgam1 <- rgamma(50, shape = 2)
xgam2 <- rgamma(50, shape = 2)
bhrsim <- xgam1 / (xgam1 + xgam2)
mean(bhrsim)
## [1] 0.5579772