模拟 SAS 与 R 中伽马分布的概率-解网

问：

我有一些代码正在尝试翻译成.在某一时刻，代码基本上模拟了以下概率。我不确定如何将其翻译成，但请在下面展示我的尝试。我想我以前从未尝试过以这种方式模拟概率。首先代码：SASRSASbhrsimRSAS

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;

我有点惊讶该函数没有同时使用 the 和 parmeters。call rangamalphabeta

以下是上述代码的输出，格式为：SASR

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

这是我试图将其翻译成 .我希望有人可以检查我的代码是否正确并指出任何错误。特别是，我没想到必须使用这一行：in ，但这可能是因为我以前从未这样做过。RRbhrsim <- xgam1 / (xgam1 + xgam2)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)

R SAS 概率伽玛

您能再解释一下SAS代码吗？rangam]（documentation.sas.com/doc/da/vdmmlcdc/8.1/ds2ref/...）的文档说它只需要两个参数，并且（== 形状参数;显然，如果你想要比例（或速率）不等于 1，你需要自己缩放偏差。我不确定如何工作（这里的文档，但也许？调用并将结果分配给？seedacall()call rangam(seed, a, varname)rangam(seed, a)varname

0赞 Ben Bolker 8/12/2023

好的，这个文档（我不知道我在看什么版本，只是在找到一些东西之前四处闲逛）证实了我的猜测。

0赞 Mark Miller 8/12/2023

从我在某处读到的内容来看，该语句的部分用于允许 50 个重复中的每一个重复的种子发生变化。直到昨天，我才看到这个。call

0赞 Mark Miller 8/12/2023

但是我无法重新定位该SAS文档链接。

答：

2赞 Ben Bolker 8/12/2023 #1

根据 SAS 文档，调用 rangam（），使用随机种子生成形状等于的 Gamma 随机偏差，并将结果分配给（并在内部更新种子 - 使用相同参数第二次调用使用 RNG 流中的新值（这个习语对我来说似乎很奇怪，但我知道什么？所以 R 等价物是（一次，在循环的开头）; （R 允许您通过指定其他参数来设置速率或小数位数！= 1;第一个参数是请求的偏差数。call rangam(seed, a, x)aseedxcall rangam(...)seedset.seed(seed)x <- rgamma(1, a)

所以我相信 R 等价物将是

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

值得一提的是，这是生成 Beta（alpha， beta） - 分布式随机数的标准算法（尽管显然与 R 内部使用的任何算法都不相同......

set.seed(111111)
bhrsim <- rbeta(50, shape1 = 2, shape2 = 2)
mean(bhrsim)
## [1] 0.5324306

我有点困惑，因为代码定义了和 .其他地方称为参数，称为参数。但是在代码中，它们在以下几行中都被视为参数：换句话说，看似被定义为代码顶部的参数，但似乎被视为代码本身中的参数。SASbhrbeta = beta parameter for beta distributionbhralpha = alpha parameter for beta distributionαshapeβscaleSASbhralphabhrbetashapecall rangam(seed1,bhralpha,xgam1); call rangam(seed2,bhrbeta,xgam2);bhrbetascaleSASshapeSAS

1赞 Limey 8/12/2023

我同意@MarkMiller。SAS代码在我看来很可疑。我担心它可能没有做它打算做的事情。如果可以的话，我会咨询原始开发人员并检查您可能必须执行的任何文档，以将它应该做什么与实际正在做什么进行比较。我意识到这可能是不可能的。

0赞 Severin Pappadeux 8/15/2023

@MarkMiller只是为了添加一个链接，这里是 Beta（a， b）的生成方式，您可以将其与您的 SAS 源进行比较 en.wikipedia.org/wiki/...

上一个：manyglm（软件包 mvabund）在使用 family = Gamma（link = log）时出现错误

下一个：在 Unity 中更改了颜色

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

simulate probability with gamma distribution in SAS versus R

评论

评论