تحلیل فراوانی سیل

1. تحلیل فراوانی سیل موسوی ندوشنی مهر 1390

2. تحلیل فراوانی • این تحلیل شامل موارد زیر است: • داشتن یک سری آماری مشخص بر حسب تعریف متغیر مورد نظر • یافتن یک توزیع مناسب آماری برای برازش بر سری آماری مورد نظر • محاسبه دبی یا بارندگی، با دوره بازگشت

3. سری‌های آماری در دبی‌ها به دو بخش عمده تقسیم می‌شوند: دبی‌های لحظه‌ای سری دبی‌های لحظه‌ای حداکثر دبی‌های متوسط سری دبی‌های متوسط حداکثر روزانه سری دبی‌های متوسط ماهانه روزانه سری دبی‌های متوسط سالانه روزانه سری‌های آماری

4. انواع دبي‌هاي حداكثر در هيدرولوژي

5. روش ساختن سری‌های آماری • برای تحلیل فراوانی سیل از داده‌های لحظه‌ای استفاده می‌شود و در صورت عدم دسترسی به آنها از دبی‌های متوسط روزانه استفاده می‌گردد. • معمولاً از تمام داده‌های موجود برای ساخت سری‌های آماری استفاده نمی‌شود (به استثنای سری‌های زمانی) و باید به روش گزینشی عمل نمود. برای این کار دو روش موجود است. • بلوک حداکثرها B‌l‌o‌c‌kM‌a‌x‌i‌m‌a • روش ساختن بلوک‌ها به این صورت است که کل داده‌ها به اندازه‌های یکسان شامل توده‌ای از داده‌ها می‌باشند. مثلاً نمونه‌های حداکثر سالانه (سال آبی و یا سال تقویمی) • در این روش برای هر سال یک داده حداکثر انتخاب می‌شود. • استفاده از داده‌های بالاتر از یک آستانه معین Peak over threshold (POT) • در این روش برای هر سال یک یا بیش از یک داده اختیار می‌گردد. البته در سال‌های بسیار خشک ممکن است داده‌ای اخذ نشود.

6. استفاده از داده‌های بالاتر از یک آستانه معین

7. مقایسه روش AM و POT

8. مفاهیم اولیه در روش POT • فراترها Exceedences: • اگر سری زمانی x1,x2,… را در نظر بگیرید و مقدار آستانه u باشد، آنگاه x-u را برای xi>u مقادیر فراتر یافته نامند. • خوشه cluster: • به گروهی از مقادیر فراتر یافته متوالی خوشه گویند. • دوره cyclelength: • زمان بین دو فرود متوالی را دوره گویند. • طول اجرا runlength: • پریودی از دوره است که مقادیر بالاتر از آستانه وجود ندارد.

9. شکل شماتیک تعاریف اولیه روش POT

10. توزیع‌ مقادیر حدی • آشنایی با نظریه مقادیر حدی • در تحلیل فراوانی اگر از توزیع‌های کلاسیک استفاده شود، برای برآورد پارامترهایی نظیر  و  از تمام داده‌های مشاهده شده استفاده می‌شود و نه از داده‌های حدی. مقادیر مرکزی روی برآورد پارامترها تاثیر می‌گذارند. • اما اگر تحلیل داده‌های حدی مورد توجه است لزومی به استفاده از تمام داده‌ها نیست و فقط از مقادیر حدی استفاده می‌گردد. نظریه داده‌های حدی، توزیع‌های احتمال مورد نظر را به‌دست می‌دهد.

11. توزیع مقادیر حدی • فرض کنید که X1,X2,…,Xn متغیرهای تصادفی مستقل و دارای توزیع یکسان هستند و Xها یک تابع توزیع F می‌باشند. • فرض کنید که Mn=max{X1,X2,…,Xn} است. • در اینجا باید مقدار F برآورد گردد.

