המטלה זכתה לציון של 99/100
. מובאת כאן לסיוע לצרכי למידה אישית לאחר פרסום הפתרון הרשמי, ולא כתחליף לו. אין להעתיק או להסתמך על הדברים.
1 - ריבועים פחותים ללא חותך
הוכחה שהשיפוע האופטימלי $a$ עבור רגרסיה לינארית ללא חותך הוא:
\[a = \frac{\sum_{i=1}^{n} x_iy_i}{\sum_{i=1}^{n} x_i^2}\]הערה: נניח כי $\sum_{i=1}^{n} x_i^2 \neq 0$.
1(א)
הנוסחה הרגילה שהוצגה בכתה לסכום הריבועים הפחותים הייתה:
\[E = \sum_{i=1}^{n} e_i^2 = \sum_{i=1}^{n} (y_i - \hat{y}_i)^2 = \sum_{i=1}^{n} (y_i - (ax_i + b))^2\]נציב $b=0$ ונקבל:
\[\sum_{i=1}^{n} (y_i - (ax_i + b))^2 = \boxed{\sum_{i=1}^{n} (y_i - ax_i)^2}\]כאשר:
- $e_i$ היא השארית (הטעות) בתצפית $i$
- $y_i$ הוא הערך המדוד של המשתנה התלוי
- $x_i$ הוא הערך המדוד של המשתנה הבלתי תלוי
- $a$ הוא השיפוע שאנחנו מנסים לאמוד
1(ב)
נניח עד לסעיף ה שנתונות לנו שתי נקודות נתונים $(x_1, y_1)$ ו-$(x_2, y_2)$.
נחשב את הביטוי המפורש של $E$ עבור $n = 2$:
\[\begin{aligned} E &= \sum_{i=1}^{2} (y_i - ax_i)^2 \\ &= (y_1 - ax_1)^2 + (y_2 - ax_2)^2 \\ &= (y_1^2 - 2ax_1y_1 + a^2x_1^2) + (y_2^2 - 2ax_2y_2 + a^2x_2^2) \\ &= \boxed{y_1^2 - 2ax_1y_1 + a^2x_1^2 + y_2^2 - 2ax_2y_2 + a^2x_2^2} \end{aligned}\]1(ג)
לאחר כינוס איברים דומים, נקבל:
\[E = \boxed{a^2(x_1^2 + x_2^2) - 2a(x_1y_1 + x_2y_2) + (y_1^2 + y_2^2)}\]1(ד)
הפונקציה $E(a) = pa^2 + qa + r$ היא פרבולה במשתנה $a$. המקדמים הם:
\[\begin{aligned} p &= x_1^2 + x_2^2 = \sum_{i=1}^{2} x_i^2 \\ q &= -2(x_1y_1 + x_2y_2) = -2\sum_{i=1}^{2} x_iy_i \\ r &= y_1^2 + y_2^2 = \sum_{i=1}^{2} y_i^2 \end{aligned}\]1(ה)
עבור פרבולה $px^2 + qx + r$, הנגזרת היא $2px + q$. הנקודה שבה הנגזרת מתאפסת (נקודת מינימום) היא כאשר $2pa + q = 0$, ולכן $a = -q/(2p)$.
מכאן שהערך של $a$ שממזער את $E(a)$ נמצא בנקודה $a = -q/(2p)$. נציב את המקדמים שמצאנו בסעיף ד ונקבל:
\[\begin{aligned} a &= -\frac{q}{2p} \\ &= \frac{-(-2\sum_{i=1}^{2} x_iy_i)}{2(x_1^2 + x_2^2)} \\ &= \frac{\sum_{i=1}^{2} x_iy_i}{x_1^2 + x_2^2} \\ &= \boxed{\frac{\sum_{i=1}^{2} x_iy_i}{\sum_{i=1}^{2} x_i^2}} \end{aligned}\]כאשר $x \neq 0$.
1(ו)
באופן כללי, עבור $n$ נקודות נתונים, השיפוע האופטימלי $a$ כאשר $b=0$ הוא:
\[a = \frac{\sum_{i=1}^{n} x_iy_i}{\sum_{i=1}^{n} x_i^2}\]אפשר להוכיח את זה על ידי אינדוקציה על $n$, להבנתנו זה מעבר לדרישות הקורס אבל בקצרה:
(1) בסיס האינדוקציה: $n=2$ הוכח בסעיף ד. (2) נניח שזה נכון עבור $n=k$. (3) נראה שזה נכון עבור $n=k+1$
2 - statsmodels
2(א)
ניתן להציג את הבעיה של רגרסיה לינארית עם ובלי חותך בעזרת מודל רגרסיה לינארית של statsmodles
. למשל עבור הדוגמה שהוצגה בכתה:
import numpy as np
import statsmodels.api as sm
# Data from the original problem
x = np.array([1, 2, 3, 4])
y = np.array([2.5, 4, 4, 7])
# With intercept (standard model)
X_with_intercept = sm.add_constant(x)
model_with_intercept = sm.OLS(y, X_with_intercept)
results_with_intercept = model_with_intercept.fit()
# Without intercept (as in exercise 1)
model_without_intercept = sm.OLS(y, x)
results_without_intercept = model_without_intercept.fit()
print(results_with_intercept.params) # [1. 1.35]
print(results_without_intercept.params) # [1.68333333]
2(ב)
מימוש של רגרסיה לינארית עם חותך:
# Manual calculation with formula from class
x_mean = np.mean(x) # = 2.5
y_mean = np.mean(y) # = 4.375
# Calculate a = Σ[(x_i - x̄)(y_i - ȳ)]/Σ[(x_i - x̄)²]
numerator = sum((x - x_mean) * (y - y_mean)) # = 6.75
denominator = sum((x - x_mean)**2) # = 5
a = numerator / denominator # = 1.35
# Calculate b = ȳ - a·x̄
b = y_mean - a * x_mean # = 4.375 - 1.35*2.5 = 1.0
מימוש של רגרסיה לינארית ללא חותך:
# Manual calculation with formula from exercise 1
numerator = sum(x * y) # = 50.5
denominator = sum(x**2) # = 30
a_without_intercept = numerator / denominator # = 1.68333
2(ג)
על בסיס החישובים והשוואה עם statsmodels
:
(1) כאשר משתמשים בפונקציה add_contant(x)
, statsmodels OLS
משתמש בנוסחה שהוצגה בכתה (עם חותך). הפלט: $a = 1.35$, $b = 1.0$ תואם לנוסחה:
(2) כאשר לא משתמשים ב-add_constant(x)
, statsmodels OLS
משתמש בנוסחה ללא חותך. הפלט: $a = 1.68333$ תואם לנוסחה:
2(ד)
מצאנו מספר פתרונות:
(1) שימוש במטודות המובנות של statsmodels
: להיעזר בפונקציה add_constant(x)
כדי להוסיף עמודת חותך במקרה שהמשתמש מעוניין במודל עם חותך. לחלופין, ניתן גם להשתמש במודל ישירות אם מעונינים בנוסחה ללא חותך:
model_without_intercept = sm.OLS(y, x)
results_without_intercept = model_without_intercept.fit()
(2) מניפולציות על על הנתונים: - ניתן להוסיף עמודת חותך למודל על ידי הוספת עמודת אחדות ל-X. - ניתן גם לבטל את החותך - אפילו אם המשתמש קרא לפונציקה add_constant(x)
, דרך איפוס עמודות האחדות ויצירת $b=0$ באופן מלאכותי:
X_zeros_intercept = sm.add_constant(x)
X_zeros_intercept[:, 0] = 0 # Replace the constant column with zeros
model_zero_intercept = sm.OLS(y, X_zeros_intercept)
results_zero_intercept = model_zero_intercept.fit()
שתי השיטות הדפיסו פלט זהה ($slope = 1.68333$) עבור רגרסיה לינארית ללא חותך, שגם תאם את החישוב הידני שלנו בעזרת הנוסחה מתרגיל 1.
לשלמות התמונה, להלן התיעוד הרשמי של הפונקציה הרלוונטית:
def add_constant(data, prepend=True, has_constant='skip'):
"""
Add a column of ones to an array.
Parameters
----------
data : array_like
A column-ordered design matrix.
prepend : bool
If true, the constant is in the first column. Else the constant is
appended (last column).
has_constant : str {'\text{raise}', 'add', 'skip'}
Behavior if ``data`` already has a constant. The default will return
data without adding another constant. If '\text{raise}', will \text{raise} an
error if any column has a constant value. Using 'add' will add a
column of 1s if a constant column is present.
Returns
-------
array_like
The original values with a constant (column of ones) as the first or
last column. Returned value type depends on input type.
Notes
-----
When the input is a pandas Series or DataFrame, the added column's name
is 'const'.
"""
if _is_using_pandas(data, None):
from statsmodels.tsa.tsatools import add_trend
return add_trend(data, trend='c', prepend=prepend, has_constant=has_constant)
# Special case for NumPy
x = np.asarray(data)
ndim = x.ndim
if ndim == 1:
x = x[:, None]
elif x.ndim > 2:
\text{raise} ValueError('Only implemented for 2-dimensional arrays')
is_nonzero_const = np.ptp(x, axis=0) == 0
is_nonzero_const &= np.all(x != 0.0, axis=0)
if is_nonzero_const.any():
if has_constant == 'skip':
return x
elif has_constant == '\text{raise}':
if ndim == 1:
\text{raise} ValueError("data is constant.")
else:
columns = np.arange(x.shape[1])
cols = ",".join([str(c) for c in columns[is_nonzero_const]])
\text{raise} ValueError(f"Column(s) {cols} are constant.")
x = [np.ones(x.shape[0]), x]
x = x if prepend else x[::-1]
return np.column_stack(x)
הבדיקות נעשו על הגסאות הבאות:
Statsmodels version: 0.14.4
NumPy version: 2.2.3
Pandas version: 2.2.3
SciPy version: 1.15.2
Matplotlib version: 3.10.0
3 - שקית תרופות
3(א)
הסתברות לשלוף כל משכך כאבים שקולה למאורע המשלים ״לא לשלוף אנטיביוטיקה. הסיכוי לשלוף אנטיביוטיקה הוא 5 מתוך 50 כדורים, כלומר:
\[P(\text{not antibiotic}) = 1 - P(\text{antibiotic}) = 1 - \frac{5}{50} = \frac{45}{50} = 0.9\]כלומר ההסתברות לשלוף משכך כאבים היא $\boxed{90\%}$.
3(ב)
מר קשישא מת אם הוא לוקח אופטלגין. הסתברות לשלוף אופטלגין ביום כלשהו = $\frac{10}{50} = 0.2$ הסתברות לא לשלוף אופטלגין ביום כלשהו = $1 - 0.2 = 0.8$ הסתברות להישרדות למשך 7 ימים לפי כלל המכפלה = $0.8^7 \approx 0.21$ או $\boxed{21\%}$
3(ג)
כאשר השקית אינה מתמלאת מחדש בשלושת הימים הראשונים, ההסתברויות משתנות בכל יום:
- יום 1: הסתברות לא לשלוף אופטלגין = $\frac{40}{50} = 0.8$
- יום 2: הסתברות לא לשלוף אופטלגין = $\frac{39}{49} \approx 0.796$ (נשארו 39 כדורים שאינם אופטלגין מתוך 49)
- יום 3: הסתברות לא לשלוף אופטלגין = $\frac{38}{48} \approx 0.792$ (נשארו 38 כדורים שאינם אופטלגין מתוך 48)
- ימים 4-7: הסתברות לא לשלוף אופטלגין = 0.8 לכל יום (שקית מתמלאת מחדש)
הסתברות להישרדות = $0.8 \times 0.796 \times 0.792 \times 0.8^4 \approx 0.206$ או $\boxed{20.6\%}$
3(ד)
לאחר לקיחת אנטיביוטיקה, מר קשישא לוקח כדורים כל 10 דקות למשך 30 דקות (מקסימום 3 כדורים). הוא מת אם הוא לוקח אופטלגין ומפסיק אם הוא לוקח פרצטמול או אספירין.
נחשב את ההסתברות שמר קשישא ישרוד עד הגעת האמבולנס באמצעות עץ הסתברויות.
מצב התחלתי של שקית התרופות:
- $20$ פרצטמול
- $15$ אספירין
- $10$ אופטלגין
- $4$ אנטיביוטיקה (כי הוא כבר לקח אחד בבוקר)
- סה”כ: $49$ כדורים
בניית עץ ההסתברויות המלא
כדור ראשון (10 דקות אחרי הקריאה לאמבולנס):
- לוקח פרצטמול: $\frac{20}{49} \approx 0.408$ - חי ומפסיק
- לוקח אספירין: $\frac{15}{49} \approx 0.306$ - חי ומפסיק
- לוקח אופטלגין: $\frac{10}{49} \approx 0.204$ - מת מיד
- לוקח אנטיביוטיקה: $\frac{4}{49} \approx 0.082$ - ממשיך לכדור שני
כדור שני (20 דקות אחרי הקריאה לאמבולנס):
אם לקח אנטיביוטיקה ראשונה, בשקית נשארו $48$ כדורים:
- לוקח פרצטמול: $0.082 \times \frac{20}{48} = 0.082 \times 0.417 \approx 0.0342$ - חי ומפסיק
- לוקח אספירין: $0.082 \times \frac{15}{48} = 0.082 \times 0.313 \approx 0.0257$ - חי ומפסיק
- לוקח אופטלגין: $0.082 \times \frac{10}{48} = 0.082 \times 0.208 \approx 0.0171$ - מת מיד
- לוקח אנטיביוטיקה: $0.082 \times \frac{3}{48} = 0.082 \times 0.063 \approx 0.0052$ - ממשיך לכדור שלישי
כדור שלישי (30 דקות אחרי הקריאה לאמבולנס):
אם לקח אנטיביוטיקה פעמיים, בשקית נשארו $47$ כדורים:
- לוקח פרצטמול: $0.0052 \times \frac{20}{47} = 0.0052 \times 0.426 \approx 0.0022$ - חי
- לוקח אספירין: $0.0052 \times \frac{15}{47} = 0.0052 \times 0.319 \approx 0.0017$ - חי
- לוקח אופטלגין: $0.0052 \times \frac{10}{47} = 0.0052 \times 0.213 \approx 0.0011$ - מת מיד
- לוקח אנטיביוטיקה: $0.0052 \times \frac{2}{47} = 0.0052 \times 0.043 \approx 0.0002$ - חי (האמבולנס מגיע)
סיכום כל מסלולי השרידות מהעץ:
\[\begin{aligned} P(\text{survive}) &= P(\text{paracetamol first}) \\ &+ P(\text{aspirin firest}) \\ &+ P(\text{antibiotic first, paracetamol second}) \\ &+ P(\text{antibiotic first, aspirin second}) \\ &+ P(\text{antibiotic first and twice, paracetamol third}) \\ &+ P(\text{antibiotic first and twice, aspirin third}) \\ &+ P(\text{antibiotics 3 times}) \end{aligned}\] \[P(\text{survive}) = 0.408 + 0.306 + 0.0342 + 0.0257 + 0.0022 + 0.0017 + 0.0002\] \[P(\text{survive}) = 0.778 \approx 77.8\%\]אפשר לאמת את התוצאה על ידי חישוב הסתברות המוות (לקיחת אופטלגין בשלב כלשהו):
\[P(\text{death}) = P(\text{optalgin first}) + P(\text{antibiotic first than optalgin}) + P(\text{antibiotic first, antibiotic second, optalgin third})\] \[P(\text{death}) = 0.204 + 0.0171 + 0.0011 = 0.222 \approx 22.2\%\]מכיוון ש-$P(\text{survive}) + P(\text{death}) = 0.778 + 0.222 = 1$, החישוב שלנו מאומת.
לכן, למר קשישא יש הסתברות של $77.8\%$ לשרוד עד הגעת האמבולנס.
3(ה)
אם נניח שרק אספירין מטפל בקריש הדם (ולא פרצטמול) והאמבולנס אינו מטפל בקריש, אז רק מי שלקח אספירין לפני הגעת האמבולנס ישרוד:
הסתברות ללקיחת אספירין לפני הגעת האמבולנס:
\[P(\text{took aspirin}) = 0.306 + 0.0257 + 0.0017 = 0.3334\]ולכן:
\[P(\text{survive one hour} \mid \text{hospital}) = \frac{0.3334}{0.778} \approx 0.429\]תחת הנחה זו, הסיכוי של מר קשישא לשרוד שעה לאחר תחילת הכאב יהיה כ-$42.9\%$.
3(ו)
עם התנאים החדשים:
- נכדו נותן לו אספירין כל בוקר
- השקית מכילה רק 20 פרצטמול ו-5 אנטיביוטיקה (סה”כ 25)
- לוקח כדור בהסתברות 0.9
- הכבד כושל לאחר 5 מנות פרצטמול רצופות
הסתברות ללקיחת פרצטמול = $0.9 \times \frac{20}{25} = 0.9 \times 0.8 = 0.72$ הסתברות לא לקיחת פרצטמול = $1 - 0.72 = 0.28$
הסתברות שהכבד לא כושל לאחר 5 ימים = $1 - 0.72^5 \approx 0.81$ או $\boxed{81\%}$
3(ז)
כשל כבד יקרה אם מר קשישא לוקח פרצטמול (פרצטמול) 5 ימים ברציפות. בתקופה של 6 ימים, רצף כזה יכול להתחיל ביום ה-1 או ביום ה-2.
נגדיר:
- A = המאורע שבימים 1-5 הוא לוקח פרצטמול
- B = המאורע שבימים 2-6 הוא לוקח פרצטמול
הסתברויות:
- $P(A) = 0.72^5 ≈ 0.1935$
- $P(B) = 0.72^5 ≈ 0.1935$
- $P(A ∩ B) = 0.72^6 ≈ 0.1393$ (המקרה שהוא לוקח פרצטמול בכל ששת הימים)
לפי עקרון ההכלה וההפרדה:
\[P(A ∪ B) = P(A) + P(B) - P(A ∩ B)\] \[P(A ∪ B) = 0.1935 + 0.1935 - 0.1393 = 0.2477\]זוהי ההסתברות שהכבד ייכשל. לכן, ההסתברות שהכבד לא ייכשל היא:
\[P(\text{liver does not fail}) = 1 - 0.2477 = 0.7523 ≈ \boxed{75.2\%}\]4 - strike
4(א)
אם שביתות והעלאות שכר הן בלתי תלויות:
\[P(\text{strike} \cap \text{raise}) = P(\text{raise}) \times P(\text{strike})\] \[\frac{1}{4} = r \times q\]מכאן:
\[\boxed{r = \frac{1}{4q}}\] \[\boxed{q = \frac{1}{4r}}\]4(ב)
מכיוון ש-$r$ ו-$q$ הן הסתברויות (בין 0 ל-1) ו-$r \times q = \frac{1}{4}$:
- ערך מינימלי עבור $r$: $\frac{1}{4}$ (כאשר $q = 1$)
- ערך מינימלי עבור $q$: $\frac{1}{4}$ (כאשר $r = 1$)
סך הכל:
\[\boxed{r \geq \frac{1}{4} \text{ and } q \geq \frac{1}{4}}\]4(ג)
\[P(\text{strike} | \text{raise}) = \frac{P(\text{strike} \cap \text{raise})}{P(\text{raise})} = \frac{\frac{1}{4}}{r} = \boxed{\frac{1}{4r}}\]הסבר: זוהי נוסחת ההסתברות המותנית. אנו מחשבים את ההסתברות לstrike בהינתן שהייתה $\text{raise}$.
4(ד)
\[P(\text{raise} | \text{not strike}) = \frac{P(\text{raise} \cap \text{not strike})}{P(\text{not strike})}\]ההסתברות ל$\text{raise}$ ו$\text{not strike}$ היא $P(\text{raise}) - P(\text{raise} \cap \text{strike}) = r - \frac{1}{4}$
ההסתברות ל$\text{not strike}$ היא $1 - P(\text{strike}) = 1 - q$
לכן:
\[P(\text{raise} | \text{not strike}) = \frac{r - \frac{1}{4}}{1 - q}\] דור פסקל