خوارزمية الزقورة

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

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

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

تُستخدم خوارزمية الزقورة لتوليد قيم عينة ذات توزيع طبيعي . (تم عرض القيم الموجبة فقط للتبسيط). النقاط الوردية هي في البداية أرقام عشوائية موزعة توزيعًا منتظمًا. يتم أولًا تقسيم دالة التوزيع المطلوبة إلى مناطق متساوية "A". يتم اختيار طبقة واحدة i عشوائيًا من المصدر المنتظم على اليسار. ثم تُضرب قيمة عشوائية من المصدر العلوي في عرض الطبقة المختارة، ويتم اختبار النتيجة x لمعرفة المنطقة التي تقع فيها داخل الطبقة، مع وجود 3 نتائج محتملة: 1) (اليسار، المنطقة السوداء الصلبة) تقع العينة بوضوح أسفل المنحنى ويمكن إخراجها فورًا، 2) (اليمين، المنطقة المخططة عموديًا) قد تقع قيمة العينة أسفل المنحنى، ويجب اختبارها بشكل أكبر. في هذه الحالة، يتم توليد قيمة y عشوائية داخل الطبقة المختارة ومقارنتها بـ f(x) . إذا كانت أقل، فإن النقطة تقع أسفل المنحنى ويتم إخراج القيمة x . إذا لم تكن كذلك (الحالة الثالثة)، يتم رفض النقطة المختارة x ويتم إعادة تشغيل الخوارزمية من البداية.

نظرية التشغيل

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

يتكون التوزيع الذي تختاره خوارزمية الزقورة من n منطقة متساوية المساحة؛ن-1{\displaystyle n-1}مستطيلات تغطي الجزء الأكبر من التوزيع المطلوب، فوق قاعدة غير مستطيلة تشمل ذيل التوزيع.

بافتراض دالة كثافة احتمالية متناقصة بشكل رتيبو(x){\displaystyle f(x)}، محددة للجميعx0{\displaystyle x\geq 0}تُعرَّف قاعدة الزقورة بأنها جميع النقاط داخل التوزيع وتحته.y1=و(x1){\displaystyle y_{1}=f(x_{1})}يتكون هذا من منطقة مستطيلة من(0،0){\displaystyle (0,0)}ل(x1،y1){\displaystyle (x_{1},y_{1})}والذيل (الذي يكون عادةً لانهائيًا) للتوزيع، حيثxx1(و yy1){\displaystyle x\geq x_{1}({\text{and }}y\leq y_{1})}.

تبلغ مساحة هذه الطبقة (لنسميها الطبقة 0) A. أضف فوقها طبقة مستطيلة بعرضx1{\displaystyle x_{1}}والارتفاعأ/x1{\displaystyle A/x_{1}}لذلك فهي تمتلك مساحة أيضًاأ{\displaystyle A}أعلى هذه الطبقة على ارتفاعy2=y1+أ/x1{\displaystyle y_{2}=y_{1}+A/x_{1}}ويتقاطع مع دالة الكثافة عند نقطة(x2،y2){\displaystyle (x_{2},y_{2})}، أينy2=و(x2){\displaystyle y_{2}=f(x_{2})}تتضمن هذه الطبقة كل نقطة في دالة الكثافة بينy1{\displaystyle y_{1}}وy2{\displaystyle y_{2}}ولكنها (على عكس الطبقة الأساسية) تتضمن أيضًا نقاطًا مثل(x1،y2){\displaystyle (x_{1},y_{2})}والتي لا تقع ضمن التوزيع المطلوب.

ثم تُضاف طبقات أخرى فوقها. لاستخدام جدول مُعدّ مسبقًا بحجمن{\displaystyle n}( ن  =  256 هو العدد النموذجي)، يختار المرءx1{\displaystyle x_{1}}بحيثxن=0{\displaystyle x_{n}=0}وهذا يعني أن الصندوق العلوي، الطبقةن-1{\displaystyle n-1}، يصل إلى ذروة التوزيع عند(0،و(0)){\displaystyle (0,f(0))}بالضبط.

