۳ گام ساده در R برای فیت کردن مناسب‌ترین تابع توزیع احتمال

بسیار خوب. به مبحث بسیار شیرین آمار می‌رسیم. کسانی که تحلیل داده را بعنوان شغل اصلی خود انتخاب می‌کنند، از اهمیت بی‌ چون و چرای دانش آماری مطلع هستند. بدون سواد آماری، یک تحلیل‌گر داده ممکن نیست بتواند دید صحیحی از آنچه که انجام می‌دهد، داشته باشد.

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

• اگر درباره توزیع‌های احتمال هیچ اطلاعاتی از پیش ندارید، پیشنهاد می‌کنم این ویدیو را تماشا کنید.
• من بیشتر اوقات از پایتون استفاده می‌کنم ولی توابع از پیش آماده‌ای که در R وجود دارند، کار را برای تحلیل‌های آماری بسیار آسوده‌تر می‌کنند. (همچنین بد نیست به این بهانه سری به R بزنیم تا بببینم آنجا چه خبر است!)

گام اول: آماده سازی مجموعه داده

اگر قبلا R را نصب نکرده‌اید می‌توانید از اینجا به راحتی نصب و با ویرایشگر پیش‌فرض از آن استفاده کنید و یا از برنامه‌ی R Studio بهره ببرید.

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

data = rnorm(100, 20, 5)

مجموعه‌ داده‌ی بدست آمده، به شکل اعداد زیر است:

همین داده را بصورت یک هیستوگرام می‌توانیم نمایش دهیم:

hist(data)

حالا همین داده را به شکل چگالی احتمال تغییر می‌دهیم. به بردار طول‌ها دقت کنید که تغییر پیدا می‌کند و احتمال وقوعِ آن مقدار می‌شود. در واقع در این حالت مجموع اعداد روی نمودار یک است. توجه کنید که ما توزیع‌‌های احتمال مختلف را روی این نمودار فیت می‌کنیم.

plot(density(data))

گام دوم: فیت کردن توزیع‌های احتمال مختلف

برای فیت‌ کردن توزیع‌های احتمال مختلف می‌توانیم از پکیج fitdistrplus در R استفاده کنیم. ابتدا این پکیج را نصب و بارگذاری می‌کنیم.

install.packages(&quotfitdistrplus&quot)

library(&quotfitdistrplus&quot)

دستور زیر یک دید کلی و جامع درباره‌ی میزان نزدیکی داده‌ی ما با انواع توابع توزیع‌های احتمال ارائه می‌دهد.

descdist(data, boot=1000)

علائم مختلف روی نمودار، هرکدام نشانگر نوع خاصی از تابع توزیع هستند. برای مثال مثلث به معنای توزیع uniform، ستاره توزیع Normal، مثبت توزیع Logistic، خط چین درشت توزیع Gamma و خط چین ریز توزیع Lognormal است. نقطه‌ی آبی به هر نشانه‌ای که نزدیک‌تر باشد، احتمالا داده‌ی ما آن توزیع را دنبال می‌کند. همانطور که می‌بینید، نمودار Cullen and Frey پیش‌بینی می‌کند که توزیع داده ما نزدیک به تابع توزیع نرمال است. در واقع این نمودار یک دید کلی (intuition) درباره شبیه‌ترین توزیع به داده‌ی ما ارائه می‌دهد.

حالا ۴ توزیع مختلف نرمال (normal)، گاما (Gamma)، لاگ نرمال (Lognormal) و نمایی (exponential) را روی داده، فیت می‌کنیم تا نتایج را ببینیم. کتابخانه fitdistrplus بصورت پیشفرض از شیوه‌ی برآورد درست‌نمایی بیشینه (Maximum Liklihood Estimation) برای محاسبه مناسب‌ترین پارامترها استفاده می‌کند. با این حال شما می‌توانید از دیگر شیوه‌ها هم استفاده کنید.

