G-H Filter

  |   Source

مرشح g-h:

مقاربة حدسية عبر تجربة فكرية

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

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

..ترددت للحظة وقلت لزميلتك أن القياس كان مختلفاً ، فردت أنها لم تقل أبدً أنه دقيق

الحساسات غير دقيقة . وهذا هو الدافع لعمليات الترشيح بالمجمل . وكان بالإمكان مباشرة تقديم الحلول المطورة بالنصف الأخير من هذا القرن ولكن هذه الحلول تم الوصول لها بطرح أسئلة بديهية حول ما نعرفه وكيف نعرفه

جرب مقياساً آخر

أول فكرة لحل المسألة السابقة , هي تجريب حساس أدق ، ولسوء الحظ فزميلتك تخبرك أنها جربت 10 حساسات كلها بنفس الدقة ، ولذلك تطلب منها إعطاءك حساساً آخر (B) أعطى وزناً 170 باوند , بينما الحساس الأول (A) أعطى قراءة 160 باوند. وبذلك , ماهي خياراتنا:

  • يمكننا تصديق الحساس الأول وأخذ القيمة 160

  • يمكننا تصديق الحساس الثاني واعتبار الوزن 170

  • يمكننا اختيار قيمة أقل من القيمتين

  • يمكننا اختيار قيمة أكبر من القيمتين

  • يمكننا اختيار قيمة بين القيمتين

الخياران الأولان ممكنان ولكن ليس لدينا سبب للأخذ بقراءة أحد الحساسين دون الآخر، ما داما بنفس مستوى الدقة

الخياران الثالث والرابع ، غير منطقيين ،أما الخيار الأخير فهو الأفضل ، لأن احتمال أن يكون الوزن الحقيقي أكبر من قراءة الحساس مساو لاحتمال أن يكون أقل وهكذا فمن المحتمل أكثر وقوع الوزن الحقيقي بين القرائتين

بالرياضيات ، يدعى هذا المقدار بالقيمة المتوقعة وهو موضوع سنعود له لاحقاً .

الآن أسأل نفسك : ماذا لو قسنا أنفسنا مليون مرة ؟ ..القيم ستتراوح بين القيم المتطرفة للأعلى أو الأدنى .. ولكن باختيار رقم بين قرائتي الحساسين نضمن عدم حصول خطأ كبير ، حتى مع حصول قراءات متطرفة .

من الواضح الآن أن التخمين الأفضل هو متوسط A و B أي $\frac{160+170}{2}=165$

بالنظر لهذا رسومياً ، وبافتراض أن الخطأ المفترض هو 8 باوند زيادة أو نقصان ، فإن التقاطع بين المجالين سيكون بين 160 و 170 باوند ولذلك فالقيمة المنطقية ستكون هناك

In [18]:
from kf_book.book_plots import figsize, plot_errorbars
with figsize(y=2):
    plot_errorbars([(160,8, 'A'), (170,8, 'B')],
                  xlims=(145,185), ylims=(-1,1));

إذا 165 تبدو قيمة جيداً ولكن بالواقع هناك معلومات إضافية ، فقيمة ك 161 غير ممكنة لأن المقياس B لا يمكنه إعطاء هكذا قيمة لأن مجال خطأه هو 8 فقط ,وكذلك بالنسبة ل 171

وبالتالي فأنسب قيمة ستكون بين 162 و 168 كما الشكل السابق

ولكن ماذا لو تم اخبارنا أن المقياس A ادق ب5 مرات من المقياس B , ولذلك قد نفكر بأخذ قراءة A لوحدها ، ولكن السؤال هنا: هل تفيدنا قراءة B بأي معلومات إضافية ؟

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

A = 160 , B = 170

ولكن ضمن مجال خطأ هو 3 للمقياس A و الخطأ ل B هو 9 باوند أي أكبر بثلاثة أضعاف

In [19]:
with figsize(y=2):
    plot_errorbars([(160,3,'A'), (170,9,'B')],
                  xlims=(145,185), ylims=(-1,1))

التقاطع هنا يظهر بوضوح المجال الممكن للقيمة الحقيقة ، وهذا ما تمت الإشارة إليه فلو أخذنا قراءة A لوحدها أي 160 فهناك خطأ إذ أنها تقع خارج المجال الممكن ل B ولذلك فالقراءة من B أدت لتحسين تقديراتنا رغم أنها أقل دقة

وبأخذ مثال آخر حيث المقياس A له مجال خطأ 1 باوند ، أما المقياس B فمجاله 9 باوند

إذا كانت القراءات كما السابق أي : A=160, B = 170 فما هو تقديرنا للوزن، فلنلاحظ رسومياً في التالي