طبقةأنا{\displaystyle i}يمتد عموديًا في المدى[yأنا،yأنا+1]{\displaystyle [y_{i},y_{i+1}]}ويمكن تقسيمها أفقيًا إلى منطقتين: الجزء (الأكبر عمومًا) في النطاق[0،xأنا+1]{\displaystyle [0,x_{i+1}]}والتي تقع بالكامل ضمن التوزيع المطلوب، والجزء (الصغير) في النطاق[xأنا+1،x]{\displaystyle [x_{i+1},x]}، وهو ما يتم احتواؤه جزئياً فقط.

مع تجاهل مشكلة الطبقة 0 للحظة، وبافتراض وجود متغيرات عشوائية منتظمةيو0{\displaystyle U_{0}}ويو1[0،1){\displaystyle U_{1}\in [0,1)}يمكن وصف خوارزمية الزقورة على النحو التالي:

  1. اختر طبقة عشوائية0أنا<ن{\displaystyle 0\leq i<n}.
  2. يتركx=يو0xأنا{\displaystyle x=U_{0}x_{i}}.
  3. لوx<xأنا+1{\displaystyle x<x_{i+1}}، يعودx{\displaystyle x}.
  4. يتركy=yأنا+يو1(yأنا+1-yأنا){\displaystyle y=y_{i}+U_{1}(y_{i+1}-y_{i})}.
  5. الحوسبةو(x){\displaystyle f(x)}. لوy<و(x){\displaystyle y<f(x)}، يعود x{\displaystyle x}.
  6. وإلا، فاختر أرقامًا عشوائية جديدة وارجع إلى الخطوة 1.

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

مع الطبقات المتقاربة، تتوقف الخوارزمية عند الخطوة 3 في نسبة كبيرة جدًا من الحالات. بالنسبة للطبقة العليان-1{\displaystyle n-1}لكن هذا الاختبار يفشل دائمًا، لأنxن=0{\displaystyle x_{n}=0}.

يمكن تقسيم الطبقة 0 أيضًا إلى منطقة مركزية وحافة، لكن الحافة عبارة عن ذيل لانهائي. لاستخدام نفس الخوارزمية للتحقق مما إذا كانت النقطة تقع في المنطقة المركزية، قم بإنشاء قيمة وهمية.x0=أ/y1{\displaystyle x_{0}=A/y_{1}}سيؤدي هذا إلى توليد نقاط معx<x1{\displaystyle x<x_{1}}بالتردد الصحيح، وفي الحالة النادرة التي يتم فيها اختيار الطبقة 0 وxx1{\displaystyle x\geq x_{1}}استخدم خوارزمية احتياطية خاصة لاختيار نقطة عشوائيًا من الطرف. ولأن هذه الخوارزمية تُستخدم أقل من مرة واحدة في الألف، فإن السرعة ليست ضرورية.

وبالتالي، فإن خوارزمية الزقورة الكاملة للتوزيعات أحادية الجانب هي:

  1. اختر طبقة عشوائية0أنا<ن{\displaystyle 0\leq i<n}.
  2. يتركx=يو0xأنا{\displaystyle x=U_{0}x_{i}}.
  3. لوx<xأنا+1{\displaystyle x<x_{i+1}}، يعودx{\displaystyle x}.
  4. لوأنا=0{\displaystyle i=0}قم بإنشاء نقطة من الذيل باستخدام خوارزمية التراجع.
  5. يتركy=yأنا+يو1(yأنا+1-yأنا){\displaystyle y=y_{i}+U_{1}(y_{i+1}-y_{i})}.
  6. الحوسبةو(x){\displaystyle f(x)}. لوy<و(x){\displaystyle y<f(x)}، يعود x{\displaystyle x}.
  7. وإلا، فاختر أرقامًا عشوائية جديدة وارجع إلى الخطوة 1.

