מספרים אקראיים ב-R: rnorm, runif, sample ו-set.seed

set.seed() קודם לכול

הגרלות שמשתנות בכל הרצה חסרות תועלת להוראה, לבדיקה, לדיבוג או למדע: אם הסימולציה שלכם אומרת היום 0.146 ומחר 0.153, איזה מספר נכנס לדוח? set.seed() מקבעת את נקודת ההתחלה של המחולל, וכך כל זרם ה"אקראיות" שבא אחריה ניתן לשחזור מדויק:

set.seed(42)
print(rnorm(3))

set.seed(42)          # same seed...
print(rnorm(3))       # ...same "random" numbers

set.seed(7)           # different seed, different stream
print(rnorm(3))

אותו seed, הגרלות זהות, בכל פעם, בכל מחשב. המספרים הם פסאודו-אקראיים: סדרה דטרמיניסטית שתוכננה לעבור כל מבחן סטטיסטי של אקראיות, וה-seed בוחר מאיפה בסדרה מתחילים. לערך ה-seed עצמו אין משמעות: 42, 7, 20260807, בחרו מה שתרצו; רק תעדו אותו. ההרגל שכדאי לבנות: set.seed() אחד בראש כל סקריפט שמשתמש באקראיות. (שימו לב שכל הגרלה מקדמת את המצב, ולכן גם סדר הקריאות משנה לשחזור.)

שיטת d/p/q/r: טבלה אחת שמפענחת את כל הספרייה

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

p ו-q הופכיים זה לזה, והם הזוג שפגשתם ב-רווחי סמך בתור qt(0.975, df):

print(pnorm(1.96))      # P(Z <= 1.96), about 0.975
print(qnorm(0.975))     # which z has 97.5% below it? About 1.96
print(dnorm(0))         # height of the bell curve at its peak

החליפו את שם המשפחה והכול עובר איתו: runif/dunif/punif/qunif, rbinom/..., rpois/..., rt/..., rexp/.... לומדים ארבע קידומות, מקבלים עשרות התפלגויות.

rnorm(): הגרלות מהתפלגות נורמלית

rnorm(n, mean, sd) מגרילה מעקומת פעמון, סוס העבודה לסימולציה של נתונים שדומים למדידות:

set.seed(42)
iq <- rnorm(1000, mean = 100, sd = 15)

print(mean(iq))
print(sd(iq))
print(mean(iq > 130))   # share more than 2 sd above the mean

הממוצע וסטיית התקן של המדגם נוחתים ליד 100 ו-15, לא בדיוק עליהם: הפער הזה הוא רעש הדגימה, והוא מצטמצם ככל ש-n גדל. השורה האחרונה היא טריק ששווה לאמץ: mean() של וקטור לוגי הוא שיעור ה-TRUE, והחלק שמעל 130 יוצא קרוב לערך התיאורטי 1 - pnorm(130, 100, 15) ≈ 2.3%. ברירות המחדל הן mean = 0, sd = 1 (ההתפלגות הנורמלית הסטנדרטית). היסטוגרמה של iq מראה את הפעמון המוכר.

runif(), rbinom(), rpois()

עוד שלוש משפחות מכסות את רוב צורכי הסימולציה:

set.seed(42)

print(runif(5))                    # uniform on [0, 1]
print(runif(3, min = 10, max = 20))

print(rbinom(8, size = 10, prob = 0.5))   # successes in 10 coin flips, 8 times

print(rpois(8, lambda = 3))        # event counts averaging 3 per interval

runif() מפזרת הגרלות באופן אחיד על פני טווח ("uniform", לא "run if": כולם קוראים את זה לא נכון פעם אחת). rbinom(n, size, prob) מדמה n ניסויים של size חזרות כל אחד, ומחזירה את מספר ההצלחות בכל אחד, כך שכל ערך למעלה הוא מספר ה"עץ" מתוך 10 הטלות. rpois(n, lambda) מייצרת ספירות של אירועים שקורים באופן בלתי תלוי בקצב ממוצע ידוע: פניות תמיכה לשעה, שגיאות הקלדה לעמוד.

sample(): דגימה וערבוב

בעוד שפונקציות r* ממציאות ערכים מתוך התפלגות, sample() מגרילה מתוך ערכים שכבר יש לכם:

set.seed(42)

print(sample(1:10, 3))                              # 3 values, no repeats
print(sample(c("win", "lose"), 5, replace = TRUE))  # 5 draws, repeats allowed
print(sample(1:5))                                  # no size = shuffle everything

