पायथन में बायेसियन विश्लेषण

यह पोस्ट बायेसियन विधियों पर मेरी पहली पोस्ट की तार्किक निरंतरता है, जो यहां पाई जा सकती है ।
मैं व्यवहार में विश्लेषण करने के तरीके के बारे में विस्तार से बात करना चाहूंगा।

जैसा कि मैंने उल्लेख किया है, बेयसियन विश्लेषण के लिए उपयोग किया जाने वाला सबसे लोकप्रिय उपकरण जेएजीएस और / या बीयूजीएस पैकेज के साथ आर भाषा है। औसत उपयोगकर्ता के लिए, पैकेजों में अंतर JAGS की क्रॉस-प्लेटफ़ॉर्म प्रकृति है (मैंने अपने अनुभव से देखा है कि उबंटू और BUGS के बीच एक संघर्ष है), और यह भी कि आप JAGS में अपने स्वयं के कार्यों और वितरणों को बना सकते हैं (हालांकि, औसत उपयोगकर्ता को इसकी आवश्यकता है) अक्सर), अगर सब पर)। वैसे, R के लिए एक उत्कृष्ट और सुविधाजनक IDE है, उदाहरण के लिए, RStudio ।
लेकिन इस पोस्ट में मैं pymc मॉड्यूल के साथ R - Python के विकल्प के बारे में लिखूंगा ।
मैं स्पाइडर को पायथन के लिए सुविधाजनक आईडीई के रूप में सुझा सकता हूं।
मैं पायथन को प्राथमिकता देता हूं, क्योंकि, सबसे पहले, मुझे आर के रूप में इस तरह की एक भी व्यापक भाषा नहीं सीखने का कोई कारण नहीं दिखता है, बस आंकड़ों से संबंधित किसी प्रकार की समस्या को गिनने के लिए। दूसरी बात, पायथन अपने मॉड्यूल के साथ, मेरे दृष्टिकोण से, सादगी, समझदारी और कार्यक्षमता में किसी भी चीज़ में R से नीच नहीं है।

मैं डेटा की रैखिक निर्भरता के गुणांक खोजने की सरल समस्या को हल करने का प्रस्ताव करता हूं। ज्ञान के विभिन्न क्षेत्रों में पैरामीटर अनुकूलन की एक समान समस्या काफी आम है, इसलिए मैं कहूंगा कि यह बहुत खुलासा है। अंत में, हम बायेसियन विश्लेषण के परिणाम और सबसे कम वर्ग विधि की तुलना कर सकते हैं।

सॉफ्टवेयर स्थापना

सबसे पहले, हमें पायथन को स्थापित करने की आवश्यकता है (उन लोगों के लिए जिन्होंने इसे अभी तक नहीं किया है)। मैंने पायथन 3.3 का उपयोग नहीं किया, लेकिन 2.7 के साथ सब कुछ निश्चित रूप से काम करता है।
फिर आपको पायथन के लिए सभी आवश्यक मॉड्यूल स्थापित करने की आवश्यकता है: सुन्न , Matplotlib ।
यदि वांछित है, तो आप अतिरिक्त मॉड्यूल भी स्थापित कर सकते हैं: स्कैपी , पाईटेबल्स , पीडॉट , आईपीथॉन और नाक ।
इन सभी प्रतिष्ठानों (पायथन को छोड़कर) को सेटपूल के माध्यम से करना आसान है।
अब आप स्वयं pymc स्थापित कर सकते हैं (इसे सेटप्टूल के माध्यम से भी स्थापित किया जा सकता है)।

एक विस्तृत pymc गाइड यहां पाया जा सकता है ।

कार्य विवरण

हमारे पास एक काल्पनिक प्रयोग के दौरान प्राप्त डेटा है, जो x के कुछ मान पर रैखिक रूप से निर्भर करता है। डेटा शोर के साथ आता है जिसका फैलाव अज्ञात है। रैखिक निर्भरता के गुणांक को खोजने के लिए आवश्यक है।

निर्णय

सबसे पहले, हम उन मॉड्यूल को आयात करते हैं जिन्हें हमने स्थापित किया था और जिनकी हमें आवश्यकता है:

import numpy import pymc 