في التوزيع ثنائي الجانب، يجب عكس النتيجة بنسبة 50% من الوقت. ويمكن القيام بذلك بسهولة في كثير من الأحيان عن طريق اختياريو0[-1،1]{\displaystyle U_{0}\in [-1,1]}وفي الخطوة الثالثة، يتم اختبار ما إذا|x| <xأنا+1{\displaystyle \mid x\mid <x_{i+1}}.

خوارزميات احتياطية للذيل

لأن خوارزمية الزقورة لا تُنتج معظم المخرجات إلا بسرعة كبيرة، وتتطلب خوارزمية احتياطية كلماx>x1{\displaystyle x>x_{1}}يكون الأمر دائمًا أكثر تعقيدًا من التنفيذ المباشر. وتعتمد خوارزمية التراجع المحددة على التوزيع.

في التوزيع الأسي، يكون ذيل التوزيع مطابقًا لجسمه. إحدى الطرق هي اللجوء إلى أبسط خوارزمية: E =  −ln  ( U1 ) وتعيين x = x1 − ln( U1 ) . طريقة أخرى هي استدعاء خوارزمية الزقورة بشكل متكرر وإضافة x1 إلى النتيجة.    

بالنسبة للتوزيع الطبيعي، يقترح مارساجليا خوارزمية مختصرة:

  1. ليكن x = −ln( U 1 )/ x 1 .
  2. ليكن y = −ln( U 2 ).
  3. إذا كان 2y > x2 ، فأرجع x  + x1 . 
  4. وإلا، فارجع إلى الخطوة 1.

بما أن x1 3.5 لأحجام الطاولات النموذجية، فإن الاختبار في الخطوة 3 ينجح في أغلب الأحيان. وبما أن −ln( U1 ) متغير يتبع التوزيع الأسي، فيمكن استخدام تطبيق للتوزيع الأسي.

التحسينات

يمكن تنفيذ الخوارزمية بكفاءة باستخدام جداول محسوبة مسبقًا لـ x i و y i = f ( x i )، ولكن هناك بعض التعديلات لجعلها أسرع:

  • لا يعتمد أي شيء في خوارزمية الزقورة على كون دالة توزيع الاحتمالية مُعَيَّرة (التكامل تحت المنحنى يساوي 1)، ويمكن أن يؤدي إزالة ثوابت التطبيع إلى تسريع حساب f ( x ).
  • تعتمد معظم مولدات الأرقام العشوائية المنتظمة على مولدات الأرقام العشوائية الصحيحة التي تُرجع عددًا صحيحًا في النطاق [0،  2 ^32 - 1]. يسمح لنا جدول 2 ^32 x i باستخدام هذه الأرقام مباشرةً لـ U 0 .
  • عند حساب التوزيعات ذات الجانبين باستخدام U 0 ذات الجانبين كما هو موضح سابقًا، يمكن تفسير العدد الصحيح العشوائي على أنه عدد موقع في النطاق [−2 31 ،  2 31 − 1]، ويمكن استخدام عامل قياس 2 −31 .
  • بدلاً من مقارنة U₀xᵢ بـ xᵢ + 1 في الخطوة 3 ، يمكن حساب xᵢ + 1 / xᵢ مسبقًا ومقارنة U₀ بهذه القيمة مباشرةً. إذا كان U₀ مولد أرقام عشوائية صحيحة، فيمكن ضرب هذه الحدود مسبقًا في 2³² ( أو 2³¹ ، حسب الاقتضاء) بحيث يمكن استخدام مقارنة الأعداد الصحيحة .
  • مع التغييرين المذكورين أعلاه، لم يعد جدول قيم x i غير المعدلة مطلوبًا ويمكن حذفه.
  • عند توليد قيم الفاصلة العائمة أحادية الدقة وفقًا لمعيار IEEE 754 ، والتي تحتوي على جزء كسري بطول 24 بت فقط (بما في ذلك الرقم 1 الضمني في البداية)، لا تُستخدم البتات الأقل أهمية من عدد صحيح عشوائي بطول 32 بت. يمكن استخدام هذه البتات لتحديد رقم الطبقة. (انظر المراجع أدناه لمناقشة مفصلة لهذا الموضوع).
  • يمكن وضع الخطوات الثلاث الأولى في دالة مضمنة ، والتي يمكنها استدعاء تنفيذ خارج الخط للخطوات الأقل استخدامًا.