In [20]:
with figsize(y=2):
    plot_errorbars([(160,1,'A'),(170,9,'B')],
                   xlims=(145,185), ylims=(-1,1))

يمكننا أن نرى هنا أن القياس الممكن الوحيد هو 161. وهذا مهم , إذ أنه من حساسين غير دقيقين نسبياً تمكننا من اشتقاق قيمة بدقة متناهية

ولذلك فحساسان , حتى لو كانا غير دقيقين , دائماً أفضل من واحد

ولكن رغم هذا ، ما من مستخدم سيرغب بشراء عدة مقاييس ليزن نفسه , ولذلك يمكن التفكير باستخدام الحساس الواحد أكثر من مرة,مثلاً لو قمنا باستخدام الحساس نفسه 1000 مرة , حيث استنتجنا أنه يمكن الوصول لقيم دقيقة باخذ المتوسط الحسابي لاعداد كبيرة من القراءات

مثلاً لنأخذ مثال المحاكاة التالي:

In [21]:
import numpy as np
measurments = np.random.uniform(160,170, size=10000)

print('Average of measurement is {:.4f}'.format(
    measurments.mean()))
Average of measurement is 164.9742

الرقم الناتج يجب أن يكون قريباً جداً من 165

البرنامج السابق يقوم بافتراض وحيد غير دقيق تماماً ، وهو أن احتمال قراءة 160 مساو لاحتمال قراءة 165، وهذا غير دقيق , إذ أن الحساسات بالعالم الواقعي تميل لاعطاء قيم قريبة من القيم الحقيقة أكثر من القيم البعيدة عنها, وهذا سيتم شرحه بفصول لاحقة .

الآن سنستخدم بدون شرح التابع numpy.random.normal() لتوليد قيم عشوائية لقراءات الحساسات مشابهة للقيم الواقعية , وسنقبل ذلك بدون برهان الآن.

In [22]:
measurments = np.random.normal(160,5, size=10000)

print('Average of measurement is {:.4f}'.format(
    measurments.mean()))
Average of measurement is 160.0414

هذا جيد ولكن هل يعقل أن يزن المستخدم نفسه ألف مرة للحصول على قراءة دقيقة ؟

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

In [23]:
import kf_book.book_plots as book_plots
book_plots.plot_hypothesis1()

كما نرى , هناك أنماط كثيرة للاحتمالات الممكنة للوزن بالمثال السابق، ولكن فلنتذكر أننا نتعامل مع وزن الإنسان وبالتالي فمن غير المعقول اكتسابة ل 20 باوند في يوم واحد فقط ليخسرها اليوم التالي، وهكذا فالافتراضات حول النظام الذي نتعامل معه تساهم في تحديد الاحتمالات الممكنة أكثر من غيرها

فلنأخذ مثالاً آخر ، حيث كانت القياسات لأيام متتالية كالتالي: 169,170,169,171,170,171,169,170,169,170 على ماذا يدل هذا ؟ هل يدل على أنك خسرت وزناً ب 1 باوند أم كسبت وزنا ب 1 باوند , أم بقي وزنك على حاله

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

In [24]:
book_plots.plot_hypothesis2()

وأيضاً يمكن لنا مشاهدة المثال التالي:

In [25]:
book_plots.plot_hypothesis3()

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

In [26]:
book_plots.plot_hypothesis4()

هذا الافتراض غير مقنع ، لأنه مامن قيمة مفردة يمكنها تفسير كل تلك القياسات


والآن لنفترض أننا كسبنا وزناً , بمقدار ما يحدده numpy وبحيث يبدو مناسباً للقياسات بتابع يدعى "ملائمة المربعات الأقل" . والشكل التالي يبين النتيجة:

In [27]:
book_plots.plot_hypothesis5()

هذا يبدو أفضل. وبالمقارنة مع الفرضية السابقة , فمن الواضح لأجل القياسات الحالية أنني كسبت وزناً باحتمال أكبر من أكون حافظت على وزني .
ولكن هل كسبت 13 باونداً , .. لا أحد يدري..
"أو ربما أحد ما يدري" .. قال زميل آخر
فلنفرض فرضاً ما لمقاربة هذا , فلنفرض أنني أكسب باونداً واحداً كل يوم ، ليس المهم كيف، أنما نعلم هذا الآن .. هذا يعني أنه باليوم الأول الذي كان القياس فيه 158 باوند سيكون التوقع لليوم الثاني هو 159 باوند

In [28]:
book_plots.plot_estimate_chart_1()