फिर हमें अपने काल्पनिक रैखिक डेटा प्राप्त करने की आवश्यकता है।
ऐसा करने के लिए, हम निर्धारित करते हैं कि हम कितने अंक चाहते हैं (इस मामले में 20), वे समान रूप से अंतराल [0, 10] पर वितरित किए जाते हैं, और हम रैखिक निर्भरता के वास्तविक गुणांक निर्धारित करते हैं। अगला, हम गॉसियन डेटा पर शोर थोपते हैं:

 #Generate data with noise number_points = 20 true_coefficients = [10.4, 5.5] x = numpy.linspace(0, 10, number_points) noise = numpy.random.normal(size = number_points) data = true_coefficients[0]*x + true_coefficients[1] + noise 

इसलिए, हमारे पास आंकड़े हैं, और अब हमें यह सोचने की जरूरत है कि हम विश्लेषण कैसे करना चाहते हैं।
सबसे पहले, हम जानते हैं (या मान लें) कि हमारा शोर गाऊसी है, जिसका अर्थ है कि संभावना समारोह गाऊसी होगा। उसके दो पैरामीटर हैं: औसत और विचरण। चूंकि औसत शोर शून्य है, संभावना फ़ंक्शन के लिए औसत मॉडल के मूल्य से निर्धारित किया जाएगा (और हमारे पास एक रैखिक मॉडल है, इसलिए दो पैरामीटर हैं)। जबकि विचरण हमारे लिए अज्ञात है, इसलिए, यह एक और पैरामीटर होगा।
नतीजतन, हमारे पास तीन पैरामीटर और एक गाऊसी संभावना समारोह है।
हम मापदंडों के मूल्यों के बारे में कुछ भी नहीं जानते हैं, इसलिए एक प्राथमिकता हम मानते हैं कि उन्हें समान रूप से मनमानी सीमाओं के साथ वितरित किया जाता है (इन सीमाओं को मनमाने ढंग से दूर ले जाया जा सकता है)।
प्राथमिकता वितरण को निर्दिष्ट करते समय, हमें उस लेबल को इंगित करना चाहिए जिसके द्वारा हम मापदंडों (पहले तर्क) के पीछे के मूल्यों के बारे में जानेंगे, और वितरण की सीमाओं (दूसरी और तीसरी बहस) का भी संकेत देंगे। उपरोक्त सभी तर्कों की आवश्यकता है (अभी भी अतिरिक्त तर्क हैं, जिनमें से एक विवरण प्रलेखन में पाया जा सकता है)।

 sigma = pymc.Uniform('sigma', 0., 100.) a = pymc.Uniform('a', 0., 20.) b = pymc.Uniform('b', 0., 20.) 

अब हमें अपना मॉडल सेट करना होगा। Pymc में दो सबसे अधिक उपयोग किए जाने वाले वर्ग हैं: नियतात्मक और स्टोचैस्टिक। यदि, दिए गए इनपुट डेटा के लिए, यह विशिष्ट रूप से मूल्य (ओं) को निर्धारित करने के लिए संभव है कि मॉडल वापस लौटे, तो यह एक निर्धारक मॉडल है। हमारे मामले में, किसी भी बिंदु के लिए रैखिक निर्भरता के गुणांक के लिए, हम विशिष्ट रूप से क्रमशः परिणाम निर्धारित कर सकते हैं, यह एक निर्धारक मॉडल है:

 @pymc.deterministic(plot=False) def linear_fit(a=a, b=b, x=x): return a*x + b 

और अंत में, हम संभावना फ़ंक्शन सेट करते हैं, जिसमें औसत मॉडल का मूल्य है, सिग्मा एक दिए गए प्राथमिक वितरण के साथ पैरामीटर है, और डेटा हमारा प्रयोगात्मक डेटा है:

 y = pymc.Normal('y', mu=linear_fit, tau=1.0/sigma**2, value=data, observed=True) 

