BigQueryでベータ分布に従う乱数を生成し、本当にできているのか確認する

CREATE OR REPLACE FUNCTION math.generate_beta(a FLOAT64, b FLOAT64, n INT64) AS (
# ベータ分布に従う乱数を生成する
# https://qiita.com/piyo7/items/36b19c6ea68baa2a69cd
ARRAY(
  WITH assert AS (
    SELECT
      IF(a > 0 AND b > 0, TRUE, ERROR('パラメータのa,bはいずれも正であることが必要です')) AS c1,
      IF(n > 0, TRUE, ERROR('nは1以上です')) AS c2
  )
  , ab AS (
    SELECT
      a+b AS alpha,
      IF(
        math.min([a, b]) <= 1,
        math.max([1/a, 1/b]),
        SQRT((a+b-2)/(2*a*b-(a+b)))
      ) AS beta,
    FROM assert
  )
  , abg AS (
    SELECT alpha, beta,
      a + (1 / beta) AS gamma
    FROM ab
  )
  , u AS (
  SELECT RAND() AS U1, RAND() AS U2
  # 棄却される可能性があるので多めに作る
  FROM UNNEST(GENERATE_ARRAY(0, CEIL(n * 1.5)))
  )
  , vw AS (
    SELECT
      U1, U2,
       beta * LN(U1/(1-U1)) AS V,
       a * EXP(beta * LN(U1/(1-U1))) AS W
    FROM u, abg
  )
  , sampling AS (
    SELECT
      ROW_NUMBER() OVER () AS rank,
      alpha * LN(a / (b+W)) + gamma * V - 1.3862944 < LN(U1*U1*U2) AS accept,
      W / (b+W) AS X,
    FROM vw, abg
  )
  SELECT X
  FROM sampling
  WHERE accept AND rank <= n
));

CREATE OR REPLACE FUNCTION math.histogram(arr ANY TYPE, bin INT64) AS (
# 配列を値によりビンにわけ個数をカウントする
ARRAY(
  WITH min_max AS (
    SELECT
      math.min(arr) AS x_min,
      (math.max(arr) - math.min(arr)) / bin AS d,
  )
  , bin_range AS (
    SELECT
      i,
      x_min + i * d AS range_left,      -- それぞれのビンの範囲の左端
      x_min + (i+1) * d AS range_right, -- それぞれのビンの範囲の右端
    FROM UNNEST(GENERATE_ARRAY(0, bin-1)) AS i, min_max
  )
  , hist AS (
    SELECT
      i,
      COUNTIF(i IS NOT NULL) AS n,
    FROM UNNEST(arr) AS x, bin_range
    WHERE TRUE
      AND x >= range_left
      AND IF(i=bin-1, x <= range_right, x < range_right)
    GROUP BY i
  )
  SELECT AS STRUCT
    range_left AS l,
    range_right AS r,
    IFNULL(n, 0) AS n,
  FROM bin_range LEFT JOIN hist USING (i)
));

SELECT (l + r) / 2 AS x, n
FROM UNNEST(math.histogram(math.generate_beta(2, 5, 10000), 100))