गलत उलटा समस्याओं का सांख्यिकीय नियमितीकरण। ट्यूरिन (भाग 1)

हेलो, हेब्र! आज हम आपको बताना चाहते हैं कि नाभिकीय भौतिकी के प्रयोग के तरीके , जेटब्रेन रिसर्च का हिस्सा क्या है ।

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

इस लेख में हम उस सांख्यिकीय नियमितीकरण की विधि के बारे में बात करना चाहते हैं जिसे हम लोकप्रिय करते हैं , जिसे XX सदी के 70 के दशक में V.F Turchin द्वारा प्रस्तावित किया गया था, और पायथन और जूलिया में कोड के रूप में इसका कार्यान्वयन।

प्रस्तुति काफी विस्तृत होगी, इसलिए जो उलटा समस्याओं के बारे में स्पष्ट हैं वे सीधे उदाहरणों में जा सकते हैं और इस लेख में सिद्धांत पढ़ सकते हैं ।

• समस्या की घटना: किसी को नियमित क्यों करना चाहिए?
• सांख्यिकीय नियमितीकरण का सैद्धांतिक विवरण
• पर व्यावहारिक उदाहरण:
  
  • अजगर
  • जूलिया
• निष्कर्ष

समस्या की घटना: किसी को नियमित क्यों करना चाहिए?

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

$$

$$f(y) = \int \limits_a^b dx K(x,y)\varphi(x)$$

वास्तव में, यह समीकरण निम्नलिखित का वर्णन करता है: हमारे मापन उपकरण को इसके हार्डवेयर फ़ंक्शन द्वारा यहां दिखाया गया है $$$K(x,y)$, जो अध्ययन किए गए स्पेक्ट्रम या अन्य इनपुट सिग्नल पर कार्य करता है$$$\varphi$ , जिसके तहत शोधकर्ता आउटपुट संकेत का मानना है$$$f(y)$ । शोधकर्ता का लक्ष्य सिग्नल को बहाल करना है$$$\varphi$ ज्ञात द्वारा$$$f(y)$ और$$$K(x,y)$ । आप इस अभिव्यक्ति को मैट्रिक्स के रूप में भी बना सकते हैं, जिसमें वैक्टर और मैट्रिस के साथ फ़ंक्शंस को प्रतिस्थापित किया जा सकता है:

$$

$$f_m = K_{mn}\varphi_n$$

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

$$

$$\varphi(x) = 2*N(2, 0.16) + N(4, 0.04)$$

एक उपकरण के रूप में, हम सबसे सरल इंटीग्रेटर लेते हैं - एक मैट्रिक्स जो हेडवाइज फ़ंक्शन का उपयोग करके हमारे सिग्नल को एक संचयी योग में अनुवाद करता है:

$$

$$K_{mn} = \theta(x_m-y_n)$$

मापा संकेत और हमारे डिवाइस का प्रकार, साथ ही माप परिणाम ग्राफ पर दिखाया गया है।

यह महत्वपूर्ण है कि किसी भी वास्तविक माप में एक त्रुटि है, इसलिए हम पांच प्रतिशत माप त्रुटि देते हुए, सामान्य शोर जोड़कर हमारे परिणाम को थोड़ा खराब कर देंगे।

हम कम से कम वर्गों विधि द्वारा संकेत को बहाल करेंगे:

$$

$$\varphi^{} = (K^TK)^{-1}K^Tf$$

और परिणामस्वरूप हम प्राप्त करते हैं:

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

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

क्लासिक्स को याद करते हुए: “इसका कोई मतलब नहीं है कि किसी समाधान की तलाश की जाए, यदि कोई हो। यह एक ऐसे कार्य से निपटने के बारे में है जिसका कोई समाधान नहीं है। यह एक गहन मौलिक प्रश्न है ... "- हम सही कार्यों के बारे में बात नहीं करेंगे और गलत लोगों को तुरंत उठा लेंगे। सौभाग्य से, हम पहले से ही इस से मिल चुके हैं: ऊपर लिखा फ्रेडहोम समीकरण एक गलत उलटा समस्या है - यहां तक ​​कि इनपुट डेटा में असीम रूप से छोटे उतार-चढ़ाव के साथ (और यहां तक ​​कि हमारी माप त्रुटियां भी असीम से दूर हैं), सटीक विश्लेषणात्मक तरीके से प्राप्त समीकरण का समाधान मनमाने ढंग से सच से अलग हो सकता है ।