तो, संपूर्ण मॉडल-थ्रेड फ़ाइल इस तरह दिखती है:

 import numpy import pymc #Generate data with noise number_points = 20 true_coefficients = [10.4, 5.5] x = numpy.linspace(0, 10, number_points) noise = numpy.random.normal(size = number_points) data = true_coefficients[0]*x + true_coefficients[1] + noise #PRIORs: #as sigma is unknown then we define it as a parameter: sigma = pymc.Uniform('sigma', 0., 100.) #fitting the line y = a*x+b, hence the coefficient are parameters: a = pymc.Uniform('a', 0., 20.) b = pymc.Uniform('b', 0., 20.) #define the model: if a, b and x are given the return value is determined, hence the model is deterministic: @pymc.deterministic(plot=False) def linear_fit(a=a, b=b, x=x): return a*x + b #LIKELIHOOD #normal likelihood with observed data (with noise), model value and sigma y = pymc.Normal('y', mu=linear_fit, tau=1.0/sigma**2, value=data, observed=True) 

अब मैं एक छोटा सा सैद्धांतिक विषयांतर बनाने का सुझाव देता हूं और चर्चा करता हूं कि pymc आखिर करता क्या है।

गणित की दृष्टि से, हमें निम्नलिखित समस्या को हल करने की आवश्यकता है:
पी (ए, बी, सिग्मा | डेटा) = पी (डेटा | ए, बी, सिग्मा) * पी (ए, बी, सिग्मा) / पी (डेटा)

क्योंकि ए, बी और सिग्मा स्वतंत्र हैं, फिर हम निम्नलिखित रूप में समीकरण को फिर से लिख सकते हैं:
p (a, b, sigma | Data) = p (Data | a, b, sigma) * p (a) * p (b) * p (sigma) / p (Data)

कागज पर, कार्य बहुत सरल दिखता है, लेकिन जब हम इसे संख्यात्मक रूप से हल करते हैं (हमें इस तथ्य के बारे में भी सोचना चाहिए कि हम इस वर्ग की किसी भी समस्या को संख्यात्मक रूप से हल करना चाहते हैं, और न केवल कुछ प्रकार की संभावनाओं के साथ), तो कठिनाइयां पैदा होती हैं।

पी (डेटा), जैसा कि मेरी पिछली पोस्ट में चर्चा की गई थी, एक निरंतर।
p (डेटा | a, b, sigma) निश्चित रूप से हमें दिया जाता है (जो कि a, b और sigma के साथ ज्ञात है, हम हमारे उपलब्ध डेटा की संभावनाओं की विशिष्ट गणना कर सकते हैं)
लेकिन पी (ए), पी (बी) और पी (सिग्मा) के बजाय, वास्तव में, हम केवल हमारे द्वारा निर्धारित कानून के अनुसार वितरित किए गए छद्म आयामी चर के जनरेटर हैं।
कैसे इस सब से एक पीछे वितरण पाने के लिए? यह सही है, जनरेट (चयन करें) a, b और सिग्मा, और फिर p (डेटा | a, b, sigma) की गणना करें। नतीजतन, हमें मूल्यों की एक श्रृंखला मिलती है, जो कि पश्च वितरण से एक नमूना है। लेकिन इससे यह सवाल उठता है कि हम इस चयन को सही तरीके से कैसे कर सकते हैं। यदि हमारे बाद के वितरण में कई मोड ("हिल्स") हैं, तो हम सभी मोड को कवर करने वाला एक नमूना कैसे उत्पन्न कर सकते हैं। यही है, समस्या यह है कि प्रभावी ढंग से एक चयन कैसे किया जाए जो कम से कम पुनरावृत्तियों में संपूर्ण वितरण को "कवर" करेगा। इसके लिए कई एल्गोरिदम हैं, जिनमें से सबसे अधिक उपयोग MCMC (मार्कोव चेन मोंटे कार्लो) है। मार्कोव श्रृंखला यादृच्छिक घटनाओं का एक क्रम है जिसमें प्रत्येक तत्व पिछले एक पर निर्भर करता है, लेकिन पिछले एक पर निर्भर नहीं करता है। मैं स्वयं एल्गोरिथ्म का वर्णन नहीं करूंगा (यह एक अलग पोस्ट का विषय हो सकता है), लेकिन मैं केवल यह नोट करूंगा कि pymc इस एल्गोरिथ्म को लागू करता है और परिणामस्वरूप मार्कोव श्रृंखला देता है, जो पश्च वितरण से एक नमूना है। सामान्यतया, यदि हम नहीं चाहते कि श्रृंखला मार्कोव हो, तो हमें इसे "पतला" करने की आवश्यकता है, अर्थात्। उदाहरण के लिए, हर दूसरे तत्व को लें।
तो, हम दूसरी फ़ाइल बनाते हैं, चलिए इसे रन_मॉडल एरोकोम कहते हैं, जिसमें हम मार्कोव श्रृंखला उत्पन्न करेंगे। Model.py और run_model.py फाइलें एक ही फ़ोल्डर में होनी चाहिए, अन्यथा आपको run_model.py फ़ाइल में कोड जोड़ना होगा:

 from sys import path path.append("/////model.py/") 

