למה Octave עדיף על Matlab בהגרלה מהתפלגות נורמלית

אחרי שסקרתי השבוע שלוש חבילות לתכנות מדעי שיכולות להוות תחליפי Matlab, אראה היום יתרון אחד של Octave על-פני החלופה הקניינית.

הכל התחיל לפני כחודשיים כשניגש אלי חבר לעבודה, ע.א., וביקש להבין אנומליה שקיבל באחת הסימולציות. אחרי שקיבל תוצאות לא הגיוניות, הוא החל לחקור את מקור הבעיה והגיע לסימולציה של מדידות שמזינות את המערכת. מסתבר שבמציאות המדידות מתפלגות לפי התפלגות נורמלית, ולכן השתמש בפונקציית randn של Matlab כדי לייצר את סט הנתונים שמהווה קלט לסימולציה.

ע.א. זה - ברנש יסודי. כדי שיוכל לתחקר את התוצאות במקרה שמשהו ישתבש, הוא קרא לפונקציית randn בעזרת seed (גרעין? תרשו לי להמשיך במונח הלועזי). בכל לולאה הוא יצר מספר אחד אקראי מהתפלגות נורמלית. את ה-seed הוא יצר מהשעון בעזרת פונקציית clock, ושמר אותו בתוך מערך לצורכי שחזור. ואכן, למרות שהשתמש בפונקציית randn שהממוצע שלה הוא אפס, הוא קיבל הגרלה שבה 90% מהמדידות היו חיוביות, ורק 10% שליליות.

תחילה חשדתי שהבעיה נעוצה באלגוריתם שמייצר את המספר האקראי. ב-Mathworks כותבים שהם משתמשים באלגוריתם Zigguart (מימוש שלו ב-C תוכלו למצוא פה). אלא שאז גיליתי שהבעיה לא חוזרת ב-Octave, למרות שגם ב-Octave משתמשים באותו אלגוריתם. אחרי נסיונות לבודד את הבעיה אני חושב שאני יכול לקבוע שמדובר במימוש שונה, והתוצאות מראות שהמימוש של Octave מניב תוצאות טובות יותר.

כדי לאשש את ההנחה שלי שמדובר במימוש האלגוריתם, כתבתי את הסקריפט הבא:

% Ranncheck.m
% Run over different seeds, and draw a number from a normal distribution

DrawNum = 2500; % Number of draws
Seed = (1:DrawNum)/2;
x = zeros(1, DrawNum); % the random variable
for i = 1:DrawNum
randn('state', Seed(i));
x(i) = randn;
end

figure;
plot(Seed, x);
title('Randn distribution');
xlabel('Seed');
ylabel('randn value');
grid on;

בתוך לולאה שרצה על i (מ-1 עד 1250), אני מגריל מספר לפי התפלגות נורמלית שה-seed שלהן הוא i (ה-seed חייב להיות שלם, ואם הוא Matlab מעגל אותו). כך שלמעשה אני יכול לצייר את הפונקציה שמקשרת בין הערך המוגרל לבין ה-seed שיצר אותו.

חשוב לציין שאת הקוד הרצתי על Matlab R2008b. בינתיים אני רואה ש-Matlab שינו את המימוש ואת הקריאה לפונקציית randn. אין לי גישה עדיין לגרסת Matlab חדשה יותר, אך סביר להניח שתהיה לי בשבועות הקרובים, ואז אבדוק את זה שוב.
להפתעתי קיבלתי את הגרף הבא:

בגרף ציירתי את הערך שמתקבל מקריאה ל-randn כפונקציה של ה-seed. ניתן לראות שהפונקציה רחוקה מלהיות פסאודו-רנדומלית. למעשה היא כמעט דטרמיניסטית. ב-zoom על הגרף ניתן להבחין במבנים של ממש.

הרצתי את אותו סקריפט ב-octave, וקיבלתי את הגרף הבא:
ניתן לראות שהמימוש ב-Octave של אות אלגוריתם מניב תוצאות מפוזרות בהרבה.

ואיך כל זה מסביר את הבאג של החבר שלי? ברגע שקוראים ל-randn עם seed שמסתמך על ה-clock (שהרזולוציה שלו היא מילישנייה) בתוך לולאה, והביצוע של הלולאה קצר יחסית (יחסית להתקדמות של ה-clock), צפויים לקבל ערכי seed קרובים אחד לשני, ואו אז ההתפלגות ממנה יוגרלו הערכים לא תיראה נורמלית (כי למעשה נגריל ערכים מאיזור מסויים בגרף הראשון.

מסקנות:

1. המימוש ב-Octave עדיף על זה של Matlab.
2. אם משתמשים ב-Matlab, לא כדאי להסתמך על ה-clock ליצירת ה-seed בקריאה ל-randn, אלא אם כן אתם בטוחים שמשך זמן ביצוע לולאה בודדת ארוך יחסית.
   
3. כדאי לבדוק מה עשו ב-Matlab בגרסאות חדשות יותר, הם יצרו עוד גנרטורים חוץ מ-state, וייתכן שהספק-באג הזה בכלל תוקן.

(וטיפ: מי שמריץ את הסקריפט על Octave בגרסה 3.0 ומקבל הודעת שגיאה שקשורה לפונטים כשהוא מבקש לשמור את הגרף בעזרת פקודת print - כדאי לו לעיין בפוסט הזה. בגרסה 3.2 הבאג אמור להיות מתוקן, אך עוד לא בדקתי).