आप इस कथन का प्रमाण शिक्षाविद् ए.एन. के शास्त्रीय कार्य के पहले अध्याय में पढ़ सकते हैं। Tikhonova "बीमार समस्याओं के समाधान के लिए तरीके।" इस पुस्तक में गलत कार्यों के साथ क्या करना है पर सुझाव दिए गए हैं, लेकिन वहां उल्लिखित तकनीक में कई कमियां हैं जो ट्यूरिन विधि में तय की गई हैं। लेकिन पहले, हम गलत कार्यों के साथ काम करने के सामान्य सिद्धांतों का वर्णन करते हैं: यदि आप इस तरह के काम में आते हैं तो क्या करें?

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

सांख्यिकीय नियमितीकरण का सैद्धांतिक विवरण

रणनीति

हम गणितीय आँकड़ों के संदर्भ में अपनी समस्या तैयार करते हैं: एक प्रसिद्ध कार्यान्वयन के अनुसार $$$f$(जो हम प्रयोग में मापते हैं) हमें पैरामीटर के मूल्य का अनुमान लगाने की आवश्यकता है$$$\varphi$ । कार्यात्मक$$$\hat{S}$ कंप्यूटिंग$$$\varphi$ आधारित है$$$f$ हमरणनीतिकहेंगे। यह निर्धारित करने के लिए कि कौन सी रणनीतियां अधिक इष्टतम हैं, हम एकद्विघात हानि फ़ंक्शन कापरिचय देते हैं। असली नुकसान समारोह कोई भी हो सकता है, हम द्विघात का चयन क्यों करते हैं? क्योंकि इसके न्यूनतम के पास कोई भी हानि फ़ंक्शन द्विघात कार्य द्वारा अनुमानित किया जा सकता है:

$$

$$L(\varphi,\hat{S}[f]) = ||\hat{\varphi}-\hat{S}[f])||_{L_2},$$

कहाँ पे $$$\hat{\varphi}$- सबसे अच्छा समाधान। तब हमारी चुनी हुई रणनीति के नुकसानजोखिम समारोह द्वारानिर्धारितकिएजाते हैं:

$$

$$R_{\hat{S}[f]}(\varphi) \equiv E[L(\varphi,\hat{S}[f])] = \int L(\varphi,\hat{S}[f])P(f|\varphi)df.$$

यहाँ $$$P(f|\varphi)$हमारे कलाकारों की टुकड़ी से अधिक की प्रायिकता घनत्व जो नुकसान औसत निकाला जाता है निर्धारित करता है। यह पहनावा माप के काल्पनिक दोहराव से बनता है।$$$f$ किसी दिए हुए के लिए$$$\varphi$ । इस तरह,$$$P(f|\varphi)$ संभावना हमारे लिए जाना जाता घनत्व है$$$f$ प्रयोग में प्राप्त किया।

बायेसियन दृष्टिकोण के अनुसार, इस पर विचार करना प्रस्तावित है$$ एकप्राथमिकता संभावना घनत्व केसाथ$\varphi$एकयादृच्छिक चर केरूप में $$$P(\varphi)$, व्यक्तकी विश्वसनीयताहमारी समस्या को विभिन्न समाधान।$$$P(\varphi)$ जानकारी है कि प्रयोग से पहले से मौजूद है के आधार पर निर्धारित किया जाता है। तब इष्टतम रणनीति का चुनाव एकपश्च जोखिमको कम करने पर आधारित है:

$$

$$r_{\hat{S}}(\varphi) \equiv E_{\varphi}E_{f}[L(\varphi,\hat{S}[f])|\varphi]$$

इस मामले में, इष्टतम रणनीति सर्वविदित है:

$$

$$\hat{S}[f] = E[\varphi|f] = \int \varphi P(\varphi|f)d\varphi,$$

जहां पीछे का घनत्व है $$बायस प्रमेय द्वारा निर्धारित किया जाता है:$P(\varphi|f)$

$$

$$P(\varphi|f)= \frac{P(\varphi)P(f|\varphi)}{\int d\varphi P(\varphi)P(f|\varphi)}$$