إنشاء الجداول

من الممكن تخزين الجدول بأكمله محسوبًا مسبقًا، أو تضمين القيم n و y 1 و A وتنفيذ f −1 ( y ) في الكود المصدري ، وحساب القيم المتبقية عند تهيئة مولد الأرقام العشوائية.

كما سبق شرحه ، يمكننا إيجاد xᵢ = f⁻¹ ( yᵢ ) و yᵢ⁺¹ = yᵢ + A / xᵢ . نكرر هذه العملية n⁻¹ مرة لطبقات الزقورة . في النهاية، يجب أن نحصل على yₙ = f ( 0 ) . سيكون هناك خطأ تقريبي بسيط ، لكن من المفيد التحقق من أن هذا الخطأ صغير بشكل مقبول.       

عند ملء قيم الجدول فعليًا، افترض فقط أن x n  =  0 و y n  = f (0)، واقبل الفرق الطفيف في مساحة الطبقة n − 1 كخطأ تقريب.   

إيجاد x 1 و A

بفرض قيمة ابتدائية (تخمين) x₁ ، نحتاج إلى طريقة لحساب مساحة t للذيل حيث x > x₁ . بالنسبة للتوزيع الأسي، تكون هذه المساحة هي e⁻ⁿx₁ ، بينما بالنسبة للتوزيع الطبيعي، بافتراض استخدامنا للدالة غير المعيارية f ( x ) = e⁻ⁿx² / 2 ، تكون المساحة هيπ2erfc(x2){\displaystyle {\sqrt {\frac {\pi }{2}}}{\text{erfc}}({\frac {x}{\sqrt {2}}})}√π / 2 erfc ( x /2 ). بالنسبة للتوزيعات الأكثر تعقيدًا،قد يكون التكامل العددي مطلوبًا.

مع هذا في متناول اليد، من x 1 ، يمكننا إيجاد y 1 = f ( x 1 )، والمساحة t في الذيل، ومساحة الطبقة الأساسية A = x 1 y 1  + t . 

ثم احسب المتسلسلتين yᵢ و xᵢ كما سبق. إذا كانت yᵢ > f ( 0 ) لأي قيمة لـ i < n، فإن التقدير الأولي x₁ كان منخفضًا جدًا ، مما أدى إلى مساحة A كبيرة جدًا . أما إذا كانت yₙ < f ( 0) ، فإن التقدير الأولي x₁ كان مرتفعًا جدًا.

بناءً على ذلك، استخدم خوارزمية لإيجاد الجذور (مثل طريقة التنصيف ) لإيجاد قيمة x₁ التي تجعل yₙ₋₁ أقرب ما يمكن إلى f ( 0). أو ابحث عن القيمة التي تجعل مساحة الطبقة العليا، xₙ₋₁ ( f ( 0) - yₙ₋₁ )، أقرب ما يمكن إلى القيمة المطلوبة A. هذا يوفر حسابًا واحدًا لـ fₙ₋₁( x ) ، وهو في الواقع الشرط الأكثر أهمية .  

تنويع ماكفارلاند

اقترح كريستوفر د. ماكفارلاند نسخة محسّنة بشكل أكبر. [ 1 ] وتطبق هذه النسخة ثلاثة تغييرات خوارزمية، على حساب جداول أكبر قليلاً.

أولاً، في الحالة الشائعة ، تُؤخذ الأجزاء المستطيلة فقط في الاعتبار، من (0، yᵢ - 1 ) إلى ( xᵢ ، yᵢ ). أما المناطق ذات الأشكال غير المنتظمة على يمين هذه الأجزاء (معظمها شبه مثلثية ، بالإضافة إلى الذيل) فتُعالج بشكل منفصل. هذا يُبسط ويُسرع المسار السريع للخوارزمية . 