حسناً , ولكن ما الفائدة من هذا , وما الداعي أصلاً لقياس الوزن إذا عرفنا معدل زيادة الوزن اليومي , .. فلننظر مجدداً للقياس باليوم الثاني

In [29]:
book_plots.plot_estimate_chart_2()

وهكذا يبدو الاختلاف بين التنبؤ والقياس ، وهذا جيد , فلو أن القيم المتنبأ بها مساوية دائماً للقياس فما الفائدة من التنبؤ؟

الفكرة الرئيسية في الكتاب في المقطع التالي فاقرأه بتمعن

ما الذي يجب فعله إذاً , إذا أخذنا البيانات فقط من القياس فسنتجاهل التنبؤ والعكس صحيح, ربما من الممكن أخذ مزيج من قيم التنبؤ مع القياس


والدمج بين القيم , يبدو ببساطة بأخذ قيمة وسط بين القيمتين السابقتين أي بين 159 و 164.2


ولكن هل يجب أن تكون التقديرات في المنتصف؟ ربما من الأفضل أخذ دقة التنبؤ ودقة القياس بعين الاعتبار بكل مرة. وربما تختلف الدقة وبذلك نحسب تقديراً موزوناً حسب الدقة بكل مرة. فللنظر للشكل:

In [30]:
book_plots.plot_estimate_chart_3()

فلنجرب الآن أخذ تقديرات بنسبة $\frac{4}{10}$ من القياس نسبة للضجيج وهنا نحن نعبر عن اعتقاد , بأن التنبؤ أدق من القياس و هكذا فالتقدير الجديد سيعطى بالعلاقة:

$$new_estimate = prediction + \frac{4}{10} (measurement - prediction)$$

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


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

In [31]:
from kf_book import book_plots
import kf_book.gh_internal as gh
import matplotlib.pyplot as plt


weights = [158.0, 164.2, 160.3, 159.9, 162.1, 164.6, 
           169.6, 167.4, 166.4, 171.0, 171.2, 172.6]

time_step = 1.0 # day
scale_factor = 4.0/10

def predict_using_gain_guess(weight, gain_rate, do_print=True, sim_rate=0): 
    """ sim_rate allows for interactive plots in the notebook"""
    
    with book_plots.figsize(y=4):
        # store the filtered results
        estimates, predictions = [weight], []

        # most filter literature uses 'z' for measurements
        for z in weights: 
            # predict new position
            prediction = weight + gain_rate * time_step

            # update filter 
            weight = prediction + scale_factor * (z - prediction)

            # save
            estimates.append(weight)
            predictions.append(prediction)
            if do_print:
                gh.print_results(estimates, prediction, weight)

        # plot results
        gh.plot_gh_results(weights, estimates, predictions, sim_rate)
initial_guess = 160.
predict_using_gain_guess(weight=initial_guess, gain_rate=1)     
previous: 160.00, prediction: 161.00 estimate 159.80
previous: 159.80, prediction: 160.80 estimate 162.16
previous: 162.16, prediction: 163.16 estimate 162.02
previous: 162.02, prediction: 163.02 estimate 161.77
previous: 161.77, prediction: 162.77 estimate 162.50
previous: 162.50, prediction: 163.50 estimate 163.94
previous: 163.94, prediction: 164.94 estimate 166.80
previous: 166.80, prediction: 167.80 estimate 167.64
previous: 167.64, prediction: 168.64 estimate 167.75
previous: 167.75, prediction: 168.75 estimate 169.65
previous: 169.65, prediction: 170.65 estimate 170.87
previous: 170.87, prediction: 171.87 estimate 172.16

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


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

In [32]:
from contextlib import contextmanager
from kf_book.book_plots import reset_figsize

@contextmanager
def interactive_plot():
    %matplotlib notebook
    yield
    %matplotlib inline
    reset_figsize()
In [33]:
with interactive_plot():
    predict_using_gain_guess(weight=initial_guess, gain_rate=1., do_print=False, sim_rate=.25)

يمكنك رؤية الرسم يتقدم ببطء أولاً التنبؤ ثم القياس وأخيراً التقدير.
وهذه الطريقة التفاعلية بالعرض سترد بهذا الكتاب كثيراً ،ويمكنك تغيير معدل المحاكاة sim_rate هذا
بالعودة للترشيح , قد تبدو النتائج السابقة سخيفة، لأنه من الواضح أن معرفة المعدل ستحسن القيم , ولكن ماذا لو أعطينا الطريقة معدلاً خاطئاً كما التالي

In [34]:
predict_using_gain_guess(initial_guess, -1., do_print=False)

