👩🏿‍🔬 🚤 🐣 تسوية إحصائية للمشاكل العكسية غير الصحيحة لهم. تورتشين (الجزء 1) 😏 🎾 🎅🏻

مرحبا يا هابر! نريد اليوم أن نخبرك بما يفعله مختبر طرق تجربة الفيزياء النووية ، وهو جزء من أبحاث JetBrains .

أين هي JetBrains وأين الفيزياء النووية ، تسأل. اتفقنا على أساس حب Kotlin ، على الرغم من أننا لن نتحدث عنه في هذا المنشور. تركز مجموعتنا على تطوير برامج تحليل البيانات والنمذجة وكتابة العلماء ، وبالتالي تركز على التعاون وتبادل المعرفة مع شركات تكنولوجيا المعلومات.

في هذه المقالة ، نريد أن نتحدث عن طريقة التنظيم الإحصائي التي ننشرها ، والتي اقترحها VF Turchin في السبعينيات من القرن العشرين ، وتطبيقها في شكل رمز في Python و Julia.

سيكون العرض التقديمي مفصلاً تمامًا ، لذلك يمكن لأولئك الذين هم واضحون حول المشكلات العكسية الانتقال مباشرة إلى الأمثلة وقراءة النظرية في هذه المقالة .

حدوث المشكلة: لماذا يجب على أي شخص تسوية على الإطلاق؟

إذا كان يكفي التجاهل ، يمكن وصف أي قياس في التجربة على النحو التالي: هناك جهاز معين يلتقط طيف أو إشارة عملية ويظهر بعض الأرقام وفقًا لنتائج القياس. مهمتنا كباحثين ، بالنظر إلى هذه الأرقام ومعرفة هيكل الجهاز ، هي فهم ما هو الطيف أو الإشارة المقاسة. أي في وجه ما يسمى المشكلة العكسية . إذا كنت تتخيل هذا الحساب ، نحصل على هذه المعادلة (والتي ، بالمناسبة ، تسمى معادلة فريدهولم من النوع الأول ):

f (y) = \int_{a}^{b} d x K (x, y) φ (x)

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

في الواقع ، تصف هذه المعادلة ما يلي: يتم تمثيل جهاز القياس لدينا من خلال وظيفة الأجهزة الخاصة به

K (x, y)

$K(x,y)$ ، والتي تعمل على الطيف المدروس أو إشارة دخل أخرى

φ

$\varphi$ ، حيث يلاحظ الباحث إشارة الخرج

f (y)

$f(y)$ . هدف الباحث هو استعادة الإشارة

φ

$\varphi$ من المعروف

f (y)

$f(y)$ و

K (x, y)

$K(x,y)$ . يمكنك أيضًا صياغة هذا التعبير في شكل مصفوفة ، واستبدال الدوال بالناقلات والمصفوفات:

f_{m} = K_{m n} φ_{n}

$f_m = K_{mn}\varphi_n$

يبدو أن إعادة بناء الإشارة ليست مهمة صعبة ، حيث أن لكل من معادلة فريدهولم ونظام المعادلات الخطية (حتى المفرطة التحديد) حل دقيق. لنجربها. دع الإشارة المقيسة توصف بأنها مجموع غاوسين:

φ (x) = 2 * N (2, 0.16) + N (4, 0.04)

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

كأداة ، نأخذ أبسط أداة دمج - مصفوفة تترجم إشارتنا إلى مجموع تراكمي باستخدام وظيفة Heaviside:

K_{m n} = θ (x_{m} - y_{n})

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

يظهر نوع الإشارة المقاسة وجهازنا ، وكذلك نتيجة القياس على الرسم البياني.

من المهم أن يحتوي أي قياس حقيقي على خطأ ، لذلك سنفسد نتائجنا قليلاً عن طريق إضافة ضوضاء عادية ، مع إعطاء خطأ قياس بنسبة خمسة بالمائة.

سنقوم باستعادة الإشارة باستخدام طريقة المربعات الصغرى:

φ^{М Н К} = (K^{T} K)^{- 1} K^{T} f

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

ونتيجة لذلك نحصل على: في

