Zipf dağıtılmış sayıları verimli bir şekilde nasıl oluşturulur?

Question

Şu anda C++ 'daki bazı veri yapılarını kıyaslamaktayım ve Zipf dağıtılmış sayılar üzerinde çalışırken bunları test etmek istiyorum.

Bu sitede verilen jeneratör kullanıyorum: http://www.cse.usf.edu/~christen/tools/toolpage.html

Bir Mersenne Twister jeneratör kullanmak uygulanmasını uyarlanmış.

İyi çalışıyor ama gerçekten yavaş. Benim durumumda, aralık büyük olabilir (yaklaşık bir milyon) ve rastgele sayıların sayısı birkaç milyon olabilir.

Alfa parametresi zamanla değişmez, giderilmiştir.

Tüm sum_prob'yi önceden tahsis etmeye çalıştım. Çok daha hızlı, ama hala büyük aralıkta yavaşlar.

Zipf dağıtılmış sayıları oluşturmanın daha hızlı bir yolu var mı? Daha az hassas bir şey bile memnuniyetle karşılanacaktır.

Teşekkür

Answer 1

Ön hesaplama tek başına çok fazla yardımcı olmaz. Ama açık olduğu gibi, sum_prob birikmiş ve artan düzende. Dolayısıyla, zipf_value değerini bulmak için bir ikili arama kullanırsak, O (n) 'den O'ya (log (n)) bir Zipf dağıtılmış sayı üretme sırasını azaltırdık. Verimlilikte çok fazla gelişme var. `)` Zipf (yapılan her arama için farklı bir değere sahip alpha`

int zipf(double alpha, int n) 
{ 
    static int first = TRUE;  // Static first time flag 
    static double c = 0;   // Normalization constant 
    static double *sum_probs;  // Pre-calculated sum of probabilities 
    double z;      // Uniform random number (0 < z < 1) 
    int zipf_value;    // Computed exponential value to be returned 
    int i;      // Loop counter 
    int low, high, mid;   // Binary-search bounds 

    // Compute normalization constant on first call only 
    if (first == TRUE) 
    { 
    for (i=1; i<=n; i++) 
     c = c + (1.0/pow((double) i, alpha)); 
    c = 1.0/c; 

    sum_probs = malloc((n+1)*sizeof(*sum_probs)); 
    sum_probs[0] = 0; 
    for (i=1; i<=n; i++) { 
     sum_probs[i] = sum_probs[i-1] + c/pow((double) i, alpha); 
    } 
    first = FALSE; 
    } 

    // Pull a uniform random number (0 < z < 1) 
    do 
    { 
    z = rand_val(0); 
    } 
    while ((z == 0) || (z == 1)); 

    // Map z to the value 
    low = 1, high = n, mid; 
    do { 
    mid = floor((low+high)/2); 
    if (sum_probs[mid] >= z && sum_probs[mid-1] < z) { 
     zipf_value = mid; 
     break; 
    } else if (sum_probs[mid] >= z) { 
     high = mid-1; 
    } else { 
     low = mid+1; 
    } 
    } while (low <= high); 

    // Assert that zipf_value is between 1 and N 
    assert((zipf_value >=1) && (zipf_value <= n)); 

    return(zipf_value); 
} 

Answer 2

Kodunuzdaki aşağıdaki satırı zipf() yapılan her arama için n kez yürütülür:

sum_prob = sum_prob + c/pow((double) i, alpha); 

O pow() işlevi çağırmak için gerekli olduğunu üzüntü vericidir, içten, bu işlev çünkü Bir değil iki Taylor dizisi toplar. [pow(x, alpha) == exp(alpha*log(x))]. alpha bir tamsayıysa, elbette, pow()'u basit çarpma ile değiştirerek kodu çok hızlandırabilirsiniz. alpha rasyonel bir sayı ise, o zaman iki Taylor serisinin yerini almak için bir Newton-Raphson iterasyonunu kodlayarak kodu daha az bir dereceye kadar hızlandırabilirsiniz. Son durumun tutulması halinde lütfen öneride bulunun.

Neyse ki, alpha'un değişmediğini belirttiniz. pow((double) i, alpha) tablosunu hazırlayarak kodu fazla hızlandıramazsınız, daha sonra tabloya zipf() numaranızı verebilirsiniz? Bu şekilde, zipf()'un pow()'u aramak zorunda kalmayacaktı. Bunun önemli zaman tasarrufu sağlayacağından şüpheleniyorum.

Yine de daha fazla iyileştirme mümkündür. sumprob() işlevini zipf() dışında bir faktöre çarptınysan ne olur? sumprob()'un kullanımı için daha agresif bir arama tablosu hazırlamıyor musunuz? Belki bu fikirlerden bazıları sizi doğru yönde hareket ettirecektir. Onlarla ne yapamayacağınızı görün.

Güncelleme: Şimdi tekrar gözden geçirilmiş olarak sorduğunuz sorunun bu cevabı kullanamayacağını görüyorum. Şu andan itibaren, sorunuz karmaşık değişken teorisinde bir soruya dönüşebilir. Bu, genellikle bildiğiniz gibi kolay sorular değildir. Yeterince zeki bir matematikçi, ilgili bir yineleme ilişkisini veya normal dağıtımının Box-Muller tekniği gibi bir hile olduğunu keşfetmiş olabilir, ancak eğer varsa, tekniği tanımadım. İyi şanslar. (Muhtemelen sizin için önemli değildir, ancak, N. N. Lebedev'in 1972 tarihli mükemmel kitabı Özel İşlevleri ve Uygulamaları, Rusça'dan İngilizce bir çevirimiçi baskıda mevcuttur.Gerçekten, gerçekten bu sorunu vurmak istiyorum varsa, Lebedev sonraki okuyabilir - ama tabii ki umutsuz ölçüsüdür, değil mi) Bu arada

Answer 3

dayalı hızlı bir yolu var mı? reddetme örneklemesinde, bkz. kod here.

Answer 4

Bulduğum tek C++ 11 Zipf rasgele jeneratörü olasılıkları hesapladı ve std::discrete_distribution kullandı. Bu, küçük aralıklar için iyi çalışır, ancak belleği genişlettiği için Zipf değerlerini çok geniş bir aralıkta (veritabanı testi için, benim durumumda) oluşturmanız gerekiyorsa, kullanışlı değildir. Bu yüzden, C++ 'da aşağıda belirtilen algoritmayı uyguladım.

Bu kodu titizlikle test etmedim ve bazı optimizasyonlar muhtemelen mümkün, ancak yalnızca sabit alan gerektiriyor ve iyi çalışıyor gibi görünüyor.

#include <algorithm> 
#include <cmath> 
#include <random> 

/** Zipf-like random distribution. 
* 
* "Rejection-inversion to generate variates from monotone discrete 
* distributions", Wolfgang Hörmann and Gerhard Derflinger 
* ACM TOMACS 6.3 (1996): 169-184 
*/ 
template<class IntType = unsigned long, class RealType = double> 
class zipf_distribution 
{ 
public: 
    typedef RealType input_type; 
    typedef IntType result_type; 

    static_assert(std::numeric_limits<IntType>::is_integer, ""); 
    static_assert(!std::numeric_limits<RealType>::is_integer, ""); 

    zipf_distribution(const IntType n=std::numeric_limits<IntType>::max(), 
         const RealType q=1.0) 
     : n(n) 
     , q(q) 
     , H_x1(H(1.5) - 1.0) 
     , H_n(H(n + 0.5)) 
     , dist(H_x1, H_n) 
    {} 

    IntType operator()(std::mt19937& rng) 
    { 
     while (true) { 
      const RealType u = dist(rng); 
      const RealType x = H_inv(u); 
      const IntType k = clamp<IntType>(std::round(x), 1, n); 
      if (u >= H(k + 0.5) - h(k)) { 
       return k; 
      } 
     } 
    } 

private: 
    /** Clamp x to [min, max]. */ 
    template<typename T> 
    static constexpr T clamp(const T x, const T min, const T max) 
    { 
     return std::max(min, std::min(max, x)); 
    } 

    /** exp(x) - 1/x */ 
    static double 
    expxm1bx(const double x) 
    { 
     return (std::abs(x) > epsilon) 
      ? std::expm1(x)/x 
      : (1.0 + x/2.0 * (1.0 + x/3.0 * (1.0 + x/4.0))); 
    } 

    /** H(x) = log(x) if q == 1, (x^(1-q) - 1)/(1 - q) otherwise. 
    * H(x) is an integral of h(x). 
    * 
    * Note the numerator is one less than in the paper order to work with all 
    * positive q. 
    */ 
    const RealType H(const RealType x) 
    { 
     const RealType log_x = std::log(x); 
     return expxm1bx((1.0 - q) * log_x) * log_x; 
    } 

    /** log(1 + x)/x */ 
    static RealType 
    log1pxbx(const RealType x) 
    { 
     return (std::abs(x) > epsilon) 
      ? std::log1p(x)/x 
      : 1.0 - x * ((1/2.0) - x * ((1/3.0) - x * (1/4.0))); 
    } 

    /** The inverse function of H(x) */ 
    const RealType H_inv(const RealType x) 
    { 
     const RealType t = std::max(-1.0, x * (1.0 - q)); 
     return std::exp(log1pxbx(t) * x); 
    } 

    /** That hat function h(x) = 1/(x^q) */ 
    const RealType h(const RealType x) 
    { 
     return std::exp(-q * std::log(x)); 
    } 

    static constexpr RealType epsilon = 1e-8; 

    IntType         n;  ///< Number of elements 
    RealType         q;  ///< Exponent 
    RealType         H_x1; ///< H(x_1) 
    RealType         H_n; ///< H(n) 
    std::uniform_real_distribution<RealType> dist; ///< [H(x_1), H(n)] 
};