Gamma分布模拟概率：SAS转R代码正确性验证求助

问题说明 #

需要将一段SAS实现的基于Gamma分布模拟概率bhrsim的代码转换为R语言，同时对SAS中call rangam未同时使用alpha和beta参数存在疑惑，且自行编写的R代码结果与SAS输出差异较大，需检查修正。

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;

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尝试代码 #

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代码存在两个关键错误：

• Gamma分布参数不匹配：SAS的call rangam(seed, alpha, x)生成的是形状参数为alpha，尺度参数固定为1的Gamma分布，而你在R中错误指定了scale = bhrbeta（尺度为2），导致生成的Gamma变量均值为shape*scale=2*2=4，与SAS中Gamma变量均值shape*1=2的分布完全不同。
• xgam2的形状参数错误：SAS中xgam2是用bhrbeta=2作为形状参数生成的Gamma变量，但你在R中写成了shape = bhralpha，逻辑上不符合原代码的设定（若bhrbeta为其他值，错误会直接显现）。

需要说明的是，bhrsim <- xgam1/(xgam1+xgam2)这一步的计算逻辑是正确的：当X_{Gamma(a,1)、Y}Gamma(b,1)时，X/(X+Y)服从Beta(a,b)分布，理论均值为a/(a+b)，这里a=b=2，理论均值为0.5，SAS输出的均值0.52符合模拟结果的随机性。

修正后的R代码 #

# 模拟SAS的两个独立种子逻辑
set.seed(11111111)
xgam1 <- rgamma(50, shape = 2, scale = 1) # 尺度为1，匹配SAS的rangam
set.seed(33333333)
xgam2 <- rgamma(50, shape = 2, scale = 1) # 形状参数为bhrbeta=2，尺度1

bhrsim <- xgam1/(xgam1+xgam2)
mean(bhrsim)
# 输出会接近0.5，和SAS的0.52一致（模拟随机性导致的差异）

由于R与SAS的随机数生成算法不同，无法通过种子完全复现SAS的精确结果，但修正参数后，模拟的分布特征会与SAS一致。

内容的提问来源于stack exchange，提问作者Mark Miller