الواقع ، يمكننا في هذا الأمر إنهاء المقالة ، مرة أخرى بعد أن أقنعنا أنفسنا بالعجز في الأساليب المثالية للرياضيات في مواجهة الواقع المادي القاسي والقاسي ، والذهاب إلى مكواة لحام الدخان.

لكن أولاً ، دعنا نتعرف على سبب حدوث هذا الفشل لنا؟ من الواضح أن النقطة هي أخطاء القياس ، ولكن ما تأثيرها؟ والحقيقة هي أن جاك هادامارد (وهو نفس الشخص الذي أضاف اندفاعة إلى صيغة كوشي - هادامارد) قسم جميع المهام إلى مهام صحيحة وغير صحيحة.

تذكر الكلاسيكيات: "لا معنى للبحث عن حل ، إن وجد. يتعلق الأمر بكيفية التعامل مع مهمة ليس لها حل. هذا سؤال أساسي للغاية ... "- لن نتحدث عن المهام الصحيحة ونتناول المهام غير الصحيحة على الفور. لحسن الحظ ، لقد حققنا هذا بالفعل: معادلة فريدهولم المكتوبة أعلاه هي مشكلة عكسية غير صحيحة - حتى مع التقلبات الصغيرة بلا حدود في بيانات الإدخال (وحتى أخطاء القياس لدينا بعيدة عن اللانهائية) ، يمكن أن يختلف حل المعادلة التي تم الحصول عليها بطريقة تحليلية دقيقة بشكل تعسفي عن المعادلة الحقيقية .

يمكنك قراءة دليل هذا البيان في الفصل الأول من العمل الكلاسيكي للأكاديمي A.N. تيخونوفا "طرق لحل المشاكل السيئة". يحتوي هذا الكتاب على نصائح حول ما يجب القيام به مع المهام غير الصحيحة ، ولكن التقنية الموضحة هناك لها عدد من العيوب التي تم إصلاحها في طريقة Turchin. لكن أولاً ، نصف المبادئ العامة للعمل مع المهام غير الصحيحة: ماذا تفعل إذا صادفت مثل هذه المهمة؟

نظرًا لأن المهمة نفسها لا يمكنها أن تقدم لنا أي شيء ، سيكون من الضروري ارتكاب جريمة صغيرة: قم بتكملة المهمة بالبيانات حتى تصبح صحيحة ، وبعبارة أخرى ، أدخل بعض * معلومات مسبقة إضافية حول المهمة * (تسمى هذه العملية تسوية المهمة). على النقيض من طريقة تيخونوف الكلاسيكية ، استنادًا إلى إدخال وظائف تنظيمية بارامترية ، تستدعي طريقة Turchin من ناحية أخرى باستخدام طرق بايزي.

الوصف النظري للتنظيم الإحصائي

إستراتيجية

نقوم بصياغة مشكلتنا من حيث الإحصائيات الرياضية: وفقًا لتطبيق معروف جيدًا

f

$f$ (الذي نقيسه في التجربة) نحتاج إلى تقدير قيمة المعلمة

φ

$\varphi$ . وظيفي

\hat{S}

$\hat{S}$ الحوسبة

φ

$\varphi$ القائم

f

$f$ سوف نسميالإستراتيجية. من أجل تحديد الإستراتيجيات الأفضل ، نقدمدالة خسارة تربيعية. يمكن أن تكون وظيفة الخسارة الحقيقية موجودة ، لماذا نختار الوظيفة التربيعية؟ لأنه يمكن تقريب أي دالة خسارة بالقرب من الحد الأدنى لها بواسطة دالة تربيعية:

L (φ, \hat{S} [f]) = | | \hat{φ} - \hat{S} [f]) | |_{L_{2}},

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

أين

- الحل الأفضل. ثم يتم تحديد الخسائر لإستراتيجيتنا المختارة منخلال وظيفة المخاطر:

\hat{φ}

$\hat{\varphi}$

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

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

هنا

P (f | φ)

$P(f|\varphi)$ تحدد

كثافة احتمالية مجموعتنا التي يتم حساب متوسط الخسارة عليها. تتكون هذه المجموعة من تكرار متعدد للقياسات.