यह दृष्टिकोण हमें परिणामी समाधान के विचरण (सहसंबंध समारोह) का निर्धारण करने की अनुमति देगा:

$$

$$D(x_1,x_2) = E[\varphi(x_1) - \hat{S}[f](x_1)][\varphi(x_2) - \hat{S}[f](x_2)]$$

इसलिए, हमने एक प्राथमिकता घनत्व पेश करके हमारी समस्या का इष्टतम समाधान प्राप्त किया है $$$P(\varphi)$। क्या हम कार्यों की दुनिया के बारे में कुछ भी कह सकते हैं$$$\varphi(x)$जो एक प्राथमिक घनत्व द्वारा दिया जाता है?

यदि इस प्रश्न का उत्तर नहीं है, तो हमें हर संभव स्वीकार करना होगा$$$\varphi(x)$समान रूप से संभावित और एक अनियमित समाधान पर लौटें। इसलिए, हमें इस सवाल का जवाब सकारात्मक में देना चाहिए।

यह वही है जो सांख्यिकीय नियमितीकरण की विधि में शामिल है - पर एक अतिरिक्त प्राथमिकता जानकारी को पेश करके एक समाधान को नियमित करना$$$\varphi(x)$। यदि शोधकर्ता के पास पहले से कोई पूर्व सूचना (एक प्राथमिकता घनत्व) है$$$P(\vec{\varphi})$), वह सिर्फ अभिन्न गणना कर सकता है और उत्तर प्राप्त कर सकता है।

यदि ऐसी कोई जानकारी नहीं है, तो अगला पैराग्राफ बताता है कि एक शोधकर्ता के पास न्यूनतम जानकारी क्या हो सकती है और नियमित समाधान प्राप्त करने के लिए इसका उपयोग कैसे किया जा सकता है।

एक पूर्व सूचना

जैसा कि ब्रिटिश वैज्ञानिकों ने दिखाया है, बाकी दुनिया में वे अंतर करना पसंद करते हैं। इसके अलावा, अगर गणितज्ञ इस ऑपरेशन की वैधता के बारे में सवाल पूछेंगे, तो भौतिक विज्ञानी आशावादी रूप से मानते हैं कि प्रकृति के नियमों को "अच्छे" कार्यों से वर्णित किया गया है, जो कि चिकनी है।

दूसरे शब्दों में, यह इसे चिकना बनाता है$$$\varphi(x)$एक उच्च प्राथमिकता संभावना घनत्व। तो आइए चिकनाई के आधार पर एक प्राथमिकता की संभावना शुरू करने की कोशिश करते हैं। ऐसा करने के लिए, हम याद करते हैं कि एक प्राथमिक सूचना का परिचय दुनिया के खिलाफ कुछ हिंसा है, जो प्रकृति के नियमों को हमारे लिए सुविधाजनक तरीके से देखने के लिए मजबूर करती है।

इस हिंसा कम किया जाना चाहिए, और एक एक प्रायोरी प्रायिकता घनत्व शुरू करने से, यह आवश्यक है कि शैनन की जानकारी पर$$$\varphi(x)$इसमें रखा $$$P(\vec{\varphi})$न्यूनतम था। उपरोक्त को औपचारिक बनाते हुए, हम फ़ंक्शन की चिकनाई के आधार पर एक प्राथमिकता घनत्व के रूप को प्राप्त करते हैं। ऐसा करने के लिए, हम सूचनाओं की एक सशर्त चरम सीमा की तलाश करेंगे:

$$

$$I[P(\vec{\varphi})] = \int \ln{P(\vec{\varphi})} P(\vec{\varphi}) d\vec{\varphi} \to min$$

निम्नलिखित शर्तों के तहत:

1. चिकनाई की स्थिति $$$\varphi(x)$। रहने दो$$$\Omega$फ़ंक्शन की चिकनाई को चिह्नित करने वाला एक निश्चित मैट्रिक्स है। फिर हमें आवश्यकता होती है कि क्रियात्मकता का एक निश्चित मूल्य प्राप्त किया जाए:
   

   $$

   $$\int(\vec{\varphi},\Omega\vec{\varphi}) P(\vec{\varphi}) d\vec{\varphi} = \omega$$
   
   $$$\omega$. .
   