सबसे पहले, हम कुछ मॉड्यूल आयात करते हैं जो हमारे लिए उपयोगी हैं:

 from numpy import polyfit from matplotlib.pyplot import figure, plot, show, legend import pymc import model 

पॉलीफ़िट कम से कम वर्गों की विधि को लागू करता है - इसके साथ हम बायेसियन विश्लेषण के परिणामों की तुलना करते हैं।
अंतिम ग्राफ बनाने के लिए आकृति, कथानक, शो, किंवदंती की आवश्यकता होती है।
मॉडल वास्तव में, हमारा मॉडल है।

फिर हम MCMC ऑब्जेक्ट बनाते हैं और चयन चलाते हैं:

 D = pymc.MCMC(model, db = 'pickle') D.sample(iter = 10000, burn = 1000) 

D. नमूना दो तर्क लेता है (वास्तव में आप अधिक निर्दिष्ट कर सकते हैं) - पुनरावृत्तियों और बर्न-इन की संख्या (इसे "वार्म-अप अवधि" कहते हैं)। वार्म-अप अवधि पहली पुनरावृत्तियों की संख्या है जो कि क्लिप की जाती हैं। तथ्य यह है कि एमसीएमसी पहले शुरुआती बिंदु पर निर्भर करता है (यह एक ऐसी संपत्ति है), इसलिए हमें निर्भरता के इस अवधि को काटने की आवश्यकता है।

इससे हमारा विश्लेषण पूरा होता है।
अब हमारे पास ऑब्जेक्ट डी है, जिसमें नमूना स्थित है, और जिसमें विभिन्न विधियां (फ़ंक्शन) हैं जो हमें इस नमूने के मापदंडों (औसत, सबसे संभावित मूल्य, विचरण, आदि) की गणना करने की अनुमति देते हैं।

परिणामों की तुलना करने के लिए, हम सबसे कम वर्ग विश्लेषण करते हैं:

 chisq_result = polyfit(model.x, model.data, 1) 

अब सभी परिणाम प्रिंट करें:

 print "\n\nResult of chi-square result: a= %f, b= %f" % (chisq_result[0], chisq_result[1]) print "\nResult of Bayesian analysis: a= %f, b= %f" % (Davalue, Dbvalue) print "\nThe real coefficients are: a= %f, b= %f\n" %(model.true_coefficients[0], model.true_coefficients[1]) 

हम pymc के लिए मानक ग्राफिक्स बनाते हैं:

 pymc.Matplot.plot(D) 

और अंत में, हम अपना अंतिम कार्यक्रम बना रहे हैं:

 figure() plot(model.x, model.data, marker='+', linestyle='') plot(model.x, Davalue * model.x + Dbvalue, color='g', label='Bayes') plot(model.x, chisq_result[0] * model.x + chisq_result[1], color='r', label='Chi-squared') plot(model.x, model.true_coefficients[0] * model.x + model.true_coefficients[1], color='k', label='Data') legend() show() 