fnorm = fitdist(data, &quotnorm&quot)
fgamma = fitdist(data, &quotgamma&quot)
flognorm = fitdist(data, &quotlnorm&quot)
fexp = fitdist(data, &quotexp&quot)

اول نگاهی به توزیع‌های فیت شده بیاندازیم. برای این‌کار یک لیست و یک وکتور ساده، اول از اسامی توزیع‌ها و دومی خودِ توزیع‌ها ایجاد می‌کنیم:

plot.legend = c(&quotNormal&quot, &quotGamma&quot,  &quotLogNormal&quot, &quotExponential&quot)
dist.list = list(fnorm, fgamma, flognorm, fexp)

با اجرای دستور کوتاه زیر نمودار توزیع‌های تئوریک بر روی داده واقعی خود را می‌بینیم:

denscomp(dist.list,legendtext = plot.legend)

همچنین می‌توانیم نمودار توزیع CDF یا همان تابع توزیع تجمعی (Cumulative distribution function) را هم ببینیم.

cdfcomp(dist.list, legendtext = plot.legend)

از هر دو این نمودارها مشخص است که توزیع نمایی (Expnential) توزیع مناسبی برای مجموعه داده ما نیست. ولی بین سه توزیع ما بقی کدام‌یک توزیع بهتری است یا (باصطلاح بهتر فیت شده است؟)

گام سوم: مقایسه بین توزیع‌ها و انتخاب بهترین

شیوه‌های مختلف آماری برای محاسبه توزیعِ "بهتر" وجود دارد. برای مثال آزمون کولموگروف–اسمیرنف (KS-TEST) یکی از شناخته‌شده‌ترین‌ها آزمایش‌ها است. هر قدر که عدد حاصل از آزمون KS پایین‌تر باشد، نشان‌ دهنده‌ی این است که توزیع بهتر فیت شده است.

gofstat(dist.list , fitnames=plot.legend)

نتیجه‌ی آزمون KS به ما نشان می‌دهد که توزیع نرمال بهترین فیت برای داده است. (البته انتظاری جز این نداشتیم چراکه داده را بصورت رندم از توزیع نرمال ایجاد کرده بودیم.) بقیه آزمون‌های انجام شده نیز ثابت می‌کنند که توزیع نرمال، توزیع مناسب‌تری برای این مجموعه داده است.

اگر با این اعداد قانع نمی‌شوید و برای اطمینان می‌خواهید p-value را برای هر آزمون ببینید، می‌توانید از دستور مستقیم آزمون KS بصورت مجزی برای هر توزیع فیت شده استفاده کنید. من برای توزیع نرمال و نمایی جهت مقایسه این کار را انجام می‌دهم.

ks.test(data, pexp, fexp$estimate[1], fexp$estimate[2])

در نتیجه‌ی این آزمون P-value صفر است. یعنی تابع نمایی اصلا توزیع مناسبی برای این مجموعه داده نیست و فرض پیروی از این توزیع رد می‌شود. D همان مقداری است که در جدول بالا هم نمایش داده شده است.

حال P-value را برای توزیع نرمال بدست می‌آوریم:

ks.test(data, pnorm, fnorm$estimate[1], fnorm$estimate[2])

در نتیجه‌ی این آزمون P-value معادل 0.736 است. این عدد بسیار بالاتر از مقدار مد نظر است و فرض دنبال کردن دیتا از این تابع توزیع احتمال را تایید می‌کند.

خلاصه

مباحث آماری گسترده و بی‌پایان هستند. اگر علاقه داریم در حوزه‌ی تحلیل داده کار کنیم، مجبوریم با آن‌ها سروکار داشته باشیم. در این نوشته‌ی کوتاه من تلاش کردم یک موضوع مبتدی را گام به گام توضیح بدهم. خوشحال می‌شوم اگر اشتباهی کرده باشم در کامنت‌ها یادآوری کنید و اگر سوالی داشتید بپرسید. قول می‌دهم برای پیدا کردن پاسخ تلاشم را کنم.