2. : $$$\int P(\vec{\varphi}) d\vec{\varphi} = 1$
   :
   

   $$

   $$P_{\alpha}(\vec{\varphi}) = \frac{\alpha^{Rg(\Omega)/2}\det\Omega^{1/2}}{(2\pi)^{N/2}} \exp(-\frac{1}{2} (\vec{\varphi},\alpha\Omega\vec{\varphi}))$$
   $$$\alpha$$$$\omega$, , , , . $$$\alpha$, . :
   
   
   1. $$$\alpha$
   2. () $$$\alpha$, $$$\alpha$
      
   3. $$$\alpha$$$$P(\alpha|\vec{f})$। यह दृष्टिकोण सही है क्योंकि यह अभिन्न का एक अच्छा अनुमान देता है यदि प्रयोगात्मक डेटा के बारे में पर्याप्त जानकारी है$$$\alpha$।

पहला मामला हमारे लिए कम रुचि का है। दूसरे मामले में, हमें इस तरह के एक बदसूरत अभिन्न की गणना करनी चाहिए:

$$

$$\left\langle \varphi_i \right\rangle = \frac{\int d\varphi\, \varphi_i P(f|\varphi) \int\limits d\alpha\,P(\alpha) \alpha^{\frac{Rg(\Omega)}{2}} \exp(-\frac{\alpha}{2} (\vec{\varphi},\Omega\vec{\varphi}))}{\int d\varphi P(f|\varphi) \int\limits d\alpha\,P(\alpha) \alpha^{\frac{Rg(\Omega)}{2}} \exp(-\frac{\alpha}{2} (\vec{\varphi},\Omega\vec{\varphi}))}$$

तीसरे मामले के लिए, हम प्रयोग में गॉसियन शोर के लिए विश्लेषणात्मक रूप से अभिन्न मूल्य प्राप्त कर सकते हैं (यह अनुभाग के माध्यम से माना जाएगा)।

यह भी ध्यान दिया जाना चाहिए कि हमने कहीं भी उपयोग नहीं किया है$$$\Omega$एक चिकनाई ऑपरेटर है। वास्तव में, हम यहां किसी भी अन्य ऑपरेटर (या उनके रैखिक संयोजन) का उपयोग कर सकते हैं, बस फ़ंक्शन की चिकनाई एक प्राथमिक जानकारी का सबसे स्पष्ट रूप है जिसे हम उपयोग कर सकते हैं।

सैम्पलिंग

हमने फ़ंक्शंस के बारे में बात की थी, लेकिन कोई भी वास्तविक डिवाइस न केवल एक निरंतरता को माप सकता है, बल्कि अंकों का एक गणनीय सेट भी हो सकता है। हम हमेशा बिंदुओं के एक सीमित सेट में माप लेते हैं, इसलिए हम अभिन्न समीकरण और मैट्रिक्स एक से संक्रमण प्रक्रियाओं को अंजाम देने के लिए मजबूर होते हैं। सांख्यिकीय नियमितीकरण की विधि में, हम निम्नानुसार आगे बढ़ते हैं: हम विघटित होंगे$$$\varphi(x)$ कुछ कार्यों की प्रणाली पर $$$\{T_n\}$:

$$

$$\varphi(x) = \sum \limits_n \varphi_n T_n(x).$$

इस प्रकार, इस विस्तार के गुणांक कुछ वेक्टर बनाते हैं $$$\vec{\varphi}$जो फंक्शन स्पेस में एक वेक्टर है।

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

फिर मैट्रिक्स के तत्व$$$K$ के रूप में गणना:

$$

$$K_{mn} = (\hat{K}T_n(x))(y_m),$$

कहाँ पे $$$y_m$- जिन बिंदुओं पर माप किए गए थे। मैट्रिक्स तत्व$$$\Omega$ हम सूत्र द्वारा गणना करेंगे:

$$

$$\Omega_{ij} = \int\limits_a^b \left(\frac{d^pT_i(x)}{dx}\right)\left(\frac{d^pT_j(x)}{dx}\right)dx,$$

कहाँ पे $$$a$तथा $$$b$- अंतराल की सीमाएं जिस पर फ़ंक्शन परिभाषित किया गया है $$$\varphi(x)$।