הארגומנט replace הוא כל הסיפור: FALSE (ברירת המחדל) מחלק קלפים: כל ערך יכול להופיע פעם אחת, ובקשה ליותר ממה שיש היא שגיאה; TRUE מטיל קוביות: כל הגרלה מתחילה מחדש. דגימה עם החזרה מהנתונים שלכם היא המנוע של bootstrapping. כשקוראים לה רק עם וקטור, sample(x) מחזירה תמורה אקראית: הדרך המקובלת לערבב.

דגימת שורות של data frame משתמשת ב-sample() בתוך בחירת שורות לפי אינדקס:

set.seed(42)

random_cars <- mtcars[sample(nrow(mtcars), 5), ]
print(random_cars[, c("mpg", "wt", "hp")])

sample(nrow(mtcars), 5) בוחרת 5 מספרי שורות אקראיים; בחירת האינדקס שולפת את השורות האלה. זה הצעד המקובל לבדיקה מדגמית של מאגר נתונים גדול או לחלוקה לקבוצות אימון ובדיקה.

סימולציית מונטה קרלו קטנה

הנה השכר של כל ארגז הכלים. השאלה: תהליך מפיק מדידות שמתפלגות N(100, 15); ממצעים 10 מהן. מה ההסתברות שהממוצע יעלה על 105? במקום לגזור את התשובה, עשו סימולציה: בצעו את הניסוי 10,000 פעמים וספרו:

set.seed(42)

trial <- function() {
  mean(rnorm(10, mean = 100, sd = 15)) > 105
}

results <- replicate(10000, trial())
print(mean(results))

# The exact answer, for comparison:
print(1 - pnorm(105, mean = 100, sd = 15 / sqrt(10)))

הסימולציה נוחתת במרחק של אלפיות ספורות מהערך המדויק (כ-0.146). זה מונטה קרלו בנשימה אחת: כתבו ניסוי אחד בתור פונקציה, הריצו אותו אלפי פעמים עם replicate(), וקחו mean() של ההצלחות. התשובה המדויקת הייתה קיימת כאן כי הבעיה הייתה פשוטה כמו בספר לימוד; ברגע שהשאלה מסתבכת (התפלגויות מוזרות, מקסימום של הגרלות מתואמות, משחק מבוסס כללים), הדרך האנליטית נסגרת, ומתכון הסימולציה ממשיך לעבוד בלי שינוי. שימו לב ל-15 / sqrt(10) בבדיקה: ממוצעים מתנדנדים פחות מהגרלות בודדות, אותו חוק שורש ריבועי שקובע את רוחב רווחי הסמך.

ברמת סימולציה, לא ברמת סודות

גבול אחד שצריך לכבד: מחולל ברירת המחדל של R (Mersenne Twister) בנוי לאיכות סטטיסטית ולמהירות, לא לסודיות. הפלט שלו דטרמיניסטי בהינתן ה-seed, ואפשר לשחזר את המצב הפנימי שלו מתוך פלט שנצפה: פגמים קטלניים לסיסמאות, לאסימוני אבטחה או לכל דבר שקשור לאבטחה. לסימולציה, ל-bootstrapping ולהוראה הוא מצוין; לקריפטוגרפיה השתמשו בספרייה שנבנתה לשם כך (למשל החבילה openssl), אף פעם לא ב-sample() או ב-runif().

מה לקחת מכאן

• set.seed() אחד בראש הסקריפט הופך כל תוצאה "אקראית" לניתנת לשחזור: אותו seed, אותן הגרלות.
• הקידומות d/p/q/r מפענחות את כל ספריית ההתפלגויות: הגרלה אקראית, צפיפות, הסתברות מצטברת, קוונטיל.
• rnorm(n, mean, sd), runif(n, min, max), rbinom(n, size, prob), rpois(n, lambda) מכסות את רוב צורכי הסימולציה.
• sample() מגרילה מהערכים שלכם: replace = TRUE לקוביות, ברירת המחדל לקלפים, בלי size כדי לערבב; df[sample(nrow(df), k), ] דוגמת שורות.
• מונטה קרלו = פונקציה של ניסוי אחד + replicate() + mean(): המתכון שעונה על שאלות הסתברות שהמתמטיקה לא מגיעה אליהן בנוחות.
• מחולל המספרים האקראיים של R מיועד לסימולציה, לא לקריפטוגרפיה.

הבא בתור: דיבוג. מה הודעות השגיאה של R אומרות באמת ואיך קוראים traceback.