12. انواع توزیع مقادیر حدی • توزیع عمومی مقادیر حدی GeneralizeExtremeValues (GEV) به‌صورت زیر است. • که دارای سه پارامتر است.  پارامتر مکان،  پارامتر مقیاس که نشان‌دهنده پراکندگی است و  پارامتر شکل است و رفتار انتهای توزیع را نشان می‌دهد. • اگر =0 باشد توزیع گامبل و اگر >0 باشد توزیع فرشه و بالاخره اگر <0 باشد، توزیع ویبول است. • اگر <0، xG=+ باشد، آنگاه G از طرف چپ کران دارد. • اگر =0 تابع xGR است. • اگر >0، xG<+ باشد، آنگاه تابع G دارای کران طرف چپ نیست.

13. نمودار توزیع‌های GEV

14. شبیه‌سازی توسط R • داده‌های توزیع GEV را می‌توان با نرم‌افزار R شبیه‌سازی نمود. • با توجه به N(0,1) تعداد 365عدد تصادفی تولید می‌شود. • از اعداد تولید شده حداکثر آن بدست می‌آید. • برای بدست آوردن مقادیر حدی دیگر، گام‌های اول تا دوم صد بار تکرار می‌گردد. • در خاتمه 100نمونه حداکثر دارید که توزیع GEV بر آن قابل برازش است. اکنون به کدهای آن توجه کنید. • y = numeric(100) • for(i in 1:100) { • z = rnorm(365) • y[i] = max(z) • }

15. تابع چگالی پرتو Pareto • فرض کنید که X1,X2,…,Xnمتغیرهای تصادفی مستقل و دارای توزیع یکسان هستند و Xها یک تابع توزیع F می‌باشند. • فرض کنید که Mn=max{X1,X2,…,Xn}است. • از رابطه فوق لگاریتم گرفته می‌شود.

16. دنباله تابع ... • با استفاده از بسط تیلور log F(x)-[1-F(x)]است و با جایگذاری در رابطه بالا نتیجه می‌شود. • اکنون اگر X|X>u مقدار حدی را مشخص کند، برای uها که بقدر کافی بزرگ باشند. داریم:

17. دنباله تابع ... • اگر Y=X-u باشد، به‌طوری که X>u باشد، آنگاه داریم:

18. دنباله تابع ... • که در آن u پارامتر آستانه •  پارامتر مقیاس • * پارامتر مقیاس تغییر یافته که برابر است با •  پارامتر شکل است که رفتار انتهای توزیع را مشخص می‌کند. • اگر <0 باشد، مشاهده دارای کران هستند. اگر =0 باشد، توزیع نمایی خواهد بود. اگر >0 باشد، توزیع دارای انتهای کشیده است.

19. دنباله تابع ...

20. شبیه‌سازی توسط R • با توجه به N(0,1)تعداد 5000 عدد تصادفی تولید می‌شود. • داده‌ها به‌صورت نزولی مرتب می‌گردد • 100مقدار حداکثر انتخاب می‌شود. • در اینجا 100نمونه وجود دارد که از آستانه x[101] بزرگتر می‌باشند و با تابع Generalize Pareto Distribution(GPD)برازش می‌یابد. اکنون به کدهای مورد نظر توجه کنید. • x=rnorm(5000) • x=sort(x,decreasing=T) • y=x[1:100]

21. تعیین مقدار آستانه • برای یافتن مقدار u بهینه ابزارهایی وجود دارد که شرح آن در ادامه خواهد آمد. • ابتدا باید توجه داشت که تعداد افزایش‌ها بالاتر از یک آستانه متغیر تصادفی است که باید آنقدر بزرگ باشد (آستانه بزرگ) تا تابع احتمال این متغیر تصادفی دارای تابع احتمال پواسون باشد (اریب اندک). اگر فرض اخیر محقق نگردد، آنگاه توزیع حدیG‌P‌D برای مقادیر بالاتر از آستانه رخ نخواهد داد. • از طرفی تعداد داده‌های بالاتر از آستانه باید بقدر کافی بزرگ باشد (آستانه کوچک) که بتوان پارامترهای تابع G‌P‌D را برآورد نمود (واریانس اندک). • بنابراین مقدار آستانه باید به گونه‌ای صورت گیرد که توازنی بین اریب و واریانس برقرار گردد. بهر روی رویکرد تعیین مقدار آستانه یک رویکرد چند جانبه است.