त्रुटियों को पुनर्गणना करने के लिए, यादृच्छिक चर के रैखिक संयोजन के फैलाव सूत्र का उपयोग करें:

$$

$$D[\varphi(x)] = D[\sum \limits_n \varphi_n T_n(x)] = \sum\limits_{i,j} \varphi_i\varphi_j cov(T_i(x), T_j(x)).$$

यह ध्यान में रखा जाना चाहिए कि कुछ मामलों में, परिमित आयाम के वेक्टर का उपयोग करके एक फ़ंक्शन का प्रतिनिधित्व आंशिक नुकसान या सूचना के परिवर्तन की ओर जाता है। वास्तव में, हम बीजगणित को एक प्रकार का नियमितीकरण मान सकते हैं, हालांकि एक गलत कार्य को सही में बदलने के लिए कमजोर और अपर्याप्त। लेकिन, वैसे भी, अब हम खोज से चले गए हैं$$$\varphi(x)$ वेक्टर खोज के लिए $$$\vec{\varphi}$और अगले भाग में हम इसे ढूंढते हैं।

गाऊसी शोर का मामला

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

$$

$$\vec{\varphi} = (K^T\Sigma^{-1}K + \alpha^*\Omega)^{-1}K^T\Sigma^{-1^{T}}\vec{f}$$

$$

$$\Sigma_{\vec{\varphi}} = (K^T\Sigma^{-1}K+\alpha^*\Omega)^{-1},$$

कहाँ पे $$$\Sigma$- बहुआयामी गाऊसी वितरण के सहसंयोजक मैट्रिक्स, $$$\alpha^*$- पैरामीटर का सबसे संभावित मूल्य $$$\alpha$, जो अधिकतम पोस्टीरियर प्रायिकता घनत्व की स्थिति से निर्धारित होता है:

$$

$$P(\alpha|\vec{f}) = C'\alpha^{\frac{Rg(\Omega)}{2}}\sqrt{|(K^T\Sigma^{-1}K+\alpha\Omega)^{-1}|}\exp(\frac{1}{2} \vec{f}^T\Sigma^{-1}K^{T}(K^T\Sigma^{-1}K+\alpha\Omega)^{-1}K^T\Sigma^{-1^{T}}\vec{f})$$

और अगर मेरे पास गॉसियन त्रुटियाँ नहीं हैं?

लेख का दूसरा भाग इसके लिए समर्पित होगा, लेकिन अब हम समस्या के सार को रेखांकित करते हैं।

$$

$$\left\langle \varphi_i \right\rangle = \frac{\int d\varphi\, \varphi_i P(f|\varphi) \int\limits d\alpha\,P(\alpha) \alpha^{\frac{Rg(\Omega)}{2}} \exp(-\frac{\alpha}{2} (\vec{\varphi},\Omega\vec{\varphi}))}{\int d\varphi P(f|\varphi) \int\limits d\alpha\,P(\alpha) \alpha^{\frac{Rg(\Omega)}{2}} \exp(-\frac{\alpha}{2} (\vec{\varphi},\Omega\vec{\varphi}))}$$

मुख्य समस्या यह है कि यह भयानक अभिन्न, सबसे पहले बहुआयामी, और दूसरी बात, अनंत सीमाओं में। इसके अलावा, यह अत्यधिक बहुआयामी, वेक्टर है$$$\vec{\varphi}$ आसानी से आयाम हो सकता है $$$m = 30-50$, और अभिन्न की गणना के लिए ग्रिड विधियों में प्रकार की जटिलता है $$$O(n^m)$, इसलिए, इस मामले में अनुपयुक्त है। बहुआयामी इंटीग्रल लेते समय, मोंटे कार्लो एकीकरण अच्छी तरह से काम करता है।

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

व्यावहारिक हिस्सा

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

अजगर

स्थापना

इसके माध्यम से स्थापित करें pip:

pip install statreg

या स्रोत कोड डाउनलोड करें ।

उदाहरण

एक उदाहरण के रूप में, staregएक मैट्रिक्स और अभिन्न समीकरण के लिए डेटा को पुनर्प्राप्त करने के लिए एक मॉड्यूल का उपयोग कैसे करें पर विचार करें ।

हम आवश्यक वैज्ञानिक पैकेज आयात करते हैं।

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.integrate import quad
%matplotlib inline