النتائج ليست مثيرة حقاً, .. أي نوع من المرشحات هذا الذي يتطلب منا أن نعرف معدل زيادة الوزن وبعد ذلك نفترض أنه ثابت.
ولكن ماذا لو قمنا بتغيير معدل زيادة الوزن من خطوة لأخرى بناء على التقدير السابق والقياس الحالي, مثلاً في اليوم الأول كان وزننا المقدر 159.8, بينما باليوم الثاني كان قياس الوزن هو 164.2 عندها الفرق هو 4.4

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

ومع ذلك يبدو أخذ القيمة 4.4 كمعدل الزيادة فوراً بعد أخذ 1 كمعدل قفزة كبيرة نسبياً ، لذلك تبدو هذه المسألة مشابهة لمسألة التوسيط وبذلك نأخذ معدل الزيادة الجديد بالدمج بين المعدل القديم والمعدل المحسوب وفق نسبة معينة , فلتكن $\frac{1}{3}$ .
وبذلك ستكون العلاقة:

$$ new gain = old gain + \frac{1}{3} \frac{measurment-predicted weight}{1 day}$$
In [35]:
from kf_book import book_plots
import kf_book.gh_internal as gh
import matplotlib.pyplot as plt

weight = 160.  # initial guess
gain_rate = -1.0 # initial guess

weights = [158.0, 164.2, 160.3, 159.9, 162.1, 164.6, 
           169.6, 167.4, 166.4, 171.0, 171.2, 172.6]

time_step = 5.
weight_scale = 4./10
gain_scale = 1./3
estimates = [weight]
predictions = []

for z in weights:
    # prediction step
    weight = weight + gain_rate*time_step
    gain_rate = gain_rate
    predictions.append(weight)
    
    # update step    
    residual = z - weight
    
    gain_rate = gain_rate + gain_scale   * (residual/time_step)
    weight    = weight    + weight_scale * residual
  
    estimates.append(weight)

with interactive_plot():
    gh.plot_gh_results(weights, estimates, predictions, time_step=0.1)

وهذا يبدو تحسيناً على ما سبق، حيث بداية كانت التقديرات الأولية سيئة ولكن حالما بدأت تتحسن التقديرات , ثبتت وتضائل الضجيج، بالنسبة للنسب السابقة فقد تم اختيارها عشوائياً وهي أسوأ من غيرها .
هناك نقطة أخرى تستحق التوضيح وهي العبارة:
gain_rate = gain_rate
ولقد كتبت هذه العبارة فقط لأوضح أن ما يتم تغييره بخطوة التنبؤ هو كل من القيمة التالية لكل المقادير من الوزن والمعدل أيضاً , وهذا السطر يمكن حذفه , فهنا افترضنا أن المعدل لن يتغير بالتنبؤ ولكن لاحقاً بطرق أخرى هذا قد يتغير.

مرشح g-h:

هذه الطريقة تعرف باسم مرشح g-h وهي عائدة بتسميتها للثابتين السابقين g لتعديل القياس وحساب القيمة المقدرة و h لحساب معدل الزيادة اليومي
هذا المرشح هو الأساس لمرشح كالمان , بعبارة أخرى مرشح كالمان هو صيغة من صيغ مرشح g-h وكذلك هو مرشح المربعات الأقل ومرشح Benedict-Bordner , حيث أن كل مرشح لديه طريقة خاصة باسناد قيم لكل من g و h ولكن فيما عدا ذلك الخوارزميات نفسها , مرشح كالمان مثلاً يغير بقيم g ,h مع كل خطوة


دعني أكرر الفكرة الأخيرة التي تلخص الكتاب:

  • عدة قيم دائماً أفضل من قيمة واحدة، لذلك لا تتجاهل أي قيمة مهما كانت غير دقيقة

  • دائماً اختر قيمة وسط بين القيم لإنشاء التقدير

  • تنبأ بالقيمة التالية ومعدل التغير بناء على القيمة الحالية المقدرة وكم نعتقد أنها ستتغير

  • القيمة التالية سيتم وضعها بين التنبؤ والقراءة التالية

فلننظر لعرض بصري للطريقة:

In [36]:
book_plots.create_predict_update_chart()