22. آزمونplot Mean Residual life • اگر X GPD(u0,,)باشد، آنگاه • رابطه E[X-u|X>u] بر حسب u خطی است و متوسط مقادیر بالاتر از یک آستانه را نشان می‌دهد. • ابزار آزمون • برای عملی شدن آزمون اخیر باید رابطه خطی بالا در حالت نمونه محاسبه نمود.

23. دنباله آزمون ... • پس از محاسبه، در یک دستگاه مختصات زوج‌های {u,e(u): u<xmax} رسم کنید. قسمت خطی گراف فوق را مشخص کنید. به عنوان مثال به کُدهای زیر توجه کنید. • library(POT) • data(ardieres) • ardieres <- clust(ardieres, 4, 10 / 365, clust.max = TRUE) • flows <- ardieres[, "obs"] • mrlplot(flows, u.range=c(4,15),lty=c(3,1,3),col = c("red","black","red")) • segments(5.8,4.7,10,4.7,lwd=2,col="blue",lty='dashed') • text(10.5,3.5,"linear") • arrows(9.7,3.4,8.4,4.5,angle=20,length=0.1)

24. دنباله آزمون ...

25. معیار پارامترها برای انتخاب مقدار آستانه • اگر X GPD(u0,,)باشد، آنگاه برای کلیه مقادیر u>u0 و X|X>u برای توزیع G‌P‌D داریم: • مقدار پارامتر شکل (u) مستقل از مقادیر آستانه‌ها u است. • مقدار پارامتر مقیاس *(u) مستقل از مقادیر آستانه‌ها u است. • ابزار آزمون: • نقاط زیر را نمایش دهید. • قسمت ثابت را مشخص کنید.

26. کدهای مربوط به تغییرات پارامترها بر حسب آستانه • library(POT) • data(ardieres) • ardieres <- clust(ardieres, 4, 10 / 365, clust.max = TRUE) • flows <- ardieres[, "obs"] • par(mfrow=c(1,2)) • tcplot(flows, u.range = c(0, 15), which=1 ) • segments(5,-0.05,15,-0.05, lwd=2,col='red') • tcplot(flows, u.range = c(0, 15), which=2 ) • segments(5,0.25,15,0.25, lwd=2,col='red')

27. دنباله معیار پارامترها ...

28. دنباله معیار پارامترها ...

29. اندیس پراکندگی • همان‌طور که قبلاً ذکر شد، متغیر شمارش تعداد نقاط حداکثر بالای آستانه دارای توزیع پواسون است. در این توزیع مقدار میانگین و واریانس یکسان است. بنابراین نسبت آنها برابر واحد است. البته مقدار واحد، مقداری نظری است و مقدار نمونه‌ای این نسبت حول عدد واحد در نوسان است. • فاصله اطمینان این نوسان با توزیع مربع خی با درجه آزادی M-1 تعریف می‌گردد. که M تعداد کل سال‌های داده‌های مشاهده شده است. فاصله اطمینان به‌صورت زیر تعریف می‌گردد.

30. ابزار آزمون اندیس پراکندگی • ابتدا باید نقاط زیر را رسم نمود. • library(POT) • data(ardieres) • ardieres <- clust(ardieres, 1.5, 10 / 365, clust.max = TRUE) • diplot(ardieres,xlim=c(1,22),ylim=c(0.5,2.0)) • abline(h=1,lty=2,col='red')

31. نمودار اندیس پراکندگی

32. جمع‌بندی تعیین مقدار آستانه • library(POT) • data(ardieres) • events0 <- clust(ardieres, u = 1.5, tim.cond = 8/365, clust.max = TRUE) • par(mfrow = c(2, 2)) • mrlplot(events0[, "obs"]) • abline(v = 6, col = "green") • diplot(events0) • abline(v = 6, col = "green") • tcplot(events0[, "obs"], which = 1) • abline(v = 6, col = "green") • tcplot(events0[, "obs"], which = 2) • abline(v = 6, col = "green")

33. جمع بندی نمودارها