हम सही संकेत निर्धारित करते हैं, जिसे हम पुनर्स्थापित करेंगे।

a = 0
b = 5
#  
phi = lambda x: 4*norm.pdf(x-2, scale=0.4) + 2*norm.pdf(x-4, scale = 0.5)
x = np.linspace(a, b,100)
plt.plot(x, phi(x));

कर्नेल और कार्यों के दृढ़ीकरण के संचालन को परिभाषित करें (नोट: np.convolutionविशेष रूप से सरणियों के लिए):

kernel = lambda x,y : np.heaviside(x-y, 1) #  
convolution =  np.vectorize(lambda y: quad(lambda x: kernel(x,y)*phi(x), a,b)[0])

हम मापा डेटा उत्पन्न करते हैं और सामान्य वितरण का उपयोग करके उन्हें शोर करते हैं:

y = np.linspace(a, b, 50)
ftrue = convolution(y)
sig = 0.05*ftrue +0.01 #  
f = norm.rvs(loc = ftrue, scale=sig)
plt.errorbar(y, f, yerr=sig);

हम अभिन्न समीकरण को हल करते हैं

हम विवेक के लिए सॉल्वर और सहायक वर्ग का आयात करते हैं:

from statreg.model import GaussErrorUnfolder
from statreg.basis import CubicSplines

विवेक के लिए एक कार्यात्मक आधार के रूप में, हम क्यूबिक स्प्लिन का उपयोग करते हैं, और एक अतिरिक्त स्थिति के रूप में, हम इंगित करते हैं कि फ़ंक्शन किनारों पर शून्य मान लेता है।

basis = CubicSplines(y, boundary='dirichlet')
model = GaussErrorUnfolder(basis, basis.omega(2))

हम समीकरण हल करते हैं:

phi_reconstruct = model.solve(kernel, f, sig, y)

हम एक समाधान अनुसूची का निर्माण कर रहे हैं:

plt.plot(x,phi(x))
phir = phi_reconstruct(x)
phiEr = phi_reconstruct.error(x)
plt.plot(x, phir, 'g')
plt.fill_between(x, phir-phiEr, phir + phiEr, color='g', alpha=0.3);

हम मैट्रिक्स समीकरण को हल करते हैं

हम विवेक के लिए सॉल्वर और सहायक वर्ग का आयात करते हैं:

from statreg.model import GaussErrorMatrixUnfolder
from statreg.basis import CubicSplines

मैट्रिसेस प्राप्त करने के लिए, हम अपने कार्यात्मक आधार का उपयोग करते हैं, लेकिन यह स्पष्ट है कि मैट्रिस को किसी भी तरह से प्राप्त किया जा सकता है।

cubicSplines = CubicSplines(y, boundary='dirichlet')
omega = cubicSplines.omega(2)
Kmn = cubicSplines.discretizeKernel(kernel,y)

हम मैट्रिक्स समीकरण को हल करते हैं:

model = GaussErrorMatrixUnfolder(omega)
result = model.solve(Kmn, f, sig)

चार्ट बनाएँ:

phir = lambda x: sum([p*bf(x) for p, bf in zip(result.phi,cubicSplines.basisFun)])
plt.plot(x,phir(x))
plt.plot(x,phi(x));

जूलिया

जैसा कि हमने उल्लेख किया है, तकनीक के आगे विकास के लिए उन्नत मोंटे कार्लो एकीकरण की आवश्यकता है। हम पायथन में कुछ मॉड्यूल का उपयोग कर सकते हैं (उदाहरण के लिए, हमने PyMC3 के साथ काम किया है), हालांकि, हम अन्य चीजों के अलावा, म्यूनिख में मैक्स प्लैंक इंस्टीट्यूट के साथ एक संयुक्त परियोजना में भाग ले रहे हैं।

इस परियोजना को बायेसियन विश्लेषण टूलकिट कहा जाता है । इसका लक्ष्य बेयसियन विश्लेषण के तरीकों के लिए उपकरण के साथ एक ढांचा तैयार करना है, जिसमें मुख्य रूप से एमसीएमसी के उपकरण शामिल हैं। अब टीम रूपरेखा के दूसरे संस्करण पर काम कर रही है, जो जूलिया में लिखा गया है (पहला खराब C ++ में लिखा गया है)। हमारे समूह का एक कार्य सांख्यिकीय नियमितीकरण के उदाहरण का उपयोग करके इस ढांचे की क्षमताओं का प्रदर्शन करना है, इसलिए हमने जूलिया में एक कार्यान्वयन लिखा ।