ثانيًا، تُستخدم المساحة الدقيقة للمناطق ذات الأشكال غير المنتظمة؛ ولا يتم تقريبها لتشمل المستطيل بأكمله إلى ( xᵢ - 1 , yᵢ ) . وهذا يزيد من احتمالية استخدام المسار السريع. 

إحدى النتائج الرئيسية لذلك هي أن عدد الطبقات أقل بقليل من n . على الرغم من أن مساحة الأجزاء غير المنتظمة الشكل تُحسب بدقة، إلا أن المجموع الكلي يتجاوز مساحة طبقة واحدة. يتم تعديل مساحة كل طبقة بحيث يكون عدد الطبقات المستطيلة عددًا صحيحًا. إذا تجاوزت القيمة الابتدائية 0  i < n عدد الطبقات المستطيلة، تنتقل المرحلة الثانية.   

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

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

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

إذا كانت الدالة محدبة (كما هو الحال في التوزيع الأسي في كل مكان، والتوزيع الطبيعي عندما تكون قيمة | x |  أكبر من  1)، فإنها تقع ضمن المثلث السفلي. يتم اختيار انحرافين منتظمين متساويين U1 و U2 ، وقبل تحويلهما إلى المستطيل المحيط بالمنطقة غير المنتظمة، يتم اختبار مجموعهما. إذا كان U1 + U2 > 1 ، فإن النقطة تقع في المثلث العلوي ويمكن عكسها إلى (1 - U1 , 1 - U2 ) . أما إذا كان U1 + U2 < 1 - ε ، حيث ε قيمة مناسبة للتسامح ، فإن النقطة تقع أسفل المنحنى ويمكن قبولها فورًا. فقط بالنسبة للنقاط القريبة جدًا من القطر ، يلزم حساب دالة التوزيع f ( x ) لإجراء اختبار رفض دقيق. (من الناحية النظرية، يجب أن يعتمد التسامح ε على الطبقة، ولكن يمكن استخدام قيمة قصوى واحدة على جميع الطبقات مع فقدان طفيف).         

إذا كانت الدالة مقعرة (كما هو الحال في التوزيع الطبيعي لـ | x | < 1)، فإنها تشمل جزءًا صغيرًا من المثلث العلوي بحيث يكون الانعكاس مستحيلاً، ولكن يمكن قبول النقاط التي تحقق إحداثياتها المعيارية U 1 + U 2 ≤ 1 على الفور، ويمكن رفض النقاط التي تحقق U 1 + U 2 > 1 + ε على الفور.  

في الطبقة الوحيدة التي تقع فيها | x |  =  1، يكون للتوزيع الطبيعي نقطة انعطاف ، ويجب تطبيق اختبار الرفض الدقيق إذا كان 1− ε < U 1 + U 2 < 1+ ε .     

يتم التعامل مع الذيل كما هو الحال في خوارزمية Ziggurat الأصلية، ويمكن اعتباره حالة رابعة لشكل المنطقة ذات الشكل الغريب إلى اليمين.

مراجع

  1. ماكفارلاند، كريستوفر د. (24 يونيو 2015). "خوارزمية زقورة مُعدّلة لتوليد أرقام شبه عشوائية موزعة أُسّيًا وطبيعيًا" . مجلة الحساب الإحصائي والمحاكاة . 86 (7): 1281-1294 . arXiv : 1403.6870 . doi : 10.1080 /00949655.2015.1060234 . PMC 4812161. PMID 27041780. مؤرشف من الأصل في 22 يونيو 2024. تم الاسترجاع في 22 يونيو 2024 .   تجدر الإشارة إلى أن مستودع Bitbucket المذكور في الورقة البحثية لم يعد متاحًا، ويمكن الوصول إلى الكود الآن على الرابط التالي: https://github.com/cd-mcfarland/fast_prng