آزمون نرمال بودن جارک برا (Jarque-Bera Test) — به زبان ساده

در مباحث آماری، اطلاع از توزیع داده‌ها، امری است که نوع و روش تحلیل آن‌ها را مشخص می‌کند. اغلب با توجه به نرمال بودن توزیع احتمال داده‌ها، روش‌های آمار پارامتری را برای انجام تحلیل‌ها به کار می‌بریم. در مقابل اگر از توزیع داده‌ها کمترین میزان اطلاع را داشته باشیم، روش‌های آمار ناپارامتری به کارمان خواهند آمد. در این نوشتار به کمک «آزمون نرمال بودن جارک برا» (Jarque-Bera Test) این شرط را مورد بررسی قرار می‌دهیم. در این بین برای انجام محاسبات این آزمون، از زبان‌های محاسبات آماری و برنامه‌نویسی R و پایتون کمک می‌گیریم.

برای آشنایی بیشتر با دیگر شیوه‌های بررسی نرمال بودن توزیع داده‌ها، نوشتار آزمون نرمال بودن داده (Normality Test) — پیاده سازی در پایتون و نمودار چندک چندک (Q-Q plot) — به زبان ساده را مطالعه کنید. همچنین خواندن آمار پارامتری و ناپارامتری — انتخاب روش های تحلیل خالی از لطف نیست.

آزمون نرمال بودن جارک برا (Jarque-Bera Test)

انتخاب روش تحلیل داده‌ها، یکی از مباحثی است که بسیاری از محققین در زمینه آمار را به خود مشغول کرده است. اغلب روش‌های استنباط آماری، مبتنی بر توزیع نرمال برای جامعه آماری است. چنین شیوه‌هایی اغلب به روش‌های استنباط «آمار پارامتری» (Parametric Statistics) شهرت دارند که مناسب برای داده‌هایی با توزیع نرمال هستند. ولی جامعه نرمال به عنوان یک جامعه ایده‌ال شناخته شده و کمتر داده‌های واقعی از این توزیع برخوردارند. در نتیجه ابداع روش‌هایی که بدون در نظر گرفتن توزیع (یا حداقل توزیع نرمال) برای جامعه آماری به کار گرفته می‌شوند، ضروری به نظر می‌رسد. چنین شاخه‌ای از آمار که مرتبط به روش‌های استنباط ناپارامتری هستند، «آمار ناپارامتری» (Non-Parametric Statistics) نامیده می‌شود.

فیلم آموزش شبیه سازی متغیرهای تصادفی و وابسته در متلب در فرادرس

کلیک کنید

تعیین اینکه آیا نمونه جمع‌آوری شده، مربوط به یک جامعه آماری نرمال است، توسط آزمون‌های مختلفی صورت می‌گیرد. اغلب این روش‌ها را با نام آزمون‌های نرمالیتی (Normality Test) می‌شناسند. در این قسمت به آزمون نرمال بودن جارک برا می‌پردازیم که بخصوص برای بررسی نرمال‌بودن باقی‌مانده‌ها در سری زمانی و رگرسیون به کار می‌رود.

یکی از «آزمون‌های نیکویی برازش» (Goodness of Fit test)، »آزمون جارک برا» (Jarque-Bera Test) است که براساس شاخص‌های تقارن چولگی (Skewness) و کشیدگی (Kurtosis)، مطابقت با توزیع نرمال را می‌سنجد.

این آزمون توسط «کارلوس جارک» (Carlos Jarque) اقتصاددان مکزیکی و «آنیل برا» (Anil Bera) اقتصاددان آمریکایی در سال 1987 در مقاله‌ای با عنوان «آزمون نرمال بودن مشاهدات و باقی‌مانده‌های رگرسیون» (Test for Normality of Observations and Regression Residuals) معرفی شد. به کمک این آزمون و آماره آن می‌توانیم نرمال بودن داده‌ها را مورد بررسی قرار دهیم.

آماره آزمون نرمال بودن جارک برا و خصوصیات آن