f

$f$ للحصول على معين

φ

$\varphi$ . في هذا الطريق،

P (f | φ)

$P(f|\varphi)$ هي كثافة الاحتمال المعروفة لنا

f

$f$ تم الحصول عليها في التجربة.

وفقًا لمنهج بايزي ، يُقترح النظر فيه

φ

$\varphi$ باعتبارهامتغير عشوائيمعكثافة الاحتمال مسبقة

P (φ)

$P(\varphi)$ ، معربا عنموثوقيةالحلول المختلفة لمشكلتنا.

P (φ)

$P(\varphi)$ يتم تحديد

بناءً على المعلومات الموجودة قبل التجربة. ثم يعتمد اختيار الإستراتيجية المثلى على تقليلالمخاطر اللاحقة:

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

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

في هذه الحالة ، الإستراتيجية المثلى معروفة جيدًا:

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

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

أين الكثافة الخلفية

P (φ | f)

$P(\varphi|f)$ تحددها نظرية بايز:

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

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

سيتيح لنا هذا النهج تحديد التباين (دالة الارتباط) للحل الناتج:

D (x_{1}, x_{2}) = E [φ (x_{1}) - \hat{S} [f] (x_{1})] [φ (x_{2}) - \hat{S} [f] (x_{2})]

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

لذا ، فقد حصلنا على الحل الأمثل لمشكلتنا من خلال إدخال كثافة مسبقة

P (φ)

$P(\varphi)$ . هل يمكننا أن نقول أي شيء عن هذا العالم من الوظائف

φ (x)

$\varphi(x)$ التي تعطى من كثافة مسبقة؟

إذا كانت الإجابة على هذا السؤال بالنفي ، فسنضطر إلى قبول كل ما هو ممكن

φ (x)

$\varphi(x)$ من المحتمل بنفس القدر والعودة إلى حل غير منتظم. لذلك يجب أن نجيب على هذا السؤال بالإيجاب.

هذا هو بالضبط ما تتكون منه طريقة التنظيم الإحصائي - تسوية الحل عن طريق تقديم معلومات إضافية مسبقة عن

φ (x)

$\varphi(x)$ . إذا كان الباحث لديه بالفعل أي معلومات مسبقة (كثافة مسبقة

P (\vec{φ})

$P(\vec{\varphi})$ ) ، يمكنه فقط حساب التكامل والحصول على الإجابة.

في حالة عدم وجود مثل هذه المعلومات ، تصف الفقرة التالية الحد الأدنى من المعلومات التي قد يمتلكها الباحث وكيفية استخدامها للحصول على حل منظم.

معلومات مسبقة

كما أظهر العلماء البريطانيون ، فإنهم في بقية العالم يحبون التفرقة. علاوة على ذلك ، إذا سأل عالم الرياضيات أسئلة حول شرعية هذه العملية ، يعتقد الفيزيائي بتفاؤل أن قوانين الطبيعة توصف بوظائف "جيدة" ، أي سلسة.

وبعبارة أخرى ، يجعلها أكثر سلاسة

φ (x)

$\varphi(x)$ ارتفاع كثافة الاحتمال المسبق. لذلك دعونا نحاول تقديم احتمال مسبق يعتمد على النعومة. للقيام بذلك ، نتذكر أن إدخال معلومات مسبقة هو بعض العنف ضد العالم ، مما يجبر قوانين الطبيعة على النظر بطريقة مريحة بالنسبة لنا.

يجب تقليل هذا العنف ، ومن خلال إدخال كثافة احتمالية مسبقة ، من الضروري أن تكون معلومات شانون عن

φ (x)

$\varphi(x)$ الواردة في

P (\vec{φ})

$P(\vec{\varphi})$ كان الحد الأدنى. إضفاء الطابع الرسمي على ما ورد أعلاه ، نستمد شكل كثافة مسبقة بناءً على سلاسة الوظيفة. للقيام بذلك ، سنبحث عن حد أقصى مشروط للمعلومات:

I [P (\vec{φ})] = \int \ln P (\vec{φ}) P (\vec{φ}) d \vec{φ} \to m i n

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

وفقا للشروط التالية:

