assignment14_montecarlo.ipynb

Kernel: Python 2 (system-wide)

שמות:

from IPython.display import HTML
import matplotlib.pyplot as plt

חישוב שטחים בשיטת מונטה קרלו

בתרגיל זה ניישם שיטה חישובית המכונה מונטה קרלו. בעזרת שיטה זו נחשב בהמשך שטח של אליפסה. אבל בשלב ראשון ניישם את השיטה המתוארת בהמשך לחישוב $\pi$.

###

קראו בעיון את הקטע שבהמשך

HTML('''<iframe src="https://he.wikipedia.org/wiki/%D7%A9%D7%99%D7%98%D7%AA_%D7%9E%D7%95%D7%A0%D7%98%D7%94_%D7%A7%D7%A8%D7%9C%D7%95"
     width=900 height=700"></iframe>''')

בהמשך קטע קוד המיישם את האלגוריתם המתואר למעלה למציאת $\pi$. הוסיפו ליד כל סימן סולמית הערה המסבירה את תפקיד משפט הקוד שמשמאלה.

from random import uniform #הבאת פונקציה שתביא מספר רנדומאלי לא שלם בתחום שנכתב
N=10000#מספר הנסיונות\נקודות
count=0#מספר הקליעות
for i in range(N):
    x=uniform(-1,1)#הגרלה של שיעור ה-איקס בין 1- ל1
    y=uniform(-1,1)#הגרלה של שיעור הוואי בין 1- ל1
    if x**2+y**2<=1:# בדיקה אם הוא נמצא בתוך או על המעגל
        count=count+1
float(count)/N*4

3.1408

חישוב $\pi$ בשיטה המתוארת כאן מכונה חישוב בדגימה ישירה

הריצו את הפונקציה תוך כדי הגדלת ערכי N מ-10 ועד ל- 10000. מה המסקנה? הסבירו במה במה תלוי הדיוק בחישוב $\pi$?

from random import uniform
def pi_Evaluation(N):
    count=0
    for i in range(N):
        x=uniform(0,1)
        y=uniform(0,1)
        if x**2+y**2<=1:
            count=count+1
    return float(count)/N*4
pi_Evaluation(100000)
#הדיוק בחישוב פאי תלוי במספר ההרצות

3.13716

מידת הפיזור של תוצאות "ניסוי" מכונה סטיית התקן בסידרת התרגילים שבהמשך תכתבו פונקציות שמטרתן לחשב את סטיית התקן המתקבלת בחישוב $\pi$ בדגימה ישירה. בשלב ראשון עליכם לכתוב פונקציה המריצה עבור אותו N את הפונקציה pi_Evaluation מספר פעמים המועבר אליה כפרמטר. הפונקציה צריכה להחזיר רשימה של כל ערכי ה-$\pi$ שהתקבלו.

השלימו את קטע הקוד שבהמשך.

<

סטיית תקן

def get_trails(k,N):
    pi_list=[pi_Evaluation(N) for i in range(k)]
    return pi_list
get_trails(5,10000)

[3.1468, 3.1424, 3.13, 3.12, 3.1432]

שורת הפקודות שלמטה משרטט תרשים שכיחויות של תוצאות של הרצות שונות עבור אותו גודל של N.הריצו את הקוד עבור ערכי N שונים ובחנו כיצד משתנה פיזור התוצאות.

import matplotlib.pyplot as plt
N=5000
trails=10000
o=plt.hist(get_trails(trails,N),bins=30)

העזרו בחבילה אותה פיתחתם לחישוב פונקציות סטיסטיות כדי לענות על השאלות שבהמשך.

בעזרת הפונקציה המחשבת ממוצע חשבו ערך ממוצע של $\pi$ המתקבל לאחר ביצוע של k ניסויים בהם בוחרים N נקודות.

import py_compile
py_compile.compile('my_statistics.py')

from my_statistics import mean
def mean_for_pi(k,N):
    return mean(get_trails(k,N))
mean_for_pi(5,10000)

3.14032

הקישור שלמטה הוא קישור לדף ווקיפדיה של המונח סטיית תקן. קראו על המונח . העזרו בחבילה אותה פיתחתם כדי לשרטט גרף של סטיית התקן כתלות ב- N. שנו את N מ- 10 ל-100 במרווחים של 10. קבעו את k ל-1000. השתמשו בפונקציה זו כדי לשרטט גרף של סטיית התקן כתלות ב-N. מה המסקנה?

סטיית תקן

from my_statistics import std
import matplotlib.pyplot as plt
from random import uniform
k=1000
stdL=[]
NL=[]
for N in range(10,110,10):
    NL.append(N)
    stdL.append(std(get_trails(k,N)))
plt.plot(NL,stdL,c='g',lw='3')
plt.xlabel('N',fontsize=16, color='red')
plt.ylabel('std',fontsize=16, color='red')
plt.title('$std by N$',fontsize=24,color='r')
#ככל שיש יותר מקרים, סטיית התקן הולכת ויורדת

שרטטו גרף של סטיית התקן כתלות ב: $\sqrt{N}\over{N}$. האם קיבלתם בקרוב יחס ישר? אם כן כיצד סטיית התקן תלויה ב-N?

from my_statistics import std
import matplotlib.pyplot as plt
from random import uniform
k=1000
stdL=[]
NL=[]
for N in range(10,110,10):
    NL.append(float(N**(-0.5)))
    stdL.append(std(get_trails(k,N)))
plt.plot(NL,stdL,c='g',lw='3')
plt.xlabel('N**(-1/2)',fontsize=16, color='red')
plt.ylabel('std',fontsize=16, color='red')
plt.title('$std by N**(-1/2)$',fontsize=24,color='r')
#היחס בין שניהם מתנהג כמו יחס ישר

Text(0.5,1,u'$std by N**(-1/2)$')

אליפסה היא המקום הגיאומטרי של כל הנקודות שסכום מרחקיהן משתי נקודות קבועות המכונות מוקדי האליפסה הוא קבוע ושווה לקוטר הגדול של האליפסה, 2a, הקוטר הקטן של האליפסה שווה ל-2b.
כתבו פונקציה בשם ellipse_area המקבל שלושה פרמטרים a, b ו- tol המחשבת שטח של אליפסה שחצי הקוטר הגדול שלה הוא a וחצי הקוטר הקטן הוא b בדיוק של tol. השגיאה בחישוב השטח היא מסדר גודל של סטיית התקן.

def monte_carlo_for_ellipse(a,b,N):
    rec_area = 4*a*b
    count=0
    for i in range(N):
        x=uniform(-a,a)
        y=uniform(-b,b)
        if ((x*x)/(a*a) + (y*y)/(b*b) <= 1):
            count=count+1
    ellipse_area2=(float(count)/N)*rec_area
    return ellipse_area2

import random
from math import pi
def ellipse_area(a,b,tol):
    N=10
    ellipse_area1 = pi*a*b#חישוב לפי נוסחא(חישוב )
    ellipse_area2 = monte_carlo_for_ellipse(a,b,N)#חישוב לפי שיטת מונטה קרלו
    while (abs(ellipse_area2 - ellipse_area1) > tol):#בדיוק של לפחות 0.3
        N*=10
        ellipse_area2 = monte_carlo_for_ellipse(a,b,N)
    return ellipse_area1, ellipse_area2#השוואה בין השטחים
ellipse_area(6,3,0.3)

(56.548667764616276, 56.4192)