همانطور که گفته شد، این آزمون براساس اندازه‌های مربوط به تقارن عمودی و افقی توزیع نرمال عمل می‌کند که توسط چولگی و کشیدگی اندازه‌گیری می‌شوند. فرض کنید $$S$$، «چولگی نمونه‌ای» (Sample Skewness) و $$K$$ نیز «کشیدگی نمونه‌ای» (Sample Kurtosis) برای مشاهدات $$x_1,x_2,\ldots,x_n$$ باشند. یعنی:

$$\large S={\frac {{\hat {\mu }}_{3}}{{\hat {\sigma }}^{3}}}={\frac {{\frac 1n}\sum _{{i=1}}^{n}(x_{i}-{\bar {x}})^{3}}{\left({\frac 1n}\sum _{{i=1}}^{n}(x_{i}-{\bar {x}})^{2}\right)^{{3/2}}}}$$

و

$$\large {\displaystyle K={\frac {{\hat {\mu }}_{4}}{{\hat {\sigma }}^{4}}}={\frac {{\frac {1}{n}}\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{4}}{\left({\frac {1}{n}}\sum _{i=1}^{n}(x_{i}-{\bar {x}})^{2}\right)^{2}}}}$$

آنگاه، آماره آزمون جارک برحسب این دو پارامتر، به صورت زیر تعیین می‌شود.

$$\large {\displaystyle {\mathit {JB}}={\frac {n}{6}}\left(S^{2}+{\frac {1}{4}}(K-3)^{2}\right)}$$

نکته: منظور از $$\mu_3$$ و $$\mu_4$$ گشتاورهای مرکزی سوم و چهارم و $$\hat{\mu}_3$$ و $$\hat{\mu}_4$$ نیز  برآورد این گشتاورها (Sample Central Moments) است. همچنین $$\bar{x}$$ و $$\hat{\sigma}$$ نیز میانگین و  انحراف استاندارد نمونه‌ای هستند.

اگر داده‌ها، از یک جامعه آماری با توزیع نرمال باشند، مقدار آماره $$JB$$ کوچک و تقریبا برابر با صفر است و توزیع احتمالاتی برای این آماره به صورت مجانبی (Asymptotically)، توزیع کای ۲ (Chi Square Distribution) با ۲ درجه آزادی (Degree of Freedom) خواهد بود. به این ترتیب از این آماره برای بررسی نرمال بودن یک مجموعه داده می‌توان استفاده کرد.

بر همین اساس فرض صفر (Null Hypothesis) به این صورت در نظر گرفته می‌شود که همزمان مقدار چولگی و کشیدگی اصلاح شده نمونه‌ای صفر باشند. البته منظورمان از «کشیدگی اصلاح شده» (Excess Kurtosis) همان $$K-3$$ است. در نتیجه اگر کشیدگی برابر با ۳ باشد، کشیدگی اصلاح شده نیز صفر خواهد بود.

در نتیجه هر گونه انحراف از مقادیر این شاخص‌ها نسبت به توزیع نرمال، براساس توزیع کای ۲ اندازه‌گیری شده و مقدار احتمال (p-Value)، برای آزمون نرمال بودن به کار گرفته می‌شود. با کوچک بودن مقدار آماره $$JB$$ رای به نرمال بودن داده‌ها داده و فرض صفر با توجه به سطح آزمون ($$\alpha$$) و توزیع کای ۲ رد نمی‌شود ولی اگر مقدار این آماره از مقدار بحرانی آزمون بزرگتر شود، فرض نرمال بودن جامعه آماری، توسط این نمونه رد خواهد شد.

برای اندازه نمونه‌های کوچک، مقادیر توزیع کای ۲ خیلی حساس هستند و با اندک تغییری در مقادیر نمونه، مقدار این آماره تغییر خواهد کرد، در نتیجه برای نمونه‌های کوچک، اغلب این آزمون رای به رد فرض صفر خواهد داد و نرمال بودن داده‌ها تایید نخواهد شد. ولی با افزایش حجم نمونه، قابلیت این آزمون برای تشخیص نرمال بودن داده‌ها بیشتر می‌شود.

آزمون نرمال بودن جارک برا در R

به منظور اجرای آزمون نرمال بودن جارک برا در زبان برنامه‌نویسی R، باید از کتابخانه سری زمانی (tseries) استفاده کنید. کدهای زیر به همراه یک مثال و خروجی دستورات مربوطه در ادامه به آزمون نرمال بودن یک مجموعه داده پرداخته‌اند.