حالة النعومة $\varphi(x)$ . اسمحوا ان $\Omega$ هي مصفوفة معينة تصف سلاسة الوظيفة. ثم نطلب أن يتم تحقيق قيمة معينة من وظائف النعومة:
$\int (\vec{φ}, Ω \vec{φ}) P (\vec{φ}) d \vec{φ} = ω$
$\int(\vec{\varphi},\Omega\vec{\varphi}) P(\vec{\varphi}) d\vec{\varphi} = \omega$
يجب على القارئ اليقظ طرح سؤال حول تحديد قيمة المعلمة.
$\omega$ . سيتم إعطاء الجواب أكثر في النص.
تطبيع الاحتمال لكل وحدة:
في ظل هذه الظروف ، ستوفر الوظيفة التالية الحد الأدنى من الوظائف:
$P_{α} (\vec{φ}) = \frac{α^{R g (Ω) / 2} det Ω^{1 / 2}}{(2 π)^{N / 2}} \exp (- \frac{1}{2} (\vec{φ}, α Ω \vec{φ}))$
معامل متصل مع ، ولكن نظرًا لأننا لا نمتلك معلومات حول القيم المحددة لوظيفة النعومة ، فإن معرفة كيفية توصيلها بالضبط أمر لا معنى له. ثم ما يجب القيام به، أنت تسأل. هنا يتم الكشف عن ثلاث طرق لك:
1. قيمة معلمة المطابقة $\alpha$ يدويًا ومن ثم الانتقال فعليًا إلى تسوية تيخونوف
2. احتساب (دمج) على كل ما هو ممكن $\alpha$ على افتراض كل ما هو ممكن $\alpha$ من المحتمل بنفس القدر
3. اختر على الأرجح $\alpha$ بكثافة الاحتمال الخلفية $P(\alpha|\vec{f})$ . هذا النهج صحيح لأنه يعطي تقريبًا جيدًا للتكامل إذا كانت البيانات التجريبية تحتوي على معلومات كافية حول $\alpha$ .

الحالة الأولى لا تهمنا كثيراً. في الحالة الثانية ، يجب أن نحسب مثل هذا التكامل القبيح هنا:

⟨ φ_{i} ⟩ = \frac{\int d φ φ_{i} P (f | φ) \int d α P (α) α^{\frac{R g (Ω)}{2}} \exp (- \frac{α}{2} (\vec{φ}, Ω \vec{φ}))}{\int d φ P (f | φ) \int d α P (α) α^{\frac{R g (Ω)}{2}} \exp (- \frac{α}{2} (\vec{φ}, Ω \vec{φ}))}

$\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$ عامل نعومة. في الواقع ، يمكننا استخدام أي عامل آخر (أو مجموعته الخطية) هنا ، فقط نعومة الوظيفة هي الشكل الأكثر وضوحًا للمعلومات المسبقة التي يمكننا استخدامها.

أخذ العينات

تحدثنا عن الدوال ، ولكن لا يمكن لأي جهاز حقيقي قياس سلسلة متصلة فحسب ، بل حتى مجموعة نقاط لا تحصى. نحن نأخذ دائمًا القياسات في مجموعة محدودة من النقاط ، لذلك نحن مضطرون إلى تنفيذ إجراءات التفريق والانتقال من المعادلة المتكاملة إلى المصفوفة. في طريقة التنظيم الإحصائي ، نمضي على النحو التالي: سوف نتحلل

φ (x)

$\varphi(x)$ على بعض نظام الوظائف

{T_{n}}

$\{T_n\}$ :

φ (x) = \sum_{n} φ_{n} T_{n} (x) .

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

وهكذا ، فإن معاملات هذا التوسع تشكل بعض المتجه

\vec{φ}

$\vec{\varphi}$ وهو متجه في الفضاء الوظيفي.

كمساحة وظيفية ، يمكننا أخذ مساحة هيلبرت ، أو ، على سبيل المثال ، مساحة كثيرات الحدود. علاوة على ذلك ، فإن اختيار الأساس في هذه المساحات محدود فقط من خلال خيالك (حاولنا العمل مع سلسلة فورييه المثلثية ، Poligandra و مكعبات التضفير ).