34. عدم قطعیت در مدلسازی • در مدلسازی عدم قطعیت مربوط به موارد زیر می گردد. • داده های ورودی به مدل (نمونه‌گیری، خطای اندازه‌گیری و تعداد کم داده‌ها بخصوص در حالت حدی) • عدم قطعیت در برآورد پارامترهای مدل • عدم قطعیت در مدلسازی: ساختار مدل به گونه‌ای است که ممکن است نتواند پدیده مورد نظر را مدل نماید.

35. برآورد پارامترهای مدل • اگر X متغیر تصادفی با مقادیر حقیقی باشد. تابع تجمعی آن به‌صورت زیر است. • در حالتی که F(.) یک تابع پیوسته باشد و تابع معکوس داشته باشد x(.) را تابع چندک متغیر X نامند. • برای 0<u<1 مقدار x(u) یکتا و در رابطه زیر صدق می‌کند. • در ادبیات محیط زیست و مهندسی چندک‌ها بر حسب دوره بازگشت بیان می‌شود. • در عمل وقتی توزیع یک پدیده مشخص شد، آنگاه با یک مجموعه‌ای از پارامترها روبرو هستیم. • بنابراین تابع چندک برای پارمترها به صورت زیر است.

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

37. دنباله نکویی برآوردگرها • برای بیان نکویی دو شاخص سویه و واریانس به شرح زیر فرموله می‌گردد. • سویه: • وقتی برآوردگر سویه ندارد که یعنی باشد. • اما ممکن است چند برآوردگر برای یک پارامتر همه بدون سویه باشند. لذا برای مقایسه آنها از معیار واریانس استفاده می‌شود. هر چقدر واریانس کمتر باشد، برآوردگر کاراتر (efficiency) است.

38. روش‌های برآورد پارامترهای آماری • در آمار روشهای مختلفی برای برآورد پارامترها وجود دارد. که به شرح زیر است. • روش گشتاورها Mehtod of Moments • روش حداکثر درستنمایی Method of Maximum Likelihood • روش گشتاورهای وزنی احتمال probability weighted moments • روش گشتاورهای خطی L-momentsMethod of

39. روش گشتاورها • در این روش گشتاور مقدار مرکزی • برای گشتاورهای مرتبه بالاتر داریم: • انحراف معیار که پراکندگی را اندازه می‌گیرد. • گشتاورهای بدون بعد مرتبه‌های بالاتر • چولگی skewness • کشیدگی kurtosis

40. دنباله روش گشتاورها • اگر تابع معکوس u=F(x) در نظر بگیرید. • یک تابع از یک متغیر تصادفی، خود یک متغیر تصادفی است و امید ریاضی آن به شکل زیر محاسبه می‌شود.

41. روش حداکثر درستنمایی • از منظر آماری برداری از داده‌ها y=(y1,…,yn) یک نمونه تصادفی است که از جامعه‌ی آماری نامعلوم استخراج شده است. هدف از تحلیل داده‌ها این است که جامعه‌ی مذکور معلوم گردد. در آمار هر جامعه‌ای با یک توزیع احتمال مشخص می‌شود. • هر توزیع دارای پارامتر و یا پارامترهایی است که مقدار منحصر به‌فرد دارند. • فرض کنید که f(y|w) تابع چگالی احتمالی است که احتمال هر داده از بردار y را به ازاء پارامتر w را به‌دست می‌دهد. طبیعی است که توزیع‌ها می‌توانند دارای بیش از یک پارامتر باشند، بنابراین می‌توان w=(w1,…,wk)منظور نمود. اگر داده‌ها از هم مستقل باشند می‌توان رابطه f(y|w) را به‌صورت زیر نوشت:

42. دنباله روش حداکثر درستنمایی • اکنون ایده فوق برای ساده‌ترین حالت که n=k=1 است، بیان می‌گردد. فرض کنید که y تعداد موفقیت‌ها در 10 آزمون برنولی را نشان می‌دهد و احتمال پیروزی در هر آزمون برابر 0.2 است که مقدار پارامتر w را نشان می‌دهد. • بنابراین: • و به‌طور کلی می‌توان نوشت:

