是否有一個庫根據JavaScript的beta版分佈生成隨機數？

Question

我需要從Beta probability distribution內的JavaScript生成隨機數。我谷歌搜索，但無法找到任何似乎支持這個庫。

任何人都可以建議在哪裏可以找到一個圖書館或代碼片段，這將做到這一點？

Answer 1

我的翻譯。這幾乎是一個字，所以它可能不是最慣用的JavaScript。

// javascript shim for Python's built-in 'sum' 
function sum(nums) { 
    var accumulator = 0; 
    for (var i = 0, l = nums.length; i < l; i++) 
    accumulator += nums[i]; 
    return accumulator; 
} 

// In case you were wondering, the nice functional version is slower. 
// function sum_slow(nums) { 
// return nums.reduce(function(a, b) { return a + b; }, 0); 
// } 
// var tenmil = _.range(1e7); sum(tenmil); sum_slow(tenmil); 

// like betavariate, but more like R's name 
function rbeta(alpha, beta) { 
    var alpha_gamma = rgamma(alpha, 1); 
    return alpha_gamma/(alpha_gamma + rgamma(beta, 1)); 
} 

// From Python source, so I guess it's PSF Licensed 
var SG_MAGICCONST = 1 + Math.log(4.5); 
var LOG4 = Math.log(4.0); 

function rgamma(alpha, beta) { 
    // does not check that alpha > 0 && beta > 0 
    if (alpha > 1) { 
    // Uses R.C.H. Cheng, "The generation of Gamma variables with non-integral 
    // shape parameters", Applied Statistics, (1977), 26, No. 1, p71-74 
    var ainv = Math.sqrt(2.0 * alpha - 1.0); 
    var bbb = alpha - LOG4; 
    var ccc = alpha + ainv; 

    while (true) { 
     var u1 = Math.random(); 
     if (!((1e-7 < u1) && (u1 < 0.9999999))) { 
     continue; 
     } 
     var u2 = 1.0 - Math.random(); 
     v = Math.log(u1/(1.0-u1))/ainv; 
     x = alpha*Math.exp(v); 
     var z = u1*u1*u2; 
     var r = bbb+ccc*v-x; 
     if (r + SG_MAGICCONST - 4.5*z >= 0.0 || r >= Math.log(z)) { 
     return x * beta; 
     } 
    } 
    } 
    else if (alpha == 1.0) { 
    var u = Math.random(); 
    while (u <= 1e-7) { 
     u = Math.random(); 
    } 
    return -Math.log(u) * beta; 
    } 
    else { // 0 < alpha < 1 
    // Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 
    while (true) { 
     var u3 = Math.random(); 
     var b = (Math.E + alpha)/Math.E; 
     var p = b*u3; 
     if (p <= 1.0) { 
     x = Math.pow(p, (1.0/alpha)); 
     } 
     else { 
     x = -Math.log((b-p)/alpha); 
     } 
     var u4 = Math.random(); 
     if (p > 1.0) { 
     if (u4 <= Math.pow(x, (alpha - 1.0))) { 
      break; 
     } 
     } 
     else if (u4 <= Math.exp(-x)) { 
     break; 
     } 
    } 
    return x * beta; 
    } 
} 

部分可測試的裝置，其被容易地計算出：

function testbeta(a, b, N) { 
    var sample_mean = sum(_.range(N).map(function() { return rbeta(a, b); }))/N; 
    var analytic_mean = a/(a + b); 
    console.log(sample_mean, "~", analytic_mean); 
} 
testbeta(5, 1, 100000); 

Answer 2

您可以在此Python代碼轉換爲JS：

SG_MAGICCONST = 1.0 + _log(4.5) 
LOG4 = log(4.0) 

def gamma(z, sqrt2pi=(2.0*pi)**0.5): 
    # Reflection to right half of complex plane 
    if z < 0.5: 
     return pi/sin(pi*z)/gamma(1.0-z) 
    # Lanczos approximation with g=7 
    az = z + (7.0 - 0.5) 
    return az ** (z-0.5)/exp(az) * sqrt2pi * fsum([ 
    0.9999999999995183, 
    676.5203681218835/z, 
    -1259.139216722289/(z+1.0), 
    771.3234287757674/(z+2.0), 
    -176.6150291498386/(z+3.0), 
    12.50734324009056/(z+4.0), 
    -0.1385710331296526/(z+5.0), 
    0.9934937113930748e-05/(z+6.0), 
    0.1659470187408462e-06/(z+7.0), 
    ]) 



def gammavariate(self, alpha, beta): 
    """Gamma distribution. Not the gamma function! 

    Conditions on the parameters are alpha > 0 and beta > 0. 

    The probability distribution function is: 

     x ** (alpha - 1) * math.exp(-x/beta) 
    pdf(x) = -------------------------------------- 
      math.gamma(alpha) * beta ** alpha 

    """ 

    # alpha > 0, beta > 0, mean is alpha*beta, variance is alpha*beta**2 

    # Warning: a few older sources define the gamma distribution in terms 
    # of alpha > -1.0 
    if alpha <= 0.0 or beta <= 0.0: 
    raise ValueError, 'gammavariate: alpha and beta must be > 0.0' 

    random = self.random 
    if alpha > 1.0: 

    # Uses R.C.H. Cheng, "The generation of Gamma 
    # variables with non-integral shape parameters", 
    # Applied Statistics, (1977), 26, No. 1, p71-74 

    ainv = _sqrt(2.0 * alpha - 1.0) 
    bbb = alpha - LOG4 
    ccc = alpha + ainv 

    while 1: 
     u1 = random() 
     if not 1e-7 < u1 < .9999999: 
     continue 
     u2 = 1.0 - random() 
     v = _log(u1/(1.0-u1))/ainv 
     x = alpha*_exp(v) 
     z = u1*u1*u2 
     r = bbb+ccc*v-x 
     if r + SG_MAGICCONST - 4.5*z >= 0.0 or r >= _log(z): 
     return x * beta 

    elif alpha == 1.0: 
    # expovariate(1) 
    u = random() 
    while u <= 1e-7: 
     u = random() 
    return -_log(u) * beta 

    else: # alpha is between 0 and 1 (exclusive) 

    # Uses ALGORITHM GS of Statistical Computing - Kennedy & Gentle 

    while 1: 
     u = random() 
     b = (_e + alpha)/_e 
     p = b*u 
     if p <= 1.0: 
     x = p ** (1.0/alpha) 
     else: 
     x = -_log((b-p)/alpha) 
     u1 = random() 
     if p > 1.0: 
     if u1 <= x ** (alpha - 1.0): 
      break 
     elif u1 <= _exp(-x): 
     break 
    return x * beta 



def betavariate(alpha, beta): 
    if y == 0: 
    return 0.0 
    else: 
    return y/(y + gammavariate(beta, 1.0)) 

這是直接從Python源代碼（略有修改），但它應該很容易轉換。

Answer 3

http://mbostock.github.com/d3/ex/stack.html

你必須生成使用γ分佈堆積條形圖層的一個例子。安靜接近你正在尋找我猜。

Answer 4

的jStat library具有功能從β分佈來樣，以及許多other distributions。

var random_num = jStat.beta.sample(alpha, beta); 

Answer 5

查看stdlib，其中包括許多分佈的可用的PRNG，包括beta分佈。例如，stdlib開發環境中，

var beta = require('@stdlib/math/base/random/beta'); 

var r = beta(2.0, 5.0); 
// returns <number> 

否則，看到這是在Apache許可下發布源代碼。