ثم عناصر المصفوفة

K

$K$ تحسب على أنها:

K_{m n} = (\hat{K} T_{n} (x)) (y_{m}),

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

أين

y_{m}

$y_m$ - النقاط التي أجريت فيها القياسات. عناصر المصفوفة

Ω

$\Omega$ سنحسب بالصيغة:

Ω_{i j} = \int_{a}^{b} (\frac{d^{p} T_{i} (x)}{d x}) (\frac{d^{p} T_{j} (x)}{d x}) d x,

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

أين

a

$a$ و

b

$b$ - حدود الفترة الزمنية التي تُحدد فيها الدالة

φ (x)

$\varphi(x)$ .

لإعادة حساب الأخطاء ، استخدم صيغة التشتت لمجموعة خطية من المتغيرات العشوائية:

D [φ (x)] = D [\sum_{n} φ_{n} T_{n} (x)] = \sum_{i, j} φ_{i} φ_{j} c o v (T_{i} (x), T_{j} (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)).$

يجب أن يؤخذ في الاعتبار أنه في بعض الحالات يؤدي تمثيل دالة باستخدام متجه البعد المحدود إلى فقدان جزئي أو تغيير في المعلومات. في الواقع ، يمكننا اعتبار الجبر نوعًا من التنظيم ، مهما كان ضعيفًا وغير كافٍ لتحويل مهمة غير صحيحة إلى مهمة صحيحة. ولكن ، على أي حال ، انتقلنا الآن من البحث

φ (x)

$\varphi(x)$ إلى بحث متجه

\vec{φ}

$\vec{\varphi}$ وفي القسم التالي نجدها.

حالة الضوضاء الغوسية

إن الحالة التي يتم فيها توزيع الأخطاء في التجربة وفقًا لغاوس ملحوظة
في أنه يمكن الحصول على حل تحليلي لمشكلتنا. نظرًا لأن المعلومات والأخطاء المسبقة لها شكل غوسي ، فإن منتجها يحتوي أيضًا على شكل غوسي ، ومن ثم يتم أخذ التكامل القبيح الذي كتبناه أعلاه بسهولة. سيكون الحل وخطأه كما يلي:

\vec{φ} = (K^{T} Σ^{- 1} K + α^{*} Ω)^{- 1} K^{T} Σ^{- 1^{T}} \vec{f}

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

Σ_{\vec{φ}} = (K^{T} Σ^{- 1} K + α^{*} Ω)^{- 1},

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

أين

Σ

$\Sigma$ - مصفوفة التغاير للتوزيع الغوسي متعدد الأبعاد ،

α^{*}

$\alpha^*$ - القيمة الأكثر احتمالاً للمعلمة

α

$\alpha$ ، والتي يتم تحديدها من حالة الحد الأقصى لكثافة الاحتمال الخلفي:

P (α | \vec{f}) = C^{'} α^{\frac{R g (Ω)}{2}} \sqrt{| (K^{T} Σ^{- 1} K + α Ω)^{- 1} |} \exp (\frac{1}{2} {\vec{f}}^{T} Σ^{- 1} K^{T} (K^{T} Σ^{- 1} K + α Ω)^{- 1} K^{T} Σ^{- 1^{T}} \vec{f})

$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})$

وإذا لم يكن لدي أخطاء غوسية؟

سيتم تخصيص الجزء الثاني من المقالة لهذا ، ولكن الآن دعونا نلخص جوهر المشكلة.

⟨ φ_{i} ⟩ = \frac{\int d φ φ_{i} P (f | φ) \int d α P (α) α^{\frac{R g (Ω)}{2}} \exp (- \frac{α}{2} (\vec{φ}, Ω \vec{φ}))}{\int d φ P (f | φ) \int d α P (α) α^{\frac{R g (Ω)}{2}} \exp (- \frac{α}{2} (\vec{φ}, Ω \vec{φ}))}

المشكلة الرئيسية هي أن هذا التكامل الرهيب ، أولاً متعدد الأبعاد ، وثانيًا ، في حدود لا نهائية. علاوة على ذلك ، إنه متعدد الأبعاد للغاية ، متجه