ابتدا ۱۰۰ مشاهده یا نمونه از توزیع نرمال تولید شده و به کمک آزمون جارک برا، نرمال بودن جامعه مربوط به نمونه مورد بررسی قرار می‌گیرد. در ادامه، همین عمل برای توزیع یکنواخت پیوسته (Uniform Distribution) نیز انجام می‌شود.

1# ---------------------------------------
2# Install package and call tseries library
3# ---------------------------------------
4install.packages("tseries")
5library(tseries)
6
7# ---------------------------------------
8# Create Normal data
9# --------------------------------------- 
10x <- rnorm(100)  # null
11jarque.bera.test(x)
12
13# ---------------------------------------
14# Create Uniform data
15# ---------------------------------------
16
17x <- runif(100)  # alternative
18jarque.bera.test(x)

نکته: خط اول برای نصب کتابخانه tseries لازم است ولی در اجراهای بعدی از آن صرف‌نظر کنید.

پس از اجرا این کد در زبان برنامه‌نویسی R، خروجی به صورت زیر خواهد بود.

1> jarque.bera.test(x)
2
3	Jarque Bera Test
4
5data:  x
6X-squared = 0.55372, df = 2, p-value = 0.7582
7
8> x <- runif(100)  # alternative
9> jarque.bera.test(x)
10
11	Jarque Bera Test
12
13data:  x
14X-squared = 6.432, df = 2, p-value = 0.04012

همانطور که مشاهده می‌کنید، آماره آزمون برابر با 0٫55372 بوده و همچنین مقدار احتمال (p-value) در قسمت اول که مربوط به آزمون داده‌های شبیه‌سازی شده از توزیع نرمال است برابر با ۰٫7582 بدست آمده است. در نتیجه فرض نرمال بودن توسط این نمونه برای جامعه آماری رد نخواهد شد.

ولی در قسمت دوم با یک نمونه ۱۰۰ تایی از توزیع یکنواخت، آزمون را انجام داده‌ایم که مشخص است مقدار آماره آزمون 6٫432 بوده و نسبت به توزیع کای ۲ با دو درجه آزادی مقدار بزرگی است. از طرفی مقدار احتمال (p-value) نیز کوچک شده (0٫04012) است که با در نظر گرفتن سطح آزمون $$\alpha=0.05$$، فرض صفر رد می‌شود. در نتیجه این جامعه آماری مربوط به این داده‌ها دارای توزیع نرمال نخواهد بود.

نکته: همانطور که مشاهده می‌کنید، آماره و آزمون جارک برا، فقط مطابقت توزیع داده‌ها را با توزیع نرمال می‌سنجد و نمی‌تواند مشخص کند که داده‌ها دارای چه توزیعی هستند.

در مسائل مربوط به رگرسیون و سری زمانی، تعداد پارامترهای مدل نیز در محاسبه آماره آزمون جارک برا، دخیل هستند. در این صورت آماره آزمون به صورت زیر در خواهد آمد.

$$\large {\mathit {JB}}={\frac {n-k}{6}}\left(S^{2}+{\frac 14}(K-3)^{2}\right)$$

که در آن $$n$$ تعداد مشاهدات و $$k$$ نیز تعداد پارامترهای برآورد شده است.

فرض کنید قرار است یک مدل سری زمانی ایجاد کنید. همانطور که می‌دانید، یکی از شرط‌های صحت مدل در چنین حالتی، نرمال بودن باقی‌مانده‌های آن است. در اینجا نیز با کمک یک مجموعه داده (hcip)، مدل سری زمانی (ARIMA) را برازش کرده، سپس نرمال بودن باقی‌مانده را مورد بررسی قرار می‌دهیم. این مجموعه داده در کتابخانه tseries قرار دارد و نیازی به دریافت آن از طریق اینترنت نیست. کدهای زیر به این منظور تهیه شده است. در این بین از کتابخانه tsoutliers کمک گرفته‌ایم.