using PyCall
include("../src/gauss_error.jl")
include("../src/kernels.jl")

a = 0.
b = 6.

function phi(x::Float64)
    mu1 = 1.
    mu2 = 4.
    n1 = 4.
    n2 = 2.
    sig1 = 0.3
    sig2 = 0.5

    norm(n, mu, sig, x) = n / sqrt(2 * pi*sig^2) * exp(-(x - mu)^2 / (2 * sig^2))
    return norm(n1, mu1, sig1, x) + norm(n2, mu2, sig2, x)
end
x = collect(range(a, stop=b, length=300))

import PyPlot.plot

myplot = plot(x, phi.(x))
savefig("function.png", dpi=1000)

इस बार हम एक अलग कोर का उपयोग करते हैं, हम एक एकीकृत कदम नहीं उठाएंगे, लेकिन गॉसियन के साथ एक दृढ़ विश्वास है, जो वास्तव में हमारे डेटा के लिए एक निश्चित "धुंधला" होता है:

function kernel(x::Float64, y::Float64)
    return getOpticsKernels("gaussian")(x, y)
end

convolution = y -> quadgk(x -> kernel(x,y) * phi(x), a, b, maxevals=10^7)[1]
y = collect(range(a, stop = b, length=50))
ftrue = convolution.(y)
sig = 0.05*abs.(ftrue) +[0.01 for i = 1:Base.length(ftrue)]
using Compat, Random, Distributions
noise = []
for sigma in sig
    n = rand(Normal(0., sigma), 1)[1]
    push!(noise, n)
end
f = ftrue + noise
plot(y, f)

इसी प्रकार, हम निर्धारित सिरों के साथ छींटों का आधार लेते हैं:

basis = CubicSplineBasis(y, "dirichlet")
Kmn = discretize_kernel(basis, kernel, y)
model = GaussErrorMatrixUnfolder([omega(basis, 2)], "EmpiricalBayes", nothing, [1e-5], [1.], [0.5])
result = solve(model, Kmn, f, sig)
phivec = PhiVec(result, basis)

x = collect(range(a, stop=b, length=5000))
plot(x, phi.(x))

phi_reconstructed = phivec.phi_function.(x)
phi_reconstructed_errors = phivec.error_function.(x)

plot(x, phi_reconstructed)
fill_between(x, phi_reconstructed - phi_reconstructed_errors, phi_reconstructed + phi_reconstructed_errors, alpha=0.3)

वास्तविक दुनिया उदाहरण

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

यह है कि प्रारंभिक एकीकृत स्पेक्ट्रम कैसा दिखता है:

और इसलिए - बहाली का परिणाम:

फिट के साथ विश्लेषण के तीन मुख्य नुकसान हैं:

• , .
• , .
• .

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

निष्कर्ष और अगले भाग की घोषणा

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

अगले भाग में हम बात करेंगे:

• MCMC का उपयोग करना
• चोल्स्की अपघटन
• एक व्यावहारिक उदाहरण के रूप में, हम प्रोटॉन और इलेक्ट्रॉनों के कक्षीय डिटेक्टर के एक मॉडल से एक संकेत को संसाधित करने के लिए नियमितीकरण के उपयोग पर विचार करते हैं
  

संदर्भ

• मूल काम ट्यूरिन
• वायुमंडलीय भौतिकी में नियमितीकरण के आवेदन पर शोध प्रबंध
• सांख्यिकीय नियमितीकरण के आवेदन पर हमारा प्रकाशन
• पायथन कार्यान्वयन
• जूलिया पर कार्यान्वयन
• पुरातत्व के प्रेमियों के लिए ALGOL-60 पर कार्यान्वयन
• बायेसियन विश्लेषण टूलकिट

द्वारा प्रकाशित किया गया मिखाइल Zeleny , शोधकर्ता परमाणु भौतिकी प्रयोग के तरीके की लेबोरेटरी में में जेटब्रेन्स रिसर्च ।