2011-03-13 19 views
7

मैं अंतराल [0,1] में एक गाऊसी वितरित यादृच्छिक संख्या जनरेटर को लागू करने की कोशिश कर रहा हूँ।गाऊसी यादृच्छिक संख्या जनरेटर

float rand_gauss (void) { 
    float v1,v2,s; 

    do { 
    v1 = 2.0 * ((float) rand()/RAND_MAX) - 1; 
    v2 = 2.0 * ((float) rand()/RAND_MAX) - 1; 

    s = v1*v1 + v2*v2; 
    } while (s >= 1.0); 

    if (s == 0.0) 
    return 0.0; 
    else 
    return (v1*sqrt(-2.0 * log(s)/s)); 
} 

यह TAOCP 3 संस्करण पेज के नुथ के 2 मात्रा में एल्गोरिथ्म के काफी एक सीधे आगे कार्यान्वयन है 122.

समस्या यह है कि rand_gauss() कभी कभी अंतराल [0 बाहर के मान देता है, 1]।

+11

एक गाऊशियन असंबद्ध है। क्या मैं कुछ भूल रहा हूँ? –

+0

एक भिन्नता और मतलब है, मैं 0 के रूप में माध्य और भिन्नता^2 के रूप में 1, मानक सामान्य वितरण है। –

+8

@ एनवीएम: एक सामान्य सामान्य वितरण कुछ संभावनाओं के साथ अनंतता और अनंतता के बीच कोई मूल्य ले सकता है; परिणाम पर कोई सीमा सीमा नहीं है। –

उत्तर

8

नुथ TAOCP के 2 मात्रा के पी 122 पर ध्रुवीय विधि का वर्णन है। यही कारण है कि एल्गोरिथ्म उत्पन्न साथ एक सामान्य वितरण मतलब = 0 और मानक विचलन = 1। लेकिन आप वांछित मानक विचलन और वांछित मतलब जोड़कर गुणा करके इसे समायोजित कर सकते हैं।

आप यह मजेदार the polar method in the C-FAQ का एक और कार्यान्वयन के लिए अपने कोड तुलना करने के लिए मिल सकती है।

4

(s >= 1.0 || s == 0.0) करने के लिए अपने अगर बयान बदलें। बेहतर अभी तक, break का उपयोग करें जैसा कि एक जटिल जोड़ी (यू, वी) लौटने वाले सिमड गाऊस यादृच्छिक संख्या के लिए निम्नलिखित उदाहरण में देखा गया है। यह मेर्सन ट्विस्टर यादृच्छिक संख्या जनरेटरdsfmt() का उपयोग करता है। यदि आप केवल एक सिंगल, असली, यादृच्छिक संख्या चाहते हैं, तो केवल u लौटाएं और अगले पास के लिए v सहेजें।

inline static void randn(double *u, double *v) 
{ 
double s, x, y; // SIMD Marsaglia polar version for complex u and v 
while (1){ 
    x = dsfmt_genrand_close_open(&dsfmt) - 1.; 
    y = dsfmt_genrand_close_open(&dsfmt) - 1.; 
    s = x*x + y*y; 
    if (s < 1) break; 
} 
s = sqrt(-2.0*log(s)/s); 
*u = x*s; *v = y*s; 
return; 
} 

यह एल्गोरिदम आश्चर्यजनक रूप से तेज़ है। दो यादृच्छिक संख्या की गणना के लिए निष्पादन समय (यू, वी) चार अलग अलग गाऊसी यादृच्छिक संख्या जनरेटर के लिए कर रहे हैं:

Times for delivering two Gaussian numbers (u + iv) 
i7-2600K @ 4GHz, gcc -Wall -Ofast -msse2 .. 

gsl_ziggurat      =  20.3 (ns) 
Box-Muller       =  78.8 (ns) 
Box-Muller with fast_sin fast_cos =  28.1 (ns) 
SIMD Marsaglia polar    =  35.0 (ns) 

fast_sin और fast_cos चार्ल्स के गैरेट के बहुपद दिनचर्या एक कारक 2.9 से बॉक्स मुलर गणना में तेजी लाने के कोस() और पाप() के घोंसले वाले बहुपद कार्यान्वयन का उपयोग करना। सिम बॉक्स बॉक्सर और ध्रुवीय एल्गोरिदम निश्चित रूप से प्रतिस्पर्धी हैं। इसके अलावा उन्हें आसानी से समांतर किया जा सकता है। जीसीसी-ओफ़ास्ट-एस का उपयोग करके, असेंबली कोड डंप से पता चलता है कि वर्ग रूट सिम एसएसई 2 है: एसकर्ट -> sqrtsd% xmm0,% xmm0

टिप्पणी: जीसीसी 5 के साथ सटीक समय प्राप्त करने के लिए यह वास्तव में कठिन और निराशाजनक है, लेकिन मुझे लगता है इन ठीक कर रहे हैं: 2016/02/03 के रूप में डीएलडब्लयू

[1] संबंधित लिंक: c malloc array pointer return in cython

[2] एल्गोरिदम की तुलना, लेकिन जरूरी नहीं कि SIMD संस्करणों के लिए: http://www.doc.ic.ac.uk/~wl/papers/07/csur07dt.pdf

[3] चार्ल्स के गैरेट: http://krisgarrett.net/papers/l2approx.pdf

+0

yup, धन्यवाद, बदल गया। 73 – W9DKI

+1

यहां सही होने का कठिन हिस्सा है: gcc -Wall -Ofast -msse2 -frename-registers -malign-double -fno-hard-aliasing -DHAVE_SSE2 = 1 -DDSFMT_MEXP = 19937 -o randn_test dSFMT.c randn_test.c - एलएम - – W9DKI

+1

मुझे लगता है @ डब्ल्यू 9 डीकेआई, आपने गलती से एक डुप्लिक खाता बनाया है। जब आप अगली बार लॉगिन करते हैं, तो उसी तरीके से लॉगिन करें जैसा आपने पहली बार किया है और किसी अन्य तरीके से नहीं। खातों के विलय के लिए भी ध्वज। यह http://stackoverflow.com/help/merging-accounts –

संबंधित मुद्दे