النظام system هو الشيء الذي نريد قياسه أياً يكن , وبعض المراجع تسميه المنشأة plant , هذه التسمية أتت من نظرية التحكم , وحالة النظام هي الأعدادات الحالية أو القيم التي تهمنا ,وهي في المثال السابق الوزن.
تقدير الحالة هو ناتج المرشح للقيمة الحالية للحالة ،أما نموذج العملية process model فيدل على المعادلة الرياضية التي تعطي الحالة اللاحقة بناء على الحالة السابقة ، مثلاً في المثال السابق نموذج العملية هو عملية الجمع لربح الوزن الأخير المحسوب مع الوزن الأخير و لكن نموذج العملية للسيارة المتحركة هو المسافة تساوي السرعة ضرب الزمن ,وهذا النموذج غير دقيق إذ أن السرعة قد لا تكون ثابتة ، لذلك نضيق مقداراً يدعى خطأ العملية .
خطوة التنبؤ تدعى انتشار النظام . وهي تستخدم نموذج العملية لتشكيل تقدير جديد . وبسبب خطأ العملية فهذا التقدير غير مثالي . بعض النصوص تدعو هذا أيضاً بالتطور evolution.
خطوة التحديث تدعى بتحديث القياسات، والخطوتين التنبؤ والتحديث معاً للمرشح تسميان بالخطوة أو القفزة epoch.
لتوضيح مناحي أكثر لتطبيق المرشح، لنفرض قطاراً يمر على سكة , والقطارات بطبيعتها تحتاج لعدة دقائق لتغيير سرعتها بشكل كبير لأنها كبيرة وبطيئة ، ولذلك إذا كان لدينا حساس مسافة يخبرنا أن المسافة تغيرت كثيراً عما تتوقعه علاقة السرعة الثابتة ، فنشك كثيرا بهذا الحساس, حتى مع افتراض أن السائق أغلق المكابح بتلك اللحظة , لأن القطارات تحتاج زمناً أطول لتغير سرعتها من ثانية واحدة ، وهكذا فيجب إعطاء التبؤ وزناً أكبر من القياسات بهكذا حالة.
يمكننا تمثيل الشكل الأخير لخوارزمية g-h بالتالي:

التهيئة

  1. تهيئة حالة المرشح

  2. تهيئة الاعتقاد بالحالة


التنبؤ

  1. استخدم سلوك النظام للتنبؤ بالحالة للخطوة اللاحقة

  2. عدل أعتقادك ليأخذ بالحسبان الارتياب بالتنبؤ


التحديث

  1. خذ قياساً واعتقاده المرافق حول دقته

  2. احسب الراسب بين القيمة المقدرة والقياس

  3. التقدير الجديد هو في مكان ما على خط الراسب


سنستخدم هذه الخوارزمية نفسها عبر الكتاب ، رغم وجود بعض التعديلات عليها

التسميات:

نعرض الآن التسميات هنا للمقادير , أولاً نرمز للقياسات ب Z وهذا هونفس الرمز بالكثير من المراجع , أما الدليل k فهو يرمز لترتيب الخطوة الزمنية فالقياس $z_k$ فهو القياسات في تلك الخطوة
خط عريض , يدل على مصفوفة أو شعاع , مثلاً لشعاع الحالة المكون من أكثر من عنصر نكتب x , وهذا يطبق على شعاع القياسات عند أخذ أكثر من حساس. ولمثال مقياس الوزن السابق , سيكون شعاع الحالة مكوناً من الوزن الحالي x ومعدل تغير الوزن $\dot{x}$ , كما يلي:

$$x = \begin{bmatrix} x\\ \dot{x} \end{bmatrix}$$


ولذلك الخوارزمية بسيطة. فهي تبدأ بقيمة أولية $\mathbf{x_0}$ وبعدها يتم إنشاء حلقة من التنبؤ بالحالة للخطوة اللاحقة ومن ثم أخذ القياس $z_k$ واختيار قيمة وسطى بين القياس والتنبؤ ، لإنشاء التقدير $\mathbf{x_k}$

تعميم البرنامج:

الكود السابق لبرنامج المرشح للوزن كان بأسماء متغيرات مرتبط بالمسألة ولكن إذا أردنا إنشاء تابع مستقل لحل أي مسألة , سنقوم بوضع متغيرات بتسميات عامة , لذلك نكتب:

In [ ]:
def g_h_filter(data, x0, dx, g, h, dt):
    """
    Performs g-h filter on 1 state variable with a fixed g and h.
    'data' contains the data to be filtered.
    'x0' is the initial value for our state variable
    'dx' is the initial change rate for our state variable
    'g' is the g-h's g scale factor
    'h' is the g-h's h scale factor
    'dt' is the length of the time step 
    """

التابع السابق يعيد البيانات كمصفوفة Numpy وليس قائمة . ويمكن تطبيقه على القياسات السابقة , والمطلوب كتابة البرنامج كما التالي:

In [37]:
#with interactive_plot():
#    # your solution here
#    pass

الحل والمناقشة:

In [38]:
from kf_book.gh_internal import plot_g_h_results
import matplotlib.pylab as pylab
from kf_book import book_plots

