Jak efektywnie generować numery rozproszone Zipf?

Question

Aktualnie testuję niektóre struktury danych w C++ i chcę je przetestować podczas pracy nad liczbami rozproszonymi w Zipf.

Używam generator dostępne na tej stronie: http://www.cse.usf.edu/~christen/tools/toolpage.html

I dostosowany wdrażanie do korzystania generator Mersenne Twister.

Działa dobrze, ale jest naprawdę powolny. W moim przypadku zakres może być duży (około miliona), a liczba generowanych liczb losowych może wynosić kilka milionów.

Parametr alfa nie zmienia się z czasem, jest ustalony.

Próbowałem uprzedzić wszystkie sum_prob. Jest o wiele szybszy, ale wciąż zwalnia na dużym obszarze.

Czy istnieje szybszy sposób generowania rozproszonych numerów Zipf? Nawet coś mniej precyzyjnego będzie mile widziane.

Dzięki

Answer 1

Same obliczenia wstępne niewiele pomagają. Ale jak to jest oczywiste, sum_prob ma charakter akumulacyjny i ma rosnącą kolejność. Jeśli więc użyjemy wyszukiwania binarnego do znalezienia wartości zipf, zmniejszymy kolejność generowania liczby rozproszonej Zipf z O (n) do O (log (n)). Co oznacza tyle poprawy wydajności.

to jest tutaj, po prostu zastąpić funkcję zipf() w genzipf.c z następujący:

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

następującą linię w kodzie jest realizowane n razy dla każdego wezwania do zipf():

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

Szkoda, że ​​konieczne jest, aby wywołać funkcję pow() ponieważ wewnętrznie tej funkcji sumuje nie jedną, ale dwie serie Taylora [z uwzględnieniem pow(x, alpha) == exp(alpha*log(x))]. Jeśli alpha jest liczbą całkowitą, oczywiście, możesz znacznie przyspieszyć kod, zastępując pow() prostym mnożeniem. Jeśli alpha jest liczbą wymierną, to możesz w mniejszym stopniu przyspieszyć kod, kodując iterację Newtona-Raphsona, aby zastąpić dwie serie Taylora. Jeśli ostatni warunek jest spełniony, proszę doradzić.

Na szczęście zaznaczyłeś, że alpha się nie zmienia. Czy nie możesz dużo przyspieszyć kodu, przygotowując tabelę o numerze pow((double) i, alpha), a następnie pozwalając zipf() wyszukać numery w tabeli? W ten sposób zipf() nie musiałby w ogóle dzwonić pod numer pow(). Podejrzewam, że zaoszczędziłoby to dużo czasu.

Możliwe są dalsze ulepszenia. Co się stanie, jeśli uwzględnisz funkcję sumprob() z zipf()? Czy nie możesz przygotować jeszcze bardziej agresywnego stołu przeglądowego do korzystania z sumprob()?

Może niektóre z tych pomysłów popchną Cię w dobrym kierunku. Zobacz, czego nie możesz z nimi zrobić.

Aktualizacja: Widzę, że twoje pytanie, które teraz zostało zmienione, może nie być w stanie wykorzystać tej odpowiedzi. Od tego momentu twoje pytanie może rozwiązać pytanie w skomplikowanej teorii zmiennych. Takie pytania często nie są łatwe, jak wiesz. Możliwe, że wystarczająco sprytny matematyk odkrył stosowną relację lub jakiś trik, taki jak technika Box-Muller dystrybucji, ale jeśli tak, to nie jestem zaznajomiony z tą techniką. Powodzenia. (To prawdopodobnie nie ma znaczenia dla ciebie, ale na wszelki wypadek późna książka N. N. Lebiediewa z 1972 r. Funkcje specjalne i ich zastosowania jest dostępna w rosyjskim tłumaczeniu na język angielski w niedrogiej wersji z miękką okładką.Jeśli naprawdę jesteś zainteresowany, naprawdęchciałeś rozwiązać ten problem, możesz przeczytać Lebiediewa dalej - ale, oczywiście, jest to desperacka miara, prawda?)

Answer 3

W międzyczasie istnieje szybsza metoda oparta na w przypadku próbkowania inwersji odrzucenia, patrz kod here.

Answer 4

Jedyny generator losowy C++ 11 Zipf, jaki mogłem znaleźć, obliczał prawdopodobieństwa jawnie i używał std::discrete_distribution. Działa to dobrze dla małych zakresów, ale nie jest użyteczne, jeśli potrzebujesz generować wartości Zipf z bardzo szerokim zakresem (dla testowania bazy danych, w moim przypadku), ponieważ spowoduje to wyczerpanie pamięci. Więc zaimplementowałem poniższy algorytm w C++.

Nie rygorystycznie testowałem tego kodu, a niektóre optymalizacje są prawdopodobnie możliwe, ale wymaga to tylko stałej przestrzeni i wygląda na to, że działa dobrze.

#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)] 
};