यहाँ run_model.py फ़ाइल की पूर्ण सामग्री है:

 from numpy import polyfit from matplotlib.pyplot import figure, plot, show, legend import pymc import model #Define MCMC: D = pymc.MCMC(model, db = 'pickle') #Sample MCMC: 10000 iterations, burn-in period is 1000 D.sample(iter = 10000, burn = 1000) #compute chi-squared fitting for comparison: chisq_result = polyfit(model.x, model.data, 1) #print the results: print "\n\nResult of chi-square result: a= %f, b= %f" % (chisq_result[0], chisq_result[1]) print "\nResult of Bayesian analysis: a= %f, b= %f" % (Davalue, Dbvalue) print "\nThe real coefficients are: a= %f, b= %f\n" %(model.true_coefficients[0], model.true_coefficients[1]) #plot graphs from MCMC: pymc.Matplot.plot(D) #plot noised data, true line and two fitted lines (bayes and chi-squared): figure() plot(model.x, model.data, marker='+', linestyle='') plot(model.x, Davalue * model.x + Dbvalue, color='g', label='Bayes') plot(model.x, chisq_result[0] * model.x + chisq_result[1], color='r', label='Chi-squared') plot(model.x, model.true_coefficients[0] * model.x + model.true_coefficients[1], color='k', label='Data') legend() show() 

टर्मिनल में, हम निम्नलिखित उत्तर देखते हैं:

ची-स्क्वायर परिणाम का परिणाम: a = 10.321533, b = 6.307100

बायेसियन विश्लेषण का परिणाम: a = 10.366272, b = 6.068982

वास्तविक गुणांक हैं: a = 10.400000, b = 5.500000

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

और फ़ाइल run_model.py के साथ फ़ोल्डर में हम निम्नलिखित ग्राफ़ देखेंगे।
पैरामीटर के लिए एक:

पैरामीटर बी के लिए:

सिग्मा पैरामीटर के लिए:

दाईं ओर हम पीछे के वितरण का एक हिस्टोग्राम देखते हैं, और बाईं ओर दो चित्र मार्कोव श्रृंखला के हैं।
मैं अब उन पर ध्यान केंद्रित नहीं करूंगा। मैं केवल यह कह सकता हूं कि निचला ग्राफ आटोक्लेररेशन ग्राफ है (अधिक विवरण यहां पाया जा सकता है )। यह MCMC के अभिसरण का एक विचार देता है।
और ऊपरी ग्राफ नमूना ट्रेस दिखाता है। यही है, यह दिखाता है कि समय के साथ नमूना कैसे हुआ। इस ट्रेस का औसत नमूना का औसत है (दाईं ओर हिस्टोग्राम में क्षैतिज अक्ष के साथ इस ग्राफ में ऊर्ध्वाधर अक्ष की तुलना करें)।

अंत में, मैं एक और दिलचस्प विकल्प के बारे में बात करूंगा।
यदि आप फिर भी पीडॉट मॉड्यूल डालते हैं और निम्नलिखित लाइन को run_model.py फ़ाइल में शामिल करते हैं:

 pymc.graph.dag(D).write_png('dag.png') 

फिर वह फ़ोल्डर में run_model.py फ़ाइल के साथ निम्न छवि बनाएगा:



यह हमारे मॉडल का प्रतिनिधित्व करने वाला एक सीधा चक्रीय ग्राफ है। सफेद दीर्घवृत्त स्टोकेस्टिक नोड्स दिखाते हैं (ये एक, बी और सिग्मा हैं), त्रिकोण निर्धारित नोड्स दिखाते हैं, और एक अंधेरे दीर्घवृत्त में हमारे छद्म प्रयोगात्मक डेटा शामिल हैं।
यही है, हम देखते हैं कि a और b का मान हमारे मॉडल (रैखिक_फिट) में प्रवेश करता है, जो स्वयं एक नियतात्मक नोड है, और फिर संभावना फ़ंक्शन y दर्ज करें। सिग्मा पहले एक स्टोकेस्टिक नोड द्वारा निर्दिष्ट किया जाता है, लेकिन चूंकि संभावना फ़ंक्शन में पैरामीटर सिग्मा नहीं है, लेकिन ताऊ = 1 / सिग्मा ^ 2, स्टोचैस्टिक मूल्य सिग्मा पहले चुकता (दाहिनी ओर ऊपरी त्रिकोण), और फिर ताऊ माना जाता है। और पहले से ही ताऊ संभावना समारोह में प्रवेश करता है, साथ ही साथ हमारे उत्पन्न डेटा भी।
मुझे लगता है कि यह ग्राफ मॉडल की व्याख्या करने और आत्म-जाँच तर्क के लिए बहुत उपयोगी है।
मॉडल।