weights = [158.0, 164.2, 160.3, 159.9, 162.1, 164.6, 
           169.6, 167.4, 166.4, 171.0, 171.2, 172.6]

def g_h_filter(data, x0, dx, g, h, dt=1.):
    x_est = x0
    results = []
    for z in data:
        #prediction step
        x_pred = x_est + (dx*dt)
        dx = dx

        # update step
        residual = z - x_pred
        dx = dx    + h * (residual) / dt
        x_est  = x_pred + g * residual     
        results.append(x_est)  
    return np.array(results)

book_plots.plot_track([0, 11], [160, 172], label='Actual weight')
data = g_h_filter(data=weights, x0=160., dx=1., g=6./10, h=2./3, dt=1.)
plot_g_h_results(weights, data)

اختيار كل من g و h:

مرشحات g-h بأنواع مختلفة , أحد المراجع تحدد 11 نوعاً لها و بالتأكيد هناك أكثر , والاختلاف بينها جميعاً هو كيفية اختيار قيم كل من g,h فبعضها يثبتها والآخر كما مرشح كالمان يغيرها ديناميكياً , .. مثلاً قد يكون هناك علاقة تربط بينهما .
الموضوع هنا ليس مرشحات g-h بل كالمان , وما يهمنا أكثر هو الجانب البيزياني لها الذي لم يتم التطرق له حتى الآن.
من المفيد رؤية تأثير تغيير قيم g,h على الأداء من خلال الأمثلة , لفهم حسنات وعيوب هذا النوع من المرشحات , وهذا سيساعد على توضيح السلوك الأعقد لمرشح كالمان الذي سيلي شرحه.

تمرين : أنشأ تابع قياس

لنكتب الآن تابعاً يولد قياسات بضجيج ، ونريد أن يكون الضجيج أبيضاً وهذا يعني قيماً مختلفة عشوائية بدون أي نمط محدد ، وكما سترى لاحقاً سيكون متوسطها صفراً ولها تباين منتهي.
القيم السابقة يمكن توليدها بتابع Numpy وهو numpy.random.randn وما نريد إدخاله هو القيمة البدائية و عدد الخطوات ومقدار الزيادة بكل خطوة وكذلك مستوى الضجيج ، ولنختبره ب 30 نقطة ونرى ناتج الترشيح له بالتابع g_h_filter السابق

In [39]:
# your code here

الحل:

In [40]:
from numpy.random import randn
def gen_data(x0, dx, count, noise_factor):
    return [x0 + dx*i + randn()*noise_factor for i in range(count)]

measurements = gen_data(0, 1, 30, 1)
data = g_h_filter(data=measurements, x0=0., dx=1., dt=1., g=.2, h=0.02)
plot_g_h_results(measurements, data)

تمرين: شروط أولية سيئة

الآن اكتب محاكاة لقراءات تبدأ من 5 بمعدل زيادة 2 ومعدل ضجيج 10 ، وتستخدم g=0.2 و h=0.02 , ولكن ضع القيمة الأولية للمرشح أي قيمة x = 100, مع استخدام التابعين gen_data و g_h_filter.

In [41]:
# your code here

الحل والمناقشة:

In [42]:
zs = gen_data(x0=5., dx=2., count=100, noise_factor=10)
data = g_h_filter(data=zs, x0=100., dx=2., dt=1., g=0.2, h=0.01)
plot_g_h_results(measurements=zs, filtered_data=data)

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

تمرين: ضجيج كبير

اختبر نفس المثال ، ولكن مع عامل ضجيج 100. أزل القيمة الأولية المسببة للرنين بتحويلها من 100 ل 5.

In [43]:
# your code here

الحل والمناقشة:

In [44]:
zs = gen_data(x0=5., dx=2., count=100, noise_factor=100)
data = g_h_filter(data=zs, x0=5., dx=2., g=0.2, h=0.02)
plot_g_h_results(measurements=zs, filtered_data=data)

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

تمرين: أثر التسارع

اكتب تابع توليد بيانات بتسارع ثابت وذلك لبيانات القياسات , بمعتى آخر كل نقطة بيانات له سرعة متزايدة عن النقطة التي تسبقها , أرسم النتيجة حيث 0 هو عامل الضجيج و g = 0.2 و h = 0.02وارسم النتائج :

الحل والمناقشة:

In [45]:
import numpy as np

def gen_data(x0, dx, count, noise_factor, accel=0.):
    zs = []
    for i in range(count):
        zs.append(x0 + dx*i + np.random.randn()*noise_factor)
        dx += accel
    return zs

