ग में पैदा प्वाइजन चर ++

Question

मैं एक प्वाइजन यादृच्छिक चर

typedef long unsigned int luint; 
luint poisson(luint lambda) { 
    double L = exp(-double(lambda)); 
    luint k = 0; 
    double p = 1; 
    do { 
     k++; 
     p *= mrand.rand(); 
    } while(p > L); 
    return (k-1); 
} 

जहां mrand MersenneTwister यादृच्छिक संख्या जनरेटर है उत्पन्न करने के लिए इस समारोह को लागू किया। मुझे लगता है कि, जैसा कि मैंने लैम्ब्डा को बढ़ाया है, अपेक्षित वितरण गलत होगा, जिसका मतलब है कि लगभग 750 पर संतृप्त होता है। क्या यह संख्यात्मक अनुमानों के कारण है या क्या मैंने कोई गलती की है?

Answer 1

another question I asked earlier से, यह आप कर सकते थे भी अनुमानित poisson(750)poisson(375) + poisson(375) के रूप में लगता है।

Answer 2

एक्सपी (-750) एक बहुत छोटी संख्या है, जो कि सबसे छोटी संभव डबल के करीब है, इसलिए आपकी समस्या संख्यात्मक है। किसी भी मामले में, आपकी जटिलता लैम्ब्डा में रैखिक होगी, इसलिए एल्गोरिदम उच्च लैम्ब्डा के लिए बहुत प्रभावी नहीं है। जब तक आपके पास इसे स्वयं कोड करने का कोई बड़ा कारण न हो, तब तक मौजूदा लाइब्रेरी कार्यान्वयन का उपयोग करने से शायद यह समझ में आता है, क्योंकि इन संख्यात्मक एल्गोरिदम आपके द्वारा सामना किए जा रहे सटीक मुद्दों के लिए सटीक रूप से स्पर्श किए जाते हैं।

Answer 3

चूंकि आप केवल L अभिव्यक्ति (p>L) में उपयोग करते हैं, तो आप अनिवार्य रूप से (log(p) > -lambda) के लिए परीक्षण कर रहे हैं। यह एक बहुत उपयोगी परिवर्तन नहीं है। निश्चित रूप से, आपको अब एक्स (-750) की आवश्यकता नहीं है, लेकिन आप इसके बजाय p ओवरफ़्लो कर देंगे।

अब, p सिर्फ Π (mrand.rand()) है, और लॉग (पी) लॉग है (Π (mrand.rand())) Σ (log (mrand.rand()) है। यह आपको देता है आवश्यक परिवर्तन:

double logp = 0; 
do { 
    k++; 
    logp += log(mrand.rand()); 
} while(logp > -lambda); 

double प्रतिपादक के केवल 11 बिट है, लेकिन एक 52 बिट अपूर्णांश इसलिए इस संख्यात्मक स्थिरता में भारी वृद्धि है कीमत चुकानी पड़ी है कि आप हर यात्रा पर एक log की जरूरत है, के बजाय है।। एक exp सामने।

Answer 4

इन स्थितियों में, आपको यादृच्छिक संख्या जनरेटर को एक से अधिक बार आमंत्रित करने की आवश्यकता नहीं है। सभी यो

double c[k] = // the probability that X <= k (k = 0,...) 

फिर एक यादृच्छिक संख्या 0 <= r < 1 उत्पन्न, और पहली पूर्णांक X ऐसी है कि c[X] > r ले: यू की जरूरत संचयी संभावनाओं की एक टेबल है। आप इस X को बाइनरी खोज के साथ पा सकते हैं।

इस तालिका बनाने के लिए, हम की जरूरत अलग-अलग संभावनाओं

p[k] = lambda^k/(k! e^lambda) // // the probability that X = k 

तो lambda बड़ी है, इस बेतहाशा गलत हो जाता है, जैसा कि आप मिल गया है। लेकिन हम यहां एक चाल का उपयोग कर सकते हैं: k = floor[lambda] के साथ सबसे बड़े मूल्य पर (या निकट) शुरू करें, और उस पल के लिए नाटक करें कि p[k]1 के बराबर है। फिर आवर्ती संबंध

p[i+1] = (p[i]*lambda)/(i+1) 

और i < k के लिए उपयोग कर

p[i-1] = (p[i]*i)/lambda 

यह सुनिश्चित करता है कि सबसे बड़ी संभावनाओं सबसे बड़ी संभव सटीक का उपयोग कर i > k के लिए p[i] गणना।

अब बस c[i], c[i+1] = c[i] + p[i+1] का उपयोग कर की गणना बिंदु तक जहां c[i+1]c[i] के समान है। फिर आप इस सीमित मूल्य c[i] द्वारा विभाजित करके सरणी को सामान्यीकृत कर सकते हैं; या आप सरणी को छोड़ सकते हैं, और यादृच्छिक संख्या 0 <= r < c[i] का उपयोग कर सकते हैं।

देखें: http://en.wikipedia.org/wiki/Inverse_transform_sampling

Answer 5

यदि आप "मौजूदा लाइब्रेरी" मार्ग पर जाते हैं, तो आपका कंपाइलर पहले से ही C++ 11 std :: यादृच्छिक पैकेज का समर्थन कर सकता है। यहाँ कैसे आप इसका इस्तेमाल होता है:

#include <random> 
#include <ctime> 
#include <iostream> 

std::mt19937 mrand(std::time(0)); // seed however you want 

typedef long unsigned int luint; 

luint poisson(luint lambda) 
{ 
    std::poisson_distribution<luint> d(lambda); 
    return d(mrand); 
} 

int main() 
{ 
    std::cout << poisson(750) << '\n'; 
    std::poisson_distribution<luint> d(750); 
    std::cout << d(mrand) << '\n'; 
    std::cout << d(mrand) << '\n'; 
} 

मैं इसे का उपयोग किया है इसके बाद के संस्करण के दो तरीके:

1. मैंने कोशिश की अपने मौजूदा इंटरफेस की नकल करने की।

2. यदि आप एक माध्य के साथ std :: poisson_distribution बनाते हैं, तो उसी वितरण के लिए उस वितरण का उपयोग करने के लिए और अधिक कुशल है (जैसा कि मुख्य() में किया गया है)।

यहाँ मेरे लिए नमूना उत्पादन होता है:

751 
730 
779