43. تابع درستنمایی • برای مجموعه‌ای از مقادیر پارامترها بعضی از داده‌ها محتمل‌تر از داده‌های دیگر هستند. برای مثال در مساله قبلی y=2 محتمل‌تر از y=5 برای w=0.2 (0.302 در مقابل 0.026) است. در عمل ابتداً با داده‌ها روبرو هستیم، لذا با عکس مساله مواجه خواهیم شد. یعنی به ازاء داده‌های مشاهده شده و مدل مورد نظر در بین توابع چگالی احتمال به دنبال بیشترین شباهت داده‌ها تولید شده با داده‌های • مشاهده شده هستیم. • برای حل مساله معکوس از تابع درستنمایی استفاده می‌شود. • یعنی • که در آن L(w|y) نشان‌دهنده درستنمایی پارامتر w بر حسب داده‌های مشاهده شده است.

44. تابع درستنمایی • اگر مثال قبلی را در نظر بگیرید. • اختلاف مهمی بین تابع چگالی احتمال f(y|w) و تابع درستنمایی L(w|y) وجود دارد.

45. تابع درستنمایی

46. کدهای شکل قبل • rm(list=ls());p=seq(0,1,.01) • w <- 4 • n <- 10 • y=dbinom(x=w,size=n,prob=p) • par(mfrow=c(1,2)) • plot(p,y,typ="l", lwd=2, xlab="parameter w", ylab="L(w|n=10,y=4)", main="", ylim=c(0,0.275),col='blue') • max_x <- which.max(y) • text(p[max_x],max(y)+0.01,"MLE") • plot(0:10,dbinom(0:10, size=n, prob=0.5),type='h', xlab="x", ylab="p(y|w=0.5)", col='red')

47. دنباله تابع درستنمایی • برآورد به روش M‌L‌E می‌تواند موجود و یکتا نباشد. در اینجا فرض می‌شود که شرط وجود ویکتایی برای محاسبه M‌L‌E برقرار است. • برای محاسبه، لگاریتم تابع درستنمایی حداکثر می‌گردد. لگاریتم و خود تابع درستنمایی همنوا هستند و دارای نتایج مرتبطی می‌باشند. دلیل لگاریتم گرفتن این است که • حاصل‌ضرب به دلیل توابع حاشیه‌ای مستقل داده‌ها به حاصل‌جمع تبدیل می‌گردد و مقدار احتمال غالباً کوچک هست و کارکردن با حاصل‌ضرب، همراه با اعداد کوچک در حل مسایل عددی اشکالاتی را باعث می‌شود.

48. دنباله تابع درستنمایی • فرض کنید که تابعl‌n L(w|y) مشتق‌پذیر است، اگر wM‌L‌Eموجود باشد، باید معادله دیفرانسیل زیر برقرار باشد. • که در آن wi=wi,M‌L‌Eبه ازاء تمام مقادیر i=1,…,k است. • در مشتق مرتبه اول ممکن است هم مقدار حداکثر و هم مقدار حداقل حاصل شود. برای اطمینان از حداکثر بودن مشتق مرتبه دوم l‌n L(w|y) محاسبه می‌شود که باید برای مقادیر wi=wi,M‌L‌Eبه ازاء تمام مقادیر منفی i=1,…,k باشد.

49. دنباله تابع درستنمایی • مشتق مرتبه دوم عبارتست از: • اکنون برای مثال به تابع درستنمایی L(w|n=10,y=7) توجه کنید. اگر لگاریتم منظور گردد. • اکنون مشتق آن محاسبه می‌گردد. • که پس از حل معادله اخیر مقدار wM‌L‌E=0.7 می‌شود.

50. دنباله تابع درستنمایی • برای اطمینان از این که مقدار مورد نظر حداکثر است و نه حداقل، مشتق مرتبه دوم محاسبه می‌شود و مقدار آن به ازاء w=wM‌L‌E محاسبه می‌گردد. • در عمل، برای برآورد M‌L‌E راه حل تحلیلی وجود ندارد. بالاخص وقتی که P‌D‌F دارای پارامترهای متعددی باشد. لذا با توجه به روش‌های عددی و بهینه‌سازی باید برآورد انجام گردد.