predictions = []
zs = gen_data(x0=10., dx=0., count=20, noise_factor=0, accel=2.)
#zs = acceldata(x0=10., dx=0., steps=20, ddx=2.)
data = g_h_filter(data=zs, x0=10., dx=0.,dt=1., g=0.2, h=0.02)
plt.xlim([0, 20])
plot_g_h_results(measurements=zs, filtered_data=data)

كل تنبؤ متأخر خلف الإشارة , وهذ منطقي إذ أن قيمة التسارع بمعادلة التنبؤ هي صفرية وبالتالي السرعة ثابتة , ولكن عندما يتم التحديث مع القياسات يتم تعديل السرعة فقط للسرعة الجديدة أما التسارع فلا يمكن تعديله , وهكذا تتغير السرعة بالقراءة الجديدة وبشكل مغاير لقيمة التنبؤ.


هذا الخطأ الأخير يدعى خطأ التأخير للجملة .وهذه مشكلة معروفة بمرشحات g-h , ولكن هناك الكثير من الدراسات والحلول المطورة حولها , التي سنتطرق لها هنا أيضاَ.

تمرين : تغيير g:

كبداية ما هو برأيك تأثير القيمة g والتي هي عامل الأخذ من قراءة الحساس وقيمة التنبؤ.

فلتجرب القيم التالية : عامل ضجيج 50 و سرعة 5 . وارسم النتائج من أجل g= 0.1, 0.4 ,0.8

الحل والنقاش:

In [46]:
np.random.seed(100)
zs = gen_data(x0=5., dx=5., count=50, noise_factor=50)
data1 = g_h_filter(data=zs, x0=0., dx=5., dt=1., g=0.1, h=0.01)
data2 = g_h_filter(data=zs, x0=0., dx=5., dt=1., g=0.4, h=0.01)
data3 = g_h_filter(data=zs, x0=0., dx=5., dt=1., g=0.8, h=0.01)

with book_plots.figsize(y=4):
    book_plots.plot_measurements(zs, lw=1, color='k')
    book_plots.plot_filter(data1, label='g=0.1', marker='s')
    book_plots.plot_filter(data2, label='g=0.4', marker='v')
    book_plots.plot_filter(data3, label='g=0.8', lw=2)
    plt.legend(loc=4)
    book_plots.set_limits([20,40], [50, 250])

من الواضح أن قيماً أقل ل g أعطت نتائج أفضل وذلك لأنه أهتمت أكثر بقيم التنبؤ واهملت قيم الضجيج الكبيرة على عكس قيم g الكبيرة , ولكن هذا أيضاً أهمل القياسات , ولذلك سندرس حالة تغير حقيقي بالقيم , بجعل الخطوة 1 ل 9 خطوات ومن ثم تغييرها ل0 .

In [47]:
zs = [5, 6, 7, 8, 9, 10, 11, 12, 13, 14]
for i in range(50):
    zs.append(14)

data1 = g_h_filter(data=zs, x0=4., dx=1., dt=1., g=0.1, h=0.01)
data2 = g_h_filter(data=zs, x0=4., dx=1., dt=1., g=0.5, h=0.01)
data3 = g_h_filter(data=zs, x0=4., dx=1., dt=1., g=0.9, h=0.01)

book_plots.plot_measurements(zs)
book_plots.plot_filter(data1, label='g=0.1', marker='s')
book_plots.plot_filter(data2, label='g=0.5', marker='v')
book_plots.plot_filter(data3, label='g=0.9', lw=2)
plt.legend(loc=4)
plt.ylim([6, 20])
Out[47]:
(6, 20)

وهذا يبين أثر إهمال الإشارة , حيث لم يتم فقط تجاهل الضجيج , بل التغييرات الفعلية بالقيمة , ربما من الأفضل أختيار g بقيمة مناسبة لكل مسألة .. بالواقع هناك طرق عديدة لاختيار g تختلف من مرشح لآخر , فمرشح Benedict-Bordner يعمل تحديداً على حل مسألة الانتقال العابر هذه.

تغيير h:

لدراسة أثر h والتي هي عامل الأخذ بالسرعة للقياسات , على حساب السرعة المتنبأ بها ,نقوم بأخذ عدة أمثلة فيها قيم مختلفة ل h


أول مثال يأخذ البيانات بدون ضجيج و يحدد قيمة h صغيرة مع قيمة أولية دقيقة جداً , الناتج كما الشكل , هو أداء ممتاز للمرشح , بعد ذلك نأخذ قيمة أولية غير دقيقة بالمرة , ونلاحظ حصول رنين كبير , قبل الاستقرار , وأخيراً نأخذ قيمة أكبر ل h = 0.5 بنفس التجربة السابقة , لنلاحظ رنين أكبر بالتردد ولكن أقل بالمطال , وكذلك يستقر بشكل أسرع قليلاً.