1# ---------------------------------------
2# Install package and call tseries library
3# ---------------------------------------
4install.packages("tsoutliers")
5library(tsoutliers)
6# fit an ARIMA model to the HICP 011600 series
7# ARIMA(0,1,0)(2,0,1) was chosen by forecast::auto.arima(ic = "bic")
8# normality of the residuals is rejected at the 5% significance level
9# due to an excess of kurtosis
10data("hicp")
11y <- log(hicp[["011600"]])
12fit1 <- arima(y, order = c(0, 1, 0), seasonal = list(order = c(2, 0, 1)))
13
14jarque.bera.test(residuals(fit1))
15
16# fit ARIMA model for the same series including outliers that were
17# detected by "tsoutliers" and for the model chosen by "auto.arima"
18# normality of the residuals is not rejected at the 5% significance level
19# after including the presence of outliers
20mo <- outliers(c("AO", "AO", "LS", "LS"), c(79, 210, 85, 225))
21xreg <- outliers.effects(mo, length(y))
22fit2 <- arima(y, order = c(1, 1, 0), seasonal = list(order = c(2, 0, 2)), 
23              xreg = xreg)
24jarque.bera.test(residuals(fit2))
25

در قسمت اول داده‌ها بدون شناسایی نقاط پرت (Outliers) در مدل سری زمانی ARIMA قرار گرفته‌اند. همانطور که مشخص است، باقی‌مانده‌ها با توجه به کوچک بودن مقدار p-value، نرمال نخواهند بود.

در قسمت دوم، به کمک تابع outliers، نقاط پرت از مدل خارج شده و بهبود زیادی ایجاد شده بطوری که فرض نرمال بودن باقی‌مانده رد نخواهد شد. خروجی را در ادامه مشاهده می‌کنید.

1> fit1 <- arima(y, order = c(0, 1, 0), seasonal = list(order = c(2, 0, 1)))
2> jarque.bera.test(residuals(fit1))
3
4	Jarque Bera Test
5
6data:  residuals(fit1)
7X-squared = 11.592, df = 2, p-value = 0.00304
8
9> jarque.bera.test(residuals(fit2))
10
11	Jarque Bera Test
12
13data:  residuals(fit2)
14X-squared = 0.15761, df = 2, p-value = 0.9242

آزمون نرمال بودن جارک برا در پایتون

به منظور پیاده‌سازی آزمون نرمال بودن جارک برا در زبان برنامه‌نویسی پایتون از کتابخانه statsmodels و stattools استفاده می‌شود.

فیلم آموزش درس رگرسیون 2 – جامع و کاربردی در فرادرس

کلیک کنید

تابع آزمون جارک برا در کتابخانه statsmodels به صورت زیر است.

$$\large statsmodels.stats.stattools.jarque\_bera(resids, axis=0)$$

در نتیجه کافی است که مقادیر باقی‌مانده یا داده‌ها را در پارامتر resids مشخص کنید تا آزمون جارک برا صورت گیرد.

نظیر این گونه محاسبه نیز توسط کتابخانه scipy.stats صورت می‌گیرد. کافی است تابع scipy.stats.jarque_bera را اجرا کنید.

برای مثال، کد زیر یک سری داده تصادفی از توزیع نرمال (Normal) و توزیع رایلی (Rayleight) تولید کرده و به کمک تابع stats.jarque_bera آن‌ها را مورد آزمون قرار می‌دهد.

1>>> from scipy import stats
2>>> np.random.seed(987654321)
3>>> x = np.random.normal(0, 1, 100000)
4>>> y = np.random.rayleigh(1, 100000)
5>>> stats.jarque_bera(x)
6(4.7165707989581342, 0.09458225503041906)
7>>> stats.jarque_bera(y)
8(6713.7098548143422, 0.0)

با توجه به خروجی‌ها مشخص است که داده‌هایی که در متغیر $$x$$ ذخیره شده‌اند (شبیه‌سازی شده از توزیع نرمال)، توسط آزمون جارک برا نیز نرمال تشخیص داده شده‌‌اند، زیرا مقدار آماره آزمون کوچک (4٫7165) و مقدار احتمال (0٫09458) بزرگتر از خطای نوع اول ($$\alpha=0.05$$) است. از طرفی فرض نرمال بودن برای داده‌هایی که از توزیع رایل ایجاد شده‌‌اند در سطح آزمون 0٫05، رد می‌شود.