R и байесовская статистика

Question

Следующее - со страницы 224 книги Введение в науку статистики

simulate in R to find the distribution of the posterior probability of p = p1p2.

> p1<-rbeta(10000,16,6);p2<-rbeta(10000,18,4)
> p<-p1*p2

We then give a table of deciles for the posterior distribution function and present a histogram.

> data.frame(quantile(p,d))
quantile.p..d.
0% 0.2825593
10% 0.4660896
20% 0.5094321
30% 0.5422747
40% 0.5712765
50% 0.5968341
60% 0.6209610
70% 0.6477835
80% 0.6776208
90% 0.7187307
100% 0.9234834

d не приведено в книге, и выглядит какследующим образом:

> d<-seq(0, 1, by=0.1)

Тем не менее, он получал совершенно разные результаты для каждого прогона:

> d<-seq(0, 1, by=0.1)
> p1<-rbeta(10000,16,6);p2<-rbeta(10000,18,4)
> p<-p1*p2
> data.frame(quantile(p,d))
     quantile.p..d.
0%        0.2522659
10%       0.4715691
20%       0.5141640
30%       0.5463989
40%       0.5722724
50%       0.5982149
60%       0.6230126
70%       0.6507722
80%       0.6794564
90%       0.7211774
100%      0.9136375
> 
> 
> p1<-rbeta(10000,16,6);p2<-rbeta(10000,18,4)
> p<-p1*p2
> data.frame(quantile(p,d))
     quantile.p..d.
0%        0.2492287
10%       0.4680199
20%       0.5122128
30%       0.5462639
40%       0.5740713
50%       0.5984662
60%       0.6233943
70%       0.6489084
80%       0.6800467
90%       0.7190094
100%      0.8674298
> 

Как это могло иметь такую ​​большую разницу?

Обновлено :

Имеет смысл, основываясь на документе rbeta:

The Beta Random Number Generating Function

Random generation for the beta distribution with parameters shape1 and shape2

rbeta(n, shape1, shape2)

Следующий код в python

from numpy.random import beta
import numpy as np
import pandas as pd
import sys

# pd.set_option('display.max_columns', None)
# pd.set_option('display.expand_frame_repr', False)
# pd.set_option('max_colwidth', -1)
# np.set_printoptions(threshold=sys.maxsize)

p1 = beta(16, 6, 10000)
p2 = beta(18, 4, 10000)

p = p1 * p2
d = np.linspace(0, 1, 11)

df = pd.DataFrame({
    'd': ["{0:.0f}%".format(val * 100) for val in d],
    'quantile': np.quantile(p, d)
})

print(df)

Вывод:

       d  quantile
0   0%    0.247273
1   10%   0.468004
2   20%   0.513628
3   30%   0.544105
4   40%   0.571387
5   50%   0.596308
6   60%   0.620983
7   70%   0.646542
8   80%   0.677496
9   90%   0.716053
10  100%  0.912590