In [48]:
zs = np.linspace(0, 1, 50)

data1 = g_h_filter(data=zs, x0=0, dx=0., dt=1., g=.2, h=0.05)
data2 = g_h_filter(data=zs, x0=0, dx=2., dt=1., g=.2, h=0.05)
data3 = g_h_filter(data=zs, x0=0, dx=2., dt=1., g=.2, h=0.5)

book_plots.plot_measurements(zs)
book_plots.plot_filter(data1, label='dx=0, h=0.05')
book_plots.plot_filter(data2, label='dx=2, h=0.05', marker='v')
book_plots.plot_filter(data3, label='dx=2, h=0.5',  marker='s')
plt.legend(loc=1);

لا تكذب على المرشح:

يمكنك وضع قيم g,h بأي طريقة , المثال التالي لمرشح أداؤه جيد أياً كان الضجيج , تم فيه وضع قيم الثابتين لصفر , ولكن هذا سيجعل المرشح يتجاهل كلياً القياسات ، وبالتالي لا فائدة منه , وبالواقع من الممكن جعل أداء مرشحك أمثلياً لقيم محددة من البيانات ولكنه قد يفشل بشكل كبير بالعمل على بيانات العالم الواقعي .

In [49]:
zs = gen_data(x0=5., dx=.2, count=100, noise_factor=100)
data = g_h_filter(data=zs, x0=5., dx=.2, dt=1., g=0., h=0.)

book_plots.plot_measurements(zs)
book_plots.plot_filter(data, label='filter', lw=2)
plt.legend(loc=1);

مرشحات g-h و FilterPy

FilterPy هي مكتبة مفتوحة المصدر , وفيها الكثير من التوابع الجاهزة لأنماط المرشحات , يمكنك تنصيبها بسهولة بكتابة الأمر في طرفية:

pip install filterpy

مثلاً لاستخدام صيغتها لمرشح G-H نكتب:

In [50]:
from filterpy.gh import GHFilter
f = GHFilter(x=0., dx=0., dt=1., g=.8, h=.2)

عند بناء المرشح أخذنا قيمة الدخل ومعدل التغير الأوليين, ومن ثم الخطوة الزمنية و أخيراً قيم كل من g و h.
لتشغيل المرشح نأخذ قراءة الحساس ونطبق الطريقة:

In [51]:
f.update(z=1.2)
Out[51]:
(0.96, 0.24)

الطريقة السابقة , تعيد قيم كل من x و dx ولكن يمكن لك الوصول لهم مباشرة بالتالي:

In [52]:
print (f.x,f.dx)
(0.96, 0.24)

كما يمكنك تغيير g و h ديناميكياً

In [53]:
print(f.update(z=2.1, g=.85, h=.15))
(1.965, 0.375)

يمكنك كذلك تشغيل مصفوفة قياسات دفعة واحدة كالتالي:

In [54]:
print(f.batch_filter([3., 4., 5.]))
[[ 1.965   0.375 ]
 [ 2.868   0.507 ]
 [ 3.875   0.632 ]
 [ 4.9014  0.7306]]

يمكنك أيضاً تعقب عدة متغيرات , مثلاً الموضع ثلاثي الأبعاد لطائرة , كما يلي :

In [55]:
x_0  = np.array([1., 10., 100.])
dx_0 = np.array([10., 12., .2])
               
f_air = GHFilter(x=x_0, dx=dx_0, dt=1., g=.8, h=.2)
f_air.update(z=np.array((2., 11., 102.)))
print(' x =', f_air.x)
print('dx =', f_air.dx)
(' x =', array([   3.8 ,   13.2 ,  101.64]))
('dx =', array([ 8.2 ,  9.8 ,  0.56]))

الصنف GHFilterOrder يسمح لك بإنشاء مرشحات g-h بأبعاد 0و1و2 وهناك مرشح g-h-k الذي يتعقب التسارع أيضاً وكلها لديها تطبيقات مختلفة , ليست موضوع هذا الكتاب , مثلاً (عامل تقليص التباين VRF)

ويمكنك العودة لكل منها في توثيق هذه المكتبة في https://filterpy.readthedocs.org

ملخص:

أشجعك على التجريب بهذا المرشح , واختيار قيم مختلفة , وبالنهاية ستصل لنتيجة أن أختيار قيم g-h تجريبياً حتى مع فهم أثرها لن يكون أمثلياً لكل الحالات , فيجب بالنهاية تصميم المرشح وليس وضعه تجريبياً.


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

Comments powered by Disqus