10. NumPy#
“بیایید صریح باشیم: کار علم هیچ ربطی به اجماع ندارد. اجماع کار سیاست است. برعکس، علم فقط به یک محقق نیاز دارد که اتفاقاً درست بگوید، که به این معنی است که او نتایجی دارد که با ارجاع به دنیای واقعی قابل تأیید هستند. در علم، اجماع بیربط است. آنچه مرتبط است، نتایج تکرارپذیر است.” – مایکل کرایتون
علاوه بر آنچه در Anaconda موجود است، این درس به کتابخانههای زیر نیاز دارد:
!pip install quantecon
10.1. مروری کلی#
NumPy یک کتابخانه درجه یک برای برنامهنویسی عددی است
به طور گسترده در دانشگاهها، امور مالی و صنعت استفاده میشود.
بالغ، سریع، پایدار و تحت توسعه مستمر است.
ما قبلاً در درسهای قبلی کدهایی شامل NumPy دیدهایم.
در این درس، بحث سیستماتیکتری را در مورد موارد زیر آغاز خواهیم کرد:
آرایههای NumPy و
عملیات پردازش آرایه اساسی که توسط NumPy ارائه میشوند.
(برای یک مرجع جایگزین، به مستندات رسمی NumPy مراجعه کنید.)
ما از import های زیر استفاده خواهیم کرد.
import numpy as np
import random
import quantecon as qe
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
from matplotlib import cm
10.2. آرایههای NumPy#
مشکل اساسی که NumPy حل میکند، پردازش سریع آرایه است.
مهمترین ساختاری که NumPy تعریف میکند، نوع داده آرایه است که به صورت رسمی numpy.ndarray نامیده میشود.
آرایههای NumPy بخش بسیار بزرگی از اکوسیستم علمی Python را پشتیبانی میکنند.
10.2.1. مبانی#
برای ایجاد یک آرایه NumPy که فقط شامل صفر است از np.zeros استفاده میکنیم
a = np.zeros(3)
a
array([0., 0., 0.])
type(a)
numpy.ndarray
آرایههای NumPy تا حدی شبیه لیستهای بومی Python هستند، به جز اینکه
دادهها باید همگن باشند (همه عناصر از یک نوع).
این انواع باید یکی از انواع داده (
dtypes) ارائه شده توسط NumPy باشند.
مهمترین این dtypes عبارتند از:
float64: عدد ممیز شناور 64 بیتی
int64: عدد صحیح 64 بیتی
bool: 8 بیتی True یا False
همچنین dtypes هایی برای نمایش اعداد مختلط، اعداد صحیح بدون علامت و غیره وجود دارد.
در ماشینهای مدرن، dtype پیشفرض برای آرایهها float64 است
a = np.zeros(3)
type(a[0])
numpy.float64
اگر بخواهیم از اعداد صحیح استفاده کنیم، میتوانیم به صورت زیر مشخص کنیم:
a = np.zeros(3, dtype=int)
type(a[0])
numpy.int64
10.2.2. شکل و بعد#
انتساب زیر را در نظر بگیرید
z = np.zeros(10)
در اینجا z یک آرایه مسطح است — نه بردار سطر و نه ستون.
z.shape
(10,)
در اینجا tuple شکل فقط یک عنصر دارد که طول آرایه است (tuple های با یک عنصر با کاما پایان مییابند).
برای اضافه کردن یک بعد اضافی، میتوانیم ویژگی shape را تغییر دهیم
z.shape = (10, 1) # تبدیل آرایه مسطح به بردار ستونی (دو بعدی)
z
array([[0.],
[0.],
[0.],
[0.],
[0.],
[0.],
[0.],
[0.],
[0.],
[0.]])
z = np.zeros(4) # آرایه مسطح
z.shape = (2, 2) # آرایه دو بعدی
z
array([[0., 0.],
[0., 0.]])
در مورد آخر، برای ساخت آرایه 2×2، میتوانیم یک tuple را نیز به تابع zeros() ارسال کنیم، مانند
z = np.zeros((2, 2)).
10.2.3. ایجاد آرایهها#
همانطور که دیدهایم، تابع np.zeros یک آرایه از صفرها ایجاد میکند.
احتمالاً میتوانید حدس بزنید که np.ones چه چیزی ایجاد میکند.
مرتبط با آن np.empty است که آرایههایی را در حافظه ایجاد میکند که بعداً میتوان آنها را با داده پر کرد
z = np.empty(3)
z
array([0., 0., 0.])
اعدادی که در اینجا میبینید مقادیر زباله هستند.
(Python سه قطعه متوالی 64 بیتی حافظه را اختصاص میدهد و محتویات موجود آن slot های حافظه به عنوان مقادیر float64 تفسیر میشوند)
برای راهاندازی یک شبکه از اعداد با فاصله یکسان از np.linspace استفاده کنید
z = np.linspace(2, 4, 5) # از 2 تا 4، با 5 عنصر
برای ایجاد یک ماتریس همانی از np.identity یا np.eye استفاده کنید
z = np.identity(2)
z
array([[1., 0.],
[0., 1.]])
علاوه بر این، آرایههای NumPy را میتوان از لیستهای Python، tuple ها و غیره با استفاده از np.array ایجاد کرد
z = np.array([10, 20]) # ndarray از لیست Python
z
array([10, 20])
type(z)
numpy.ndarray
z = np.array((10, 20), dtype=float) # در اینجا 'float' معادل 'np.float64' است
z
array([10., 20.])
z = np.array([[1, 2], [3, 4]]) # آرایه 2 بعدی از یک لیست از لیستها
z
array([[1, 2],
[3, 4]])
همچنین np.asarray را ببینید که عملکرد مشابهی انجام میدهد، اما یک نسخه مجزا از دادههای موجود در یک آرایه NumPy ایجاد نمیکند.
برای خواندن دادههای آرایه از یک فایل متنی حاوی دادههای عددی، از np.loadtxt استفاده کنید — برای جزئیات به مستندات مراجعه کنید.
10.2.4. نمایهگذاری آرایه#
برای یک آرایه مسطح، نمایهگذاری مانند توالیهای Python است:
z = np.linspace(1, 2, 5)
z
array([1. , 1.25, 1.5 , 1.75, 2. ])
z[0]
np.float64(1.0)
z[0:2] # دو عنصر، از عنصر 0 شروع میشود
array([1. , 1.25])
z[-1]
np.float64(2.0)
برای آرایههای دو بعدی، نحو نمایه به صورت زیر است:
z = np.array([[1, 2], [3, 4]])
z
array([[1, 2],
[3, 4]])
z[0, 0]
np.int64(1)
z[0, 1]
np.int64(2)
و الی آخر.
ستونها و سطرها را میتوان به صورت زیر استخراج کرد
z[0, :]
array([1, 2])
z[:, 1]
array([2, 4])
آرایههای NumPy از اعداد صحیح نیز میتوانند برای استخراج عناصر استفاده شوند
z = np.linspace(2, 4, 5)
z
array([2. , 2.5, 3. , 3.5, 4. ])
indices = np.array((0, 2, 3))
z[indices]
array([2. , 3. , 3.5])
در نهایت، یک آرایه از dtype bool میتواند برای استخراج عناصر استفاده شود
z
array([2. , 2.5, 3. , 3.5, 4. ])
d = np.array([0, 1, 1, 0, 0], dtype=bool)
d
array([False, True, True, False, False])
z[d]
array([2.5, 3. ])
در زیر خواهیم دید که چرا این مفید است.
یک نکته جانبی: همه عناصر یک آرایه را میتوان با استفاده از نماد slice برابر با یک عدد قرار داد
z = np.empty(3)
z
array([2. , 3. , 3.5])
z[:] = 42
z
array([42., 42., 42.])
10.2.5. متدهای آرایه#
آرایهها متدهای مفیدی دارند که همگی به دقت بهینه شدهاند
a = np.array((4, 3, 2, 1))
a
array([4, 3, 2, 1])
a.sort() # a را در محل مرتب میکند
a
array([1, 2, 3, 4])
a.sum() # مجموع
np.int64(10)
a.mean() # میانگین
np.float64(2.5)
a.max() # بیشینه
np.int64(4)
a.argmax() # ایندکس عنصر حداکثر را برمیگرداند
np.int64(3)
a.cumsum() # مجموع تجمعی عناصر a
array([ 1, 3, 6, 10])
a.cumprod() # حاصلضرب تجمعی عناصر a
array([ 1, 2, 6, 24])
a.var() # واریانس
np.float64(1.25)
a.std() # انحراف معیار
np.float64(1.118033988749895)
a.shape = (2, 2)
a.T # معادل a.transpose()
array([[1, 3],
[2, 4]])
متد دیگری که ارزش دانستن دارد searchsorted() است.
اگر z یک آرایه غیر نزولی باشد، پس z.searchsorted(a) ایندکس
اولین عنصر z که >= a است را برمیگرداند
z = np.linspace(2, 4, 5)
z
array([2. , 2.5, 3. , 3.5, 4. ])
z.searchsorted(2.2)
np.int64(1)
10.3. عملیات حسابی#
عملگرهای +، -، *، / و ** همگی به صورت عنصر به عنصر روی آرایهها عمل میکنند
a = np.array([1, 2, 3, 4])
b = np.array([5, 6, 7, 8])
a + b
array([ 6, 8, 10, 12])
a * b
array([ 5, 12, 21, 32])
میتوانیم یک اسکالر را به هر عنصر به صورت زیر اضافه کنیم
a + 10
array([11, 12, 13, 14])
ضرب اسکالر مشابه است
a * 10
array([10, 20, 30, 40])
آرایههای دوبعدی از همان قوانین کلی پیروی میکنند
A = np.ones((2, 2))
B = np.ones((2, 2))
A + B
array([[2., 2.],
[2., 2.]])
A + 10
array([[11., 11.],
[11., 11.]])
A * B
array([[1., 1.],
[1., 1.]])
به ویژه، A * B حاصلضرب ماتریسی نیست، بلکه یک حاصلضرب عنصر به عنصر است.
10.4. ضرب ماتریسی#
ما از نماد @ برای ضرب ماتریسی استفاده میکنیم، به صورت زیر:
A = np.ones((2, 2))
B = np.ones((2, 2))
A @ B
array([[2., 2.],
[2., 2.]])
نحو با آرایههای مسطح کار میکند — NumPy حدس آموزشدیدهای از آنچه شما میخواهید میزند:
A @ (0, 1)
array([1., 1.])
از آنجا که ما در حال ضرب سمت راست هستیم، tuple به عنوان یک بردار ستونی تلقی میشود.
10.5. Broadcasting#
(این بخش یک بحث عالی در مورد broadcasting ارائه شده توسط Jake VanderPlas را گسترش میدهد.)
Note
Broadcasting یک جنبه بسیار مهم از NumPy است. در عین حال، broadcasting پیشرفته نسبتاً پیچیده است و برخی از جزئیات زیر را میتوان در اولین مرور اجمالی خواند.
در عملیات عنصر به عنصر، آرایهها ممکن است شکل یکسانی نداشته باشند.
وقتی این اتفاق میافتد، NumPy به طور خودکار آرایهها را به شکل یکسان گسترش میدهد هرگاه امکانپذیر باشد.
این ویژگی مفید (اما گاهی گیجکننده) در NumPy broadcasting نامیده میشود.
ارزش broadcasting این است که
حلقههای
forمیتوانند اجتناب شوند، که به کد عددی کمک میکند تا سریع اجرا شود وbroadcasting میتواند به ما اجازه دهد عملیات را روی آرایهها پیادهسازی کنیم بدون اینکه واقعاً برخی ابعاد این آرایهها را در حافظه ایجاد کنیم، که میتواند مهم باشد وقتی آرایهها بزرگ هستند.
به عنوان مثال، فرض کنید a یک آرایه \(3 \times 3\) است (a -> (3, 3))، در حالی که b یک آرایه مسطح با سه عنصر است (b -> (3,)).
هنگام جمع کردن آنها با هم، NumPy به طور خودکار b -> (3,) را به b -> (3, 3) گسترش میدهد.
جمع عنصر به عنصر منجر به یک آرایه \(3 \times 3\) میشود
a = np.array(
[[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
b = np.array([3, 6, 9])
a + b
array([[ 4, 8, 12],
[ 7, 11, 15],
[10, 14, 18]])
در اینجا یک نمایش بصری از این عملیات broadcasting آورده شده است:
در مورد b -> (3, 1) چطور؟
در این حالت، NumPy به طور خودکار b -> (3, 1) را به b -> (3, 3) گسترش میدهد.
جمع عنصر به عنصر سپس منجر به یک ماتریس \(3 \times 3\) میشود
b.shape = (3, 1)
a + b
array([[ 4, 5, 6],
[10, 11, 12],
[16, 17, 18]])
در اینجا یک نمایش بصری از این عملیات broadcasting آورده شده است:
در برخی موارد، هر دو عملوند گسترش داده میشوند.
وقتی a -> (3,) و b -> (3, 1) داریم، a به a -> (3, 3) گسترش داده خواهد شد، و b به b -> (3, 3) گسترش داده خواهد شد.
در این حالت، جمع عنصر به عنصر منجر به یک ماتریس \(3 \times 3\) میشود
a = np.array([3, 6, 9])
b = np.array([2, 3, 4])
b.shape = (3, 1)
a + b
array([[ 5, 8, 11],
[ 6, 9, 12],
[ 7, 10, 13]])
در اینجا یک نمایش بصری از این عملیات broadcasting آورده شده است:
در حالی که broadcasting بسیار مفید است، گاهی اوقات میتواند گیجکننده به نظر برسد.
برای مثال، بیایید سعی کنیم a -> (3, 2) و b -> (3,) را جمع کنیم.
a = np.array(
[[1, 2],
[4, 5],
[7, 8]])
b = np.array([3, 6, 9])
a + b
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[62], line 7
3 [4, 5],
4 [7, 8]])
5 b = np.array([3, 6, 9])
6
----> 7 a + b
ValueError: operands could not be broadcast together with shapes (3,2) (3,)
ValueError به ما میگوید که عملوندها نمیتوانند با هم broadcast شوند.
در اینجا یک نمایش بصری برای نشان دادن اینکه چرا این broadcasting نمیتواند اجرا شود آورده شده است:
میبینیم که NumPy نمیتواند آرایهها را به یک اندازه گسترش دهد.
این به این دلیل است که، وقتی b از b -> (3,) به b -> (3, 3) گسترش مییابد، NumPy نمیتواند b را با a -> (3, 2) مطابقت دهد.
زمانی که به ابعاد بالاتر میرویم، موارد پیچیدهتر میشوند.
برای کمک به ما، میتوانیم از لیست قوانین زیر استفاده کنیم:
مرحله 1: وقتی ابعاد دو آرایه مطابقت ندارند، NumPy ابعادی را به آن یک که ابعاد کمتری دارد با اضافه کردن بعد (ابعاد) در سمت چپ ابعاد موجود گسترش میدهد.
برای مثال، اگر
a -> (3, 3)وb -> (3,)باشد، broadcasting یک بعد به سمت چپ اضافه خواهد کرد تاb -> (1, 3)شود؛اگر
a -> (2, 2, 2)وb -> (2, 2)باشد، broadcasting یک بعد به سمت چپ اضافه خواهد کرد تاb -> (1, 2, 2)شود؛اگر
a -> (3, 2, 2)وb -> (2,)باشد، broadcasting دو بعد به سمت چپ اضافه خواهد کرد تاb -> (1, 1, 2)شود (همچنین میتوانید این فرآیند را به عنوان عبور از مرحله 1 دو بار ببینید).
مرحله 2: وقتی دو آرایه بعد یکسانی دارند اما شکلهای متفاوت، NumPy سعی خواهد کرد ابعادی را که شاخص شکل آنها 1 است گسترش دهد.
برای مثال، اگر
a -> (1, 3)وb -> (3, 1)باشد، broadcasting ابعاد با شکل 1 را در هر دوaوbگسترش خواهد داد تاa -> (3, 3)وb -> (3, 3)شوند؛اگر
a -> (2, 2, 2)وb -> (1, 2, 2)باشد، broadcasting بعد اولbرا گسترش خواهد داد تاb -> (2, 2, 2)شود؛اگر
a -> (3, 2, 2)وb -> (1, 1, 2)باشد، broadcastingbرا در همه ابعاد با شکل 1 گسترش خواهد داد تاb -> (3, 2, 2)شود.
مرحله 3: پس از مرحله 1 و 2، اگر دو آرایه هنوز مطابقت نداشته باشند، یک
ValueErrorایجاد خواهد شد. برای مثال، فرض کنیدa -> (2, 2, 3)وb -> (2, 2)باشندطبق مرحله 1،
bبهb -> (1, 2, 2)گسترش خواهد یافت؛طبق مرحله 2،
bبهb -> (2, 2, 2)گسترش خواهد یافت؛میبینیم که پس از دو مرحله اول، آنها با یکدیگر مطابقت ندارند. بنابراین، یک
ValueErrorایجاد خواهد شد
10.6. قابلیت تغییر و کپی کردن آرایهها#
آرایههای NumPy انواع داده قابل تغییر هستند، مانند لیستهای Python.
به عبارت دیگر، محتویات آنها میتوانند پس از مقداردهی اولیه در حافظه تغییر یابند (جهش یابند).
این مناسب است اما، زمانی که با مدل نامگذاری و ارجاع Python ترکیب میشود، میتواند منجر به اشتباهات توسط مبتدیان NumPy شود.
در این بخش برخی از موضوعات کلیدی را بررسی میکنیم.
10.6.1. قابلیت تغییر#
قبلاً نمونههایی از قابلیت تغییر را در بالا دیدیم.
در اینجا یک مثال دیگر از جهش یک آرایه NumPy آورده شده است
a = np.array([42, 44])
a
array([42, 44])
a[-1] = 0 # عنصر آخر را به 0 تغییر دهید
a
array([42, 0])
قابلیت تغییر منجر به رفتار زیر میشود (که میتواند برای برنامهنویسان MATLAB شوکهکننده باشد…)
rng = np.random.default_rng()
a = rng.standard_normal(3)
a
array([0.06659573, 0.63505303, 1.18147425])
b = a
b[0] = 0.0
a
array([0. , 0.63505303, 1.18147425])
آنچه اتفاق افتاده این است که ما a را با تغییر دادن b تغییر دادهایم.
نام b به a متصل است و صرفاً یک ارجاع دیگر به
آرایه میشود (مدل انتساب Python با جزئیات بیشتر بعداً در دوره شرح داده شده است).
از این رو، حقوق برابری برای ایجاد تغییرات در آن آرایه دارد.
این در واقع معقولترین رفتار پیشفرض است!
این به این معنی است که ما فقط اشارهگرها را به داده منتقل میکنیم، به جای ایجاد کپی.
ایجاد کپی از نظر سرعت و حافظه گران است.
10.6.2. ایجاد کپی#
البته هنگام نیاز میتوان b را یک کپی مستقل از a ساخت.
این کار میتواند با استفاده از np.copy انجام شود
a = rng.standard_normal(3)
a
array([ 0.40324852, -0.72292693, 0.91765211])
b = np.copy(a)
b
array([ 0.40324852, -0.72292693, 0.91765211])
اکنون b یک کپی مستقل است (که deep copy نامیده میشود)
b[:] = 1
b
array([1., 1., 1.])
a
array([ 0.40324852, -0.72292693, 0.91765211])
توجه کنید که تغییر در b بر a تأثیر نگذاشته است.
10.7. ویژگیهای اضافی#
بیایید نگاهی به برخی ویژگیهای مفید دیگر NumPy بیندازیم.
10.7.1. توابع جهانی#
NumPy نسخههایی از توابع استاندارد log، exp، sin و غیره را فراهم میکند که به صورت عنصر به عنصر روی آرایهها عمل میکنند
z = np.array([1, 2, 3])
np.sin(z)
array([0.84147098, 0.90929743, 0.14112001])
این نیاز به حلقههای صریح عنصر به عنصر مانند موارد زیر را از بین میبرد
n = len(z)
y = np.empty(n)
for i in range(n):
y[i] = np.sin(z[i])
از آنجا که آنها به صورت عنصر به عنصر روی آرایهها عمل میکنند، این توابع گاهی اوقات توابع برداری شده نامیده میشوند.
در اصطلاح NumPy، آنها همچنین ufuncها، یا توابع جهانی نامیده میشوند.
همانطور که در بالا دیدیم، عملیات حسابی معمول (+، *، و غیره) نیز
به صورت عنصر به عنصر کار میکنند، و ترکیب اینها با ufunc ها مجموعه بسیار بزرگی از توابع سریع عنصر به عنصر را میدهد.
z
array([1, 2, 3])
(1 / np.sqrt(2 * np.pi)) * np.exp(- 0.5 * z**2)
array([0.24197072, 0.05399097, 0.00443185])
همه توابع تعریف شده توسط کاربر به صورت عنصر به عنصر عمل نخواهند کرد.
برای مثال، ارسال تابع f تعریف شده در زیر به یک آرایه NumPy باعث ValueError میشود
def f(x):
return 1 if x > 0 else 0
تابع NumPy np.where یک جایگزین برداری شده ارائه میدهد:
x = rng.standard_normal(4)
x
array([ 0.35612417, 1.37535143, 0.26586937, -0.05957849])
np.where(x > 0, 1, 0) # اگر x > 0 درست باشد 1 درج کنید، در غیر این صورت 0
array([1, 1, 1, 0])
همچنین میتوانید از np.vectorize برای برداری کردن یک تابع داده شده استفاده کنید
f = np.vectorize(f)
f(x) # ارسال همان بردار x مانند مثال قبلی
array([1, 1, 1, 0])
با این حال، این رویکرد همیشه همان سرعت یک تابع برداری شده با دقت بیشتر ساخته شده را به دست نمیآورد.
(بعداً خواهیم دید که JAX یک نسخه قدرتمند از np.vectorize دارد که میتواند و معمولاً کد بسیار کارآمدی تولید میکند.)
10.7.2. مقایسهها#
به عنوان یک قاعده، مقایسهها روی آرایهها به صورت عنصر به عنصر انجام میشوند
z = np.array([2, 3])
y = np.array([2, 3])
z == y
array([ True, True])
y[0] = 5
z == y
array([False, True])
z != y
array([ True, False])
وضعیت برای >، <، >= و <= مشابه است.
همچنین میتوانیم مقایسهها را با اسکالرها انجام دهیم
z = np.linspace(0, 10, 5)
z
array([ 0. , 2.5, 5. , 7.5, 10. ])
z > 3
array([False, False, True, True, True])
این به ویژه برای استخراج شرطی مفید است
b = z > 3
b
array([False, False, True, True, True])
z[b]
array([ 5. , 7.5, 10. ])
البته میتوانیم — و اغلب انجام میدهیم — این کار را در یک مرحله انجام دهیم
z[z > 3]
array([ 5. , 7.5, 10. ])
10.7.3. بستههای فرعی#
NumPy برخی از قابلیتهای اضافی مرتبط با برنامهنویسی علمی را از طریق بستههای فرعی خود ارائه میدهد.
قبلاً دیدهایم که چگونه میتوانیم با استفاده از
Generator تصادفی NumPy متغیرهای تصادفی تولید کنیم.
z = rng.standard_normal(10000) # تولید توزیع نرمال استاندارد
y = rng.binomial(10, 0.5, size=1000) # 1000 نمونه از Bin(10, 0.5)
y.mean()
np.float64(5.083)
یکی دیگر از بستههای فرعی رایج استفاده شده np.linalg است
A = np.array([[1, 2], [3, 4]])
np.linalg.det(A) # محاسبه دترمینان
np.float64(-2.0000000000000004)
np.linalg.inv(A) # محاسبه معکوس
array([[-2. , 1. ],
[ 1.5, -0.5]])
بسیاری از این قابلیتها همچنین در SciPy موجود است، مجموعهای از ماژولها که روی NumPy ساخته شدهاند.
نسخههای SciPy را به زودی با جزئیات بیشتر پوشش خواهیم داد.
برای فهرست جامعی از آنچه در NumPy موجود است به این مستندات مراجعه کنید.
10.7.4. چندنخی ضمنی#
قبلاً مفهوم موازیسازی از طریق چندنخی را مورد بحث قرار دادیم.
NumPy سعی میکند چندنخی را در بسیاری از کد کامپایل شده خود پیادهسازی کند.
بیایید نگاهی به یک مثال بیندازیم تا این را در عمل ببینیم.
قطعه کد بعدی مقادیر ویژه تعداد زیادی ماتریس تولید شده تصادفی را محاسبه میکند.
اجرای آن چند ثانیه طول میکشد.
n = 20
m = 1000
for i in range(n):
X = rng.standard_normal((m, m))
λ = np.linalg.eigvals(X)
اکنون، بیایید نگاهی به خروجی مانیتور سیستم htop در دستگاه ما در حالی که این کد در حال اجرا است بیندازیم:
میبینیم که 4 مورد از 8 CPU با سرعت کامل در حال اجرا هستند.
این به این دلیل است که روتین eigvals NumPy به زیبایی وظایف را تقسیم میکند و
آنها را به نخهای مختلف توزیع میکند.
10.8. تمرینها#
Exercise 10.1
عبارت چندجملهای زیر را در نظر بگیرید
قبلاً، یک تابع ساده p(x, coeff) برای ارزیابی (10.1) بدون در نظر گرفتن کارایی نوشتید.
اکنون یک تابع جدید بنویسید که همان کار را انجام میدهد، اما از آرایههای NumPy و عملیات آرایه برای محاسبات خود استفاده میکند، به جای هر شکلی از حلقه Python.
(چنین قابلیتی قبلاً به عنوان np.poly1d پیادهسازی شده است، اما به خاطر تمرین از این کلاس استفاده نکنید)
Hint
از np.cumprod() استفاده کنید
Solution to Exercise 10.1
این کد کار را انجام میدهد
def p(x, coef):
X = np.ones_like(coef)
X[1:] = x
y = np.cumprod(X) # y = [1, x, x**2,...]
return coef @ y
بیایید آن را تست کنیم
x = 2
coef = np.linspace(2, 4, 3)
print(coef)
print(p(x, coef))
# برای مقایسه
q = np.poly1d(np.flip(coef))
print(q(x))
[2. 3. 4.]
24.0
24.0
Exercise 10.2
فرض کنید q یک آرایه NumPy با طول n با q.sum() == 1 باشد.
فرض کنید که q یک تابع جرم احتمال را نمایش میدهد.
میخواهیم یک متغیر تصادفی گسسته \(x\) تولید کنیم به طوری که \(\mathbb P\{x = i\} = q_i\).
به عبارت دیگر، x مقادیری را در range(len(q)) میگیرد و x = i با احتمال q[i].
الگوریتم استاندارد (تبدیل معکوس) به صورت زیر است:
فاصله واحد \([0, 1]\) را به \(n\) زیر فاصله \(I_0, I_1, \ldots, I_{n-1}\) تقسیم کنید به طوری که طول \(I_i\) برابر با \(q_i\) باشد.
یک متغیر تصادفی یکنواخت \(U\) در \([0, 1]\) رسم کنید و \(i\) ای را برگردانید که \(U \in I_i\).
احتمال رسم \(i\) برابر با طول \(I_i\) است که برابر با \(q_i\) است.
میتوانیم الگوریتم را به صورت زیر پیادهسازی کنیم
from random import uniform
def sample(q):
a = 0.0
U = uniform(0, 1)
for i in range(len(q)):
if a < U <= a + q[i]:
return i
a = a + q[i]
اگر نمیتوانید ببینید این چگونه کار میکند، سعی کنید جریان را برای یک مثال ساده، مانند q = [0.25, 0.75] به خوبی بررسی کنید
رسم فواصل روی کاغذ کمک میکند.
تمرین شما این است که آن را با استفاده از NumPy سریعتر کنید و از حلقههای صریح اجتناب کنید
Hint
از np.searchsorted و np.cumsum استفاده کنید
اگر میتوانید، قابلیت را به عنوان یک کلاس به نام DiscreteRV پیادهسازی کنید، جایی که
دادههای یک نمونه از کلاس، بردار احتمالات
qاستکلاس یک متد
draw()دارد که یک نمونه را مطابق با الگوریتم توصیف شده در بالا برمیگرداند
اگر میتوانید، متد را طوری بنویسید که draw(k) k نمونه از q برگرداند.
Solution to Exercise 10.2
در اینجا اولین تلاش ما برای یک راهحل آورده شده است:
from numpy import cumsum
class DiscreteRV:
"""
یک آرایه از نمونهها را از یک متغیر تصادفی گسسته با بردار
احتمالات داده شده توسط q تولید میکند.
"""
def __init__(self, q, seed=None):
"""
آرگومان q یک آرایه NumPy است، یا شبیه آرایه، غیر منفی و جمع
به 1 میشود
"""
self.q = q
self.Q = cumsum(q)
self.rng = np.random.default_rng(seed)
def draw(self, k=1):
"""
k نمونه از q برمیگرداند. برای هر چنین نمونهای، مقدار i با
احتمال q[i] برگردانده میشود.
"""
return self.Q.searchsorted(self.rng.uniform(0, 1, size=k))
منطق واضح نیست، اما اگر وقت خود را بگذارید و آن را به آرامی بخوانید، درک خواهید کرد.
با این حال، در اینجا یک مشکل وجود دارد.
فرض کنید که q پس از ایجاد یک نمونه از discreteRV تغییر کند،
برای مثال با
q = (0.1, 0.9)
d = DiscreteRV(q)
d.q = (0.5, 0.5)
مشکل این است که Q بر این اساس تغییر نمیکند، و Q دادهای است که در متد draw استفاده میشود.
برای مقابله با این موضوع، یک گزینه محاسبه Q هر بار که متد draw فراخوانی میشود است.
اما این نسبت به محاسبه یکبار Q ناکارآمد است.
یک گزینه بهتر استفاده از descriptors است.
یک راهحل از کتابخانه quantecon با استفاده از descriptors که همانطور که میخواهیم رفتار میکند را میتوان در اینجا یافت.
Exercise 10.3
بحث قبلی ما در مورد تابع توزیع تجمعی تجربی را به خاطر بیاورید.
وظیفه شما این است که
متد
__call__را با استفاده از NumPy کارآمدتر کنید.یک متد اضافه کنید که ECDF را روی \([a, b]\) رسم میکند، جایی که \(a\) و \(b\) پارامترهای متد هستند.
Solution to Exercise 10.3
یک راهحل مثال در زیر داده شده است.
در واقع، ما فقط این کد را از QuantEcon گرفتهایم و یک متد plot به آن اضافه کردهایم
"""
ecdf.py را از QuantEcon اصلاح میکند تا یک متد plot اضافه کند
"""
class ECDF:
"""
تابع توزیع تجربی یک بعدی با توجه به یک بردار از
مشاهدات.
پارامترها
----------
observations : array_like
یک آرایه از مشاهدات
ویژگیها
----------
observations : array_like
یک آرایه از مشاهدات
"""
def __init__(self, observations):
self.observations = np.asarray(observations)
def __call__(self, x):
"""
ecdf را در x ارزیابی میکند
پارامترها
----------
x : scalar(float)
x که ecdf در آن ارزیابی میشود
برمیگرداند
-------
scalar(float)
کسری از نمونه کمتر از x
"""
return np.mean(self.observations <= x)
def plot(self, ax, a=None, b=None):
"""
ecdf را روی فاصله [a, b] رسم کنید.
پارامترها
----------
a : scalar(float), optional(default=None)
نقطه انتهای پایین فاصله رسم
b : scalar(float), optional(default=None)
نقطه انتهای بالای فاصله رسم
"""
# === انتخاب فاصله معقول اگر [a, b] مشخص نشده باشد === #
if a is None:
a = self.observations.min() - self.observations.std()
if b is None:
b = self.observations.max() + self.observations.std()
# === تولید رسم === #
x_vals = np.linspace(a, b, num=100)
f = np.vectorize(self.__call__)
ax.plot(x_vals, f(x_vals))
plt.show()
در اینجا یک مثال از استفاده آورده شده است
fig, ax = plt.subplots()
rng = np.random.default_rng()
X = rng.standard_normal(1000)
F = ECDF(X)
F.plot(ax)
Exercise 10.4
به یاد بیاورید که broadcasting در Numpy میتواند به ما کمک کند عملیات عنصر به عنصر را روی آرایهها با تعداد متفاوتی از ابعاد بدون استفاده از حلقههای for انجام دهیم.
در این تمرین، سعی کنید از حلقههای for برای تکرار نتیجه عملیات broadcasting زیر استفاده کنید.
قسمت 1: سعی کنید این مثال ساده را با استفاده از حلقههای for تکرار کنید و نتایج خود را با عملیات broadcasting زیر مقایسه کنید.
rng = np.random.default_rng(123)
x = rng.standard_normal((4, 4))
y = rng.standard_normal(4)
A = x / y
در اینجا خروجی آورده شده است
print(A)
قسمت 2: به سمت تکرار نتیجه عملیات broadcasting زیر حرکت کنید. در عین حال، سرعت broadcasting و حلقه for که پیادهسازی میکنید را مقایسه کنید.
برای این قسمت از تمرین میتوانید از توابع tic/toc از کتابخانه quantecon برای زمانسنجی اجرا استفاده کنید.
بیایید مطمئن شویم که این کتابخانه نصب شده است.
!pip install quantecon
اکنون میتوانیم بسته quantecon را import کنیم.
rng = np.random.default_rng(123)
x = rng.standard_normal((1000, 100, 100))
y = rng.standard_normal(100)
with qe.Timer("Broadcasting operation"):
B = x / y
Broadcasting operation: 0.0115 seconds elapsed
در اینجا خروجی آورده شده است
print(B)
Solution to Exercise 10.4
راهحل قسمت 1
rng = np.random.default_rng(123)
x = rng.standard_normal((4, 4))
y = rng.standard_normal(4)
C = np.empty_like(x)
n = len(x)
for i in range(n):
for j in range(n):
C[i, j] = x[i, j] / y[j]
نتایج را برای بررسی پاسخ خود مقایسه کنید
print(C)
همچنین میتوانید از array_equal() برای بررسی پاسخ خود استفاده کنید
print(np.array_equal(A, C))
True
راهحل قسمت 2
rng = np.random.default_rng(123)
x = rng.standard_normal((1000, 100, 100))
y = rng.standard_normal(100)
with qe.Timer("For loop operation"):
D = np.empty_like(x)
d1, d2, d3 = x.shape
for i in range(d1):
for j in range(d2):
for k in range(d3):
D[i, j, k] = x[i, j, k] / y[k]
For loop operation: 3.8551 seconds elapsed
توجه کنید که حلقه for مدت زمان بسیار بیشتری نسبت به عملیات broadcasting طول میکشد.
نتایج را برای بررسی پاسخ خود مقایسه کنید
print(D)
print(np.array_equal(B, D))
True