\vec{φ}

$\vec{\varphi}$ يمكن أن يكون لها أبعاد بسهولة

m = 30 - 50

$m = 30-50$ ، وأساليب الشبكة لحساب التكاملات لها تعقيد من النوع

O (n^{m})

$O(n^m)$ لذلك ، غير قابل للتطبيق في هذه الحالة. عند استخدام تكاملات متعددة الأبعاد ، يعمل تكامل مونتي كارلو بشكل جيد.

علاوة على ذلك ، نظرًا لأن حدودنا لا نهائية ، يجب علينا استخدام طريقة أخذ العينات المهمة ، ولكننا نحتاج بعد ذلك إلى تحديد وظيفة لأخذ العينات. لجعل كل شيء أكثر أتمتة ، يجب عليك استخدام سلسلة ماركوف مونت كارلو (MCMC) ، والتي يمكنها تكييف وظيفة أخذ العينات بشكل مستقل مع التكامل المحسوب. سنتحدث عن تطبيق MCMC في المقالة التالية.

جزء عملي

تمت كتابة التنفيذ الأول لطريقة التنظيم الإحصائي مرة أخرى في السبعينيات على Algol وتم استخدامه بنجاح للحسابات في فيزياء الغلاف الجوي. على الرغم من حقيقة أنه لا يزال لدينا المصادر المكتوبة يدويًا للخوارزمية ، قررنا إضافة القليل من الحداثة وتنفيذها في Python ، ثم في Julia.

بيثون

التركيب

التثبيت من خلال 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));

جوليا

كما ذكرنا ، يتطلب المزيد من التطوير لهذه التقنية تكاملاً متقدماً في مونت كارلو. يمكننا استخدام بعض الوحدات النمطية في Python (على سبيل المثال ، عملنا مع PyMC3) ، ومع ذلك ، فإننا ، من بين أمور أخرى ، نشارك في مشروع مشترك مع معهد ماكس بلانك في ميونيخ.

يسمى هذا المشروع مجموعة أدوات التحليل البايزية . هدفها هو إنشاء إطار عمل مع أدوات لطرق التحليل بايزي ، بما في ذلك أدوات MCMC في المقام الأول. يعمل الفريق الآن على الإصدار الثاني من الإطار المكتوب بلغة جوليا (الأول مكتوب بلغة سي ++ ++). تتمثل إحدى مهام مجموعتنا في إظهار إمكانات هذا الإطار باستخدام مثال التسوية الإحصائية ، لذلك قمنا بكتابة تنفيذ في جوليا .

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 eV يتم تشكيلها أثناء التشتت المزدوج ، ويتم استعادتها بواسطة نوبة ، حيث تم تحديدها بوضوح في وظيفة التركيب. قد تكون ذروة 30 eV شذوذًا إحصائيًا (الأخطاء كبيرة جدًا في هذه المرحلة) ، أو ربما تشير إلى وجود تشتت انفصالي إضافي.

استنتاجات وإعلان الجزء التالي

أخبرناك عن تقنية مفيدة يمكن تكييفها مع العديد من مهام تحليل البيانات (بما في ذلك التعلم الآلي) ، والحصول على "توافق" صادق للإجابة - الحل الأكثر عقلانية للمعادلة في مواجهة عدم اليقين الناجم عن أخطاء القياس. كمكافأة لطيفة ، نحصل على قيم لخطأ القرار. يمكن لأولئك الذين يرغبون في المشاركة في تطوير أو تطبيق طريقة التنظيم الإحصائي المساهمة في شكل رمز في Python أو Julia أو على شيء آخر.

في الجزء التالي سنتحدث عن:

باستخدام MCMC
التحلل Cholesky
كمثال عملي ، نعتبر استخدام التنظيم لمعالجة إشارة من نموذج كاشف مداري للبروتونات والإلكترونات

المراجع

أرسلت بواسطة ميخائيل زيليني ، الباحث في مختبر طرق تجربة الفيزياء النووية في البحوث JetBrains .

تسوية إحصائية للمشاكل العكسية غير الصحيحة لهم. تورتشين (الجزء 1)