Transcript Document
تحلیل فراوانی سیل
موسوی ندوشنی
مهر 1390
1
تحلیل فراوانی
این تحلیل شامل موارد زیر است:
داشتن یک سری آماری مشخص بر حسب تعریف متغیر
مورد نظر
یافتن یک توزیع مناسب آماری برای برازش بر سری آماری
مورد نظر
محاسبه دبی یا بارندگی ،با دوره بازگشت
2
سریهای آماری
سریهای آماری در دبیها به دو بخش عمده تقسیم
میشوند:
دبیهای لحظهای
سری دبیهای لحظهای حداکثر
دبیهای متوسط
سری دبیهای متوسط حداکثر روزانه
سری دبیهای متوسط ماهانه روزانه
t d / 2
سری دبیهای1
ساالنه روزانه
متوسط
av e (t )
f (u ) du
t d / 2
3
d
f (u ) du
t d / 2
t d / 2
v ol (t )
انواع دبيهاي حداكثر در هیدرولوژي
4
روش ساختن سریهای آماری
برای تحلیل فراوانی سیل از دادههای لحظهای استفاده میشود
و در صورت عدم دسترسی به آنها از دبیهای متوسط
روزانه استفاده میگردد.
معموالً از تمام دادههای موجود برای ساخت سریهای
آماری استفاده نمیشود (به استثنای سریهای زمانی) و باید
به روش گزینشی عمل نمود .برای این کار دو روش موجود
است.
بلوک حداکثرها Block Maxima
روش ساختن بلوکها به این صورت است که کل دادهها به
اندازههای یکسان شامل تودهای از دادهها میباشند .مثالً نمونههای
حداکثر ساالنه (سال آبی و یا سال تقویمی)
استفاده از دادههای باالتر از یک آستانه معین Peak over threshold
)(POT
5
در این روش برای هر سال یک داده حداکثر انتخاب میشود.
در این روش برای هر سال یک یا بیش از یک داده اختیار میگردد .البته
در سالهای بسیار خشک ممکن است دادهای اخذ نشود.
استفاده از دادههای باالتر از یک آستانه معین
6
مقایسه روش AMو POT
7
مفاهیم اولیه در روش POT
فراترها :Exceedences
اگر سری زمانی … x1,x2,را در نظر بگیرید و مقدار آستانه u
باشد ،آنگاه x-uرا برای xi>uمقادیر فراتر یافته نامند.
خوشه :cluster
به گروهی از مقادیر فراتر یافته متوالی خوشه گویند.
دوره :cycle length
زمان بین دو فرود متوالی را دوره گویند.
طول اجرا :run length
پریودی از دوره است که مقادیر باالتر از آستانه وجود ندارد.
8
شکل شماتیک تعاریف اولیه روش POT
9
توزیع مقادیر حدی
آشنایی با نظریه مقادیر حدی
در تحلیل فراوانی اگر از توزیعهای کالسیک استفاده شود،
برای برآورد پارامترهایی نظیر و از تمام دادههای
مشاهده شده استفاده میشود و نه از دادههای حدی .مقادیر
مرکزی روی برآورد پارامترها تاثیر میگذارند.
اما اگر تحلیل دادههای حدی مورد توجه است لزومی به
استفاده از تمام دادهها نیست و فقط از مقادیر حدی استفاده
میگردد .نظریه دادههای حدی ،توزیعهای احتمال مورد نظر
را بهدست میدهد.
10
توزیع مقادیر حدی
فرض کنید که X1,X2,…,Xnمتغیرهای تصادفی
مستقل و دارای توزیع یکسان هستند و Xها یک تابع
توزیع Fمیباشند.
فرض کنید که } Mn=max{X1,X2,…,Xnاست.
]x
n
,X
]x
x,
n
1
x ] P [X
P [X
x ]
n
در اینجا باید مقدار Fبرآورد گردد.
n
11
1
n
P [M
P [X
]) [ F ( x
M n bn
P
) x G (x
an
انواع توزیع مقادیر حدی
توزیع عمومی مقادیر حدی Values (GEV) Generalize Extreme
بهصورت زیر است.
1/
x
G ( x ) exp 1
که دارای سه پارامتر است .پارامتر مکان ،پارامتر مقیاس که
نشاندهنده پراکندگی است و پارامتر شکل است و رفتار انتهای
توزیع را نشان میدهد.
اگر =0باشد توزیع گامبل و اگر >0باشد توزیع فرشه و
باالخره اگر <0باشد ،توزیع ویبول است.
اگر xG=+ ،<0باشد ،آنگاه Gاز طرف چپ کران دارد.
اگر =0تابع xGRاست.
اگر xG<+ ،>0باشد ،آنگاه تابع Gدارای کران طرف چپ
نیست.
12
نمودار توزیعهای GEV
13
شبیهسازی توسط R
دادههای توزیع GEVرا میتوان با نرمافزار Rشبیهسازی
نمود.
14
با توجه به ) N(0,1تعداد 365عدد تصادفی تولید میشود.
از اعداد تولید شده حداکثر آن بدست میآید.
برای بدست آوردن مقادیر حدی دیگر ،گامهای اول تا دوم صد بار
تکرار میگردد.
در خاتمه 100نمونه حداکثر دارید که توزیع GEVبر آن قابل برازش
است .اکنون به کدهای آن توجه کنید.
)y = numeric(100
{ )for(i in 1:100
)z = rnorm(365
)y[i] = max(z
}
تابع چگالی پرتو Pareto
فرض کنید که X1,X2,…,Xnمتغیرهای تصادفی مستقل و
دارای توزیع یکسان هستند و Xها یک تابع توزیع F
میباشند.
فرض کنید که } Mn=max{X1,X2,…,Xnاست.
1/
x
n
F ( x ) exp 1
از رابطه فوق لگاریتم گرفته میشود.
1/
15
x
n log F ( x ) 1
دنباله تابع ...
با استفاده از بسط تیلور ]) log F(x)-[1-F(xاست و با
جایگذاری در رابطه باال نتیجه میشود.
1/
1
x
1 F ( x ) 1
n
اکنون اگر X|X>uمقدار حدی را مشخص کند ،برای uها
که بقدر کافی بزرگ باشند .داریم:
16
... دنباله تابع
: آنگاه داریم، باشدX>u بهطوری که، باشدY=X-u اگر
P r( X u y | X u )
y
1
(u )
1
y
1
1
(u y )
1
n
1
(u )
1
n
1/
1/
1/
1/
17
دنباله تابع ...
که در آن uپارامتر آستانه
پارامتر مقیاس
) (u
1
* پارامتر مقیاس تغییر یافته که برابر است با
پارامتر شکل است که رفتار انتهای توزیع را مشخص
میکند.
اگر <0باشد ،مشاهده دارای کران هستند .اگر =0
باشد ،توزیع نمایی خواهد بود .اگر >0باشد ،توزیع
دارای انتهای کشیده است.
18
*
دنباله تابع ...
19
شبیهسازی توسط 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
20
تعیین مقدار آستانه
21
برای یافتن مقدار uبهینه ابزارهایی وجود دارد که شرح آن در
ادامه خواهد آمد.
ابتدا باید توجه داشت که تعداد افزایشها باالتر از یک آستانه متغیر
تصادفی است که باید آنقدر بزرگ باشد (آستانه بزرگ) تا تابع
احتمال این متغیر تصادفی دارای تابع احتمال پواسون باشد (اریب
اندک) .اگر فرض اخیر محقق نگردد ،آنگاه توزیع حدی GPDبرای
مقادیر باالتر از آستانه رخ نخواهد داد.
از طرفی تعداد دادههای باالتر از آستانه باید بقدر کافی بزرگ باشد
(آستانه کوچک) که بتوان پارامترهای تابع GPDرا برآورد نمود
(واریانس اندک).
بنابراین مقدار آستانه باید به گونهای صورت گیرد که توازنی بین
اریب و واریانس برقرار گردد .بهر روی رویکرد تعیین مقدار
آستانه یک رویکرد چند جانبه است.
آزمون Mean Residual life plot
اگر ) X GPD(u0,,باشد ،آنگاه
for u u 0
u
1
E [X u | X u ]
رابطه ] E[X-u|X>uبر حسب uخطی است و متوسط مقادیر
باالتر از یک آستانه را نشان میدهد.
ابزار آزمون
برای عملی شدن آزمون اخیر باید رابطه خطی باال در حالت
نمونه محاسبه نمود.
nu
) u
i
(x
i 1
22
1
nu
eˆ (u )
... دنباله آزمون
u<xmax}
{u,e(u): در یک دستگاه مختصات زوجهای، پس از محاسبه
به عنوان مثال به. قسمت خطی گراف فوق را مشخص کنید.رسم کنید
.ُکدهای زیر توجه کنید
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='dashe
d')
text(10.5,3.5,"linear")
arrows(9.7,3.4,8.4,4.5,angle=20,length=0.1)
23
دنباله آزمون ...
24
معیار پارامترها برای انتخاب مقدار آستانه
اگر ) X GPD(u0,,باشد ،آنگاه برای کلیه مقادیر
u>u0و X|X>uبرای توزیع GPDداریم:
مقدار پارامتر شکل (u)مستقل از مقادیر آستانهها u
است.
مقدار پارامتر مقیاس ) *(uمستقل از مقادیر آستانهها u
است.
ابزار آزمون:
نمایش ) (u
زیر را :u x
نقاط
دهید* and u , .
m ax
25
قسمت ثابت را مشخص کنید.
m ax
x
u , (u ) :u
کدهای مربوط به تغییرات پارامترها بر حسب آستانه
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')
26
دنباله معیار پارامترها ...
27
دنباله معیار پارامترها ...
28
اندیس پراکندگی
همانطور که قبالً ذکر شد ،متغیر شمارش تعداد نقاط
حداکثر باالی آستانه دارای توزیع پواسون است .در این
توزیع مقدار میانگین و واریانس یکسان است .بنابراین
نسبت آنها برابر واحد است .البته مقدار واحد ،مقداری
نظری است و مقدار نمونهای این نسبت حول عدد واحد
در نوسان است.
فاصله اطمینان این نوسان با توزیع مربع خی با درجه
آزادی M-1تعریف میگردد .که Mتعداد کل سالهای
اطمینان به
دادههای مشاهده شده است .
صورت
فاصله
,
M 1
M 1
زیر تعریف میگردد.
2
1 / 2 , M 1
29
2
/ 2 , M 1
ابزار آزمون اندیس پراکندگی
.ابتدا باید نقاط زیر را رسم نمود
2
ˆ
u
u
,
:
u
x
m ax
ˆ
u
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')
30
نمودار اندیس پراکندگی
31
جمعبندی تعیین مقدار آستانه
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")
32
جمع بندی نمودارها
33
عدم قطعیت در مدلسازی
در مدلسازی عدم قطعیت مربوط به موارد زیر می
گردد.
داده های ورودی به مدل (نمونهگیری ،خطای اندازهگیری و
تعداد کم دادهها بخصوص در حالت حدی)
عدم قطعیت در برآورد پارامترهای مدل
عدم قطعیت در مدلسازی :ساختار مدل به گونهای است که
ممکن است نتواند پدیده مورد نظر را مدل نماید.
34
برآورد پارامترهای مدل
اگر Xمتغیر تصادفی با مقادیر حقیقی باشد .تابع تجمعی آن
بهصورت زیر است.
0 F (x ) 1
] F ( x ) P r[ X x
در حالتی که ) F(.یک تابع پیوسته باشد و تابع معکوس داشته باشد
) x(.را تابع چندک متغیر Xنامند.
برای 0<u<1مقدار ) x(uیکتا و در رابطه زیر صدق میکند.
F ( x (u )) u
در ادبیات محیط زیست و مهندسی چندکها بر حسب دوره بازگشت
بیان میشود.
در عمل وقتی توزیع یک پدیده مشخص شد ،آنگاه با یک
مجموعهای از پارامترها روبرو هستیم.
1 , p
بنابراین تابع چندک برای پارمترها به صورت زیر است.
) p
35
x (u ; 1 ,
نکویی برآوردگرها
پارامترهای مجهول از روی دادهها برآورد میشوند.
ها برآوردگر تابعی از دادهها
برای یک سری از داده ˆ
است.
ˆ
برآوردگر یک متغیر تصادفی است و دارای توزیع
احتمال میباشد.
ˆ
نکویی برآوردگر بستگی به این دارد که چقدر به
نزدیک شده است .انحراف این از یکدیگر را میتوان به
دو بخش تقسیم نمود.
36
اریب یا سویه :میزان تمایل برآوردگر به کوچکتر و یا
بزرگتر شدن از مقدار واقعی پارامتر
میزان پراکندگی :انحراف تصادفی تخمین از مقدار واقعی
حتی اگر برآوردگر سویه نداشته باشد.
دنباله نکویی برآوردگرها
برای بیان نکویی دو شاخص سویه و واریانس به شرح
زیر فرموله میگردد.
) bias(ˆ ) E (ˆ
سویه:
وقتی برآوردگر سویه ندارد کهbias( ˆ )=0
یعنیE (ˆ )
باشد.
اما ممکن است چند برآوردگر برای یک پارامتر همه بدون
سویه باشند .لذا برای مقایسه آنها از معیار واریانس استفاده
میشود .هر چقدر واریانس کمتر باشد ،برآوردگر کاراتر
) (efficiencyاست.
37
روشهای برآورد پارامترهای آماری
در آمار روشهای مختلفی برای برآورد پارامترها وجود
دارد .که به شرح زیر است.
روش گشتاورها Mehtod of Moments
روش حداکثر درستنمایی Method of Maximum Likelihood
روش گشتاورهای وزنی احتمال probability weighted moments
روش گشتاورهای خطی Method of L-moments
38
روش گشتاورها
مرکزی E ( x
)
در این روش گشتاور مقدار
برای گشتاورهای مرتبه باالتر داریم:
r 2, 3,
r
) r E (X
انحراف معیار که پراکندگی را اندازه میگیرد.
گشتاورهای بدون بعد مرتبههای باالتر
1/ 2
2
) E (X
r / 2
r /2
چولگی skewness
کشیدگی kurtosis
39
3 / 2
3/ 2
4 / 2
2
1/ 2
2
دنباله روش گشتاورها
اگر تابع معکوس ) u=F(xدر نظر بگیرید.
1
x (u )du
E (X )
0
یک تابع از یک متغیر تصادفی ،خود یک متغیر تصادفی
است و امید ریاضی آن به شکل زیر محاسبه میشود.
1
g ( x (u )) du
0
40
g ( x ) f ( x ) dx
g ( x ) dF ( x )
E [ g ( X )]
روش حداکثر درستنمایی
از منظر آماری برداری از دادهها ) y=(y1,…,ynیک
نمونه تصادفی است که از جامعهی آماری نامعلوم استخراج
شده است .هدف از تحلیل دادهها این است که جامعهی مذکور
معلوم گردد .در آمار هر جامعهای با یک توزیع احتمال
مشخص میشود.
هر توزیع دارای پارامتر و یا پارامترهایی است که مقدار
منحصر بهفرد دارند.
فرض کنید که ) f(y|wتابع چگالی احتمالی است که احتمال
هر داده از بردار yرا به ازاء پارامتر wرا بهدست میدهد.
طبیعی است که توزیعها میتوانند دارای بیش از یک پارامتر
باشند ،بنابراین میتوان )w=(w1,…,wkمنظور نمود .اگر
)( yf(y|wراf
رابطه ) f ( y
توان , , y
می | w
باشند) f
مستقل | ( y
هم w ) f
| wاز( y
داده)ها
41بهصورت زیر نوشت:
n
n
1
1
n
1
دنباله روش حداکثر درستنمایی
اکنون ایده فوق برای سادهترین حالت که n=k=1است،
بیان میگردد .فرض کنید که yتعداد موفقیتها در 10
آزمون برنولی را نشان میدهد و احتمال پیروزی در هر
آزمون برابر 0.2است که مقدار پارامتر wرا نشان
میدهد.
! 10
f ( y | n ,w )
)(0.2) (1 0.2
بنابراینy :
0,1, ,10
10 y
y
! ) y !(10 y
و بهطور کلی میتوان نوشت:
,n
42
0 w 1; y 0,1,
ny
) (1 w
y
w
!n
! ) y !( n y
f ( y | n ,w )
تابع درستنمایی
برای مجموعهای از مقادیر پارامترها بعضی از دادهها محتملتر از
دادههای دیگر هستند .برای مثال در مساله قبلی y=2محتملتر از
y=5برای 0.302( w=0.2در مقابل )0.026است .در عمل
ابتداً با دادهها روبرو هستیم ،لذا با عکس مساله مواجه خواهیم شد.
یعنی به ازاء دادههای مشاهده شده و مدل مورد نظر در بین توابع
چگالی احتمال به دنبال بیشترین شباهت دادهها تولید شده با دادههای
مشاهده شده هستیم.
برای حل مساله معکوس از تابع درستنمایی استفاده میشود.
یعنی
) L (w | y ) f ( y | w
که در آن ) L(w|yنشاندهنده درستنمایی پارامتر wبر حسب
دادههای مشاهده شده است.
43
تابع درستنمایی
اگر مثال قبلی را در نظر بگیرید.
w (1 w ) 0 w 1
3
7
! 10
!7 !3
L (w | n 10, y 7 ) f ( y 7 | n 10, w )
اختالف مهمی بین تابع چگالی احتمال ) f(y|wو تابع
درستنمایی) L(w|yوجود دارد.
44
تابع درستنمایی
45
کدهای شکل قبل
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')
46
دنباله تابع درستنمایی
برآورد به روش MLEمیتواند موجود و یکتا نباشد.
در اینجا فرض میشود که شرط وجود ویکتایی برای
محاسبه MLEبرقرار است.
برای محاسبه ،لگاریتم تابع درستنمایی حداکثر میگردد.
لگاریتم و خود تابع درستنمایی همنوا هستند و دارای
نتایج مرتبطی میباشند .دلیل لگاریتم گرفتن این است که
حاصلضرب به دلیل توابع حاشیهای مستقل دادهها به
حاصلجمع تبدیل میگردد و مقدار احتمال غالبا ً کوچک هست
و کارکردن با حاصلضرب ،همراه با اعداد کوچک در حل
مسایل عددی اشکاالتی را باعث میشود.
47
دنباله تابع درستنمایی
فرض کنید که تابع) ln L(w|yمشتقپذیر است ،اگر
wMLEموجود باشد ،باید معادله دیفرانسیل زیر برقرار
) ln L (w | y
باشد.
i
w
که در آن wi=wi,MLEبه ازاء تمام مقادیر i=1,…,k
است.
در مشتق مرتبه اول ممکن است هم مقدار حداکثر و هم
مقدار حداقل حاصل شود .برای اطمینان از حداکثر بودن
مشتق مرتبه دوم ) ln L(w|yمحاسبه میشود که باید
برای مقادیر wi=wi,MLEبه ازاء تمام مقادیر منفی
i=1,…,k 48باشد.
دنباله تابع درستنمایی
مشتق مرتبه دوم عبارتست از:
) ln L (w | y
2
0
i
w
اکنون برای مثال به تابع درستنمایی )L(w|n=10,y=7
منظور10گردد.
توجه کنید .اگر لگاریتم
!
ln L (w | n 10, y 7 ) ln
) 7 ln w 3 ln(1 w
!7 !3
اکنون مشتق آن محاسبه میگردد.
0
7 10w
) w (1 w
3
1 w
7
w
)d ln L (w | n 10, y 7
dw
که پس از حل معادله اخیر مقدار wMLE=0.7میشود.
49
دنباله تابع درستنمایی
برای اطمینان از این که مقدار مورد نظر حداکثر است و
نه حداقل ،مشتق مرتبه دوم محاسبه میشود و مقدار آن
به ازاء w=wMLEمحاسبه میگردد.
47.62 0
50
3
2
) (1 w
7
2
w
) d ln L (w | n 10, y 7
2
2
dw
در عمل ،برای برآورد MLEراه حل تحلیلی وجود
ندارد .باالخص وقتی که PDFدارای پارامترهای
متعددی باشد .لذا با توجه به روشهای عددی و
بهینهسازی باید برآورد
انجام گردد.
دنباله تابع درستنمایی
فرض کنید که Xیک متغیر تصادفی نرمال باشد که
میانگینش مجهول است و واریانس آن یک است .اگر
X1,X2,…,Xnیک نمونه تصادفی Xباشد تابع
درستنمایی عبارتست از:
(x i ) /2
2
e
2
2
51
2
لگاریتمش عبارتست از:
) (x i
n
1
/2
2
ln(2 )
n
(x )
i 1
i 1
) (x i
e
n
2
X
f
L ( )
n /2
1
2
K ( ) ln L ( )
دنباله تابع درستنمایی
مشتق مرتبه اول و دوم عبارتست از:
n
i
x
)
i
(x
dK
d
2
n 0
d K
2
d
فرض کنید که Xیک متغیر تصادفی یکنواخت در فاصله
) (0,که در آن پارامتر مثبت مجهولی است ،باشد.
بنابراین
1
0x
52
f X (x )
دنباله تابع درستنمایی
تابع درستنمایی آن برای یک نمونه به اندازه nعبارت است
1
از:
L ( )
منحنی نمایش آن به گونهای است که شیب در هیچ جا صفر
نمیشود .لذا احتیاجی به مشتق گرفتن و مساوی صفر قرار
دادن نیست .اما توجه داسته باشید که هر قدر به صفر
نزدیکتر شود تابع درستنمایی از نظر مقدار افزایش مییابد.
پس ) L(به ازاء کوچکترین مقداری که میتواند قبول
کند حداکثر میشود .واضح است که نمیتواند کوچکتر از
بزرگترین مقداری باشد که در نمونه ما رخ دهد (زیرا فرض
این است کهˆ x
هر عضو نمونه همان توزیعی را داردکه X
دارد).
است.
بنابراین
n
) (n
53
دنباله تابع درستنمایی
مثال :نمونههای زیر از یک توزیع یکنواخت در فاصله
) (0,اخذ شده است.
0.5,0.6,0.1,2,0.9,1.6,0.7,0.8,1,1.5
با استفاده از روش گشتاور تخمین پارامتر دو برابر
میانگین است 20.94=1.88 .اما همانطور که مشاهده
میشود .نمونه با مقدار 2از آن بزرگتر است و در
فاصله ) (0,قرار نمیگیرد.
با استفاده از روش حداکثر درستنمایی تخمین پارامتر
برابر 2است که مشکل روش گشتاورها را ندارد.
54
گشتاورهای خطی
این روش جانشینی برای بیان شکل توزیع احتمال است .به
لحاظ تاریخی از تغییر گشتاورهای وزنی احتمال حاصل شده
است.
گشتاورهای وزنی احتمال یک متغیر تصادفی با تابع توزیع
تجمعی ،توسط Greenwoodو همکاران ( )1979به
صورت کمیتهای
M
E X F ( X ) 1 F ( X )
M
خاص ،موارد ویژهی مفید،
تعریف شدند .Mبه طور
هستند.
و
گشتاورهای وزنی احتمال
برای یک توزیع که تابع چندک دارد ،روابط زیر رانتیجه
r
s
p
55
1,0 , s
s
1, r ,0
r
p , r ,s
دنباله گشتاورهای خطی
اگر احتمال تجاوز ) (Exceedanceدر نظر باشد.
1
x (u )(1 u ) s d u
s
s
1 F ( X )
p
E X
1,0 , s
s M
0
اگر احتمال عدم تجاوز ) (Non-Exceedanceدر
نظر باشد.
1
r
F
(
X
)
x
(
u
)
u
du
r
r
p
E X
1, r , 0
r M
0
گشتاورهای وزنی احتمال و به عنوان اساس
روشهای برآورد پارامترهای توزیعهای احتمال مورد
استفاده قرار گرفتهاند .با این حال تفسیر مستقیم آنها
بهعنوان شاخصهای مقیاس و شکل توزیع احتمال دشوار
56است.
دنباله گشتاورهای خطی
برای مثال ،برآوردهایی از پارامترهای مقیاس توزیعها،
یا 2
مضاربی از 2
هستند.
این ترکیبات خطی از انتگرالهای وزن یافته با
مجموعهای از چند جملهایهای متعامد پدید میآیند.
را تعریف میکنیم که در آنها
ما چند جملهایهای
و از شرایط زیر پیروی میکنند:
1
0
1
0
یک چند جملهای rدرجه بر uحسب است.
p r (1) 1
*
1
r s
p r (u ) p s (u ) du 0
*
*
0
57
دنباله گشتاورهای خطی
گشتاورهای وزنی sو rتوسط روابط زیر با هم مرتبط
میشوند.
s
( 1)
i
s
i
i
i
(
)1
i
58
s
i 0
r
r
i 0 i
r
در سال 1990گشتاور خطی که ترکیب خطی از
گشتاورهای وزنی احتمال است ،توسط هاسکینگ معرفی
شد.
گشتاورهای خطی مناسبتر از گشتاورهای وزنی احتمال
هستند ،زیرا بهصورت مستقیم میتواند اندازههایی مثل
مقیاس و شکل توزیع احتمال را تبیین کنند.
دنباله گشتاورهای خطی
گشتاورهای خطی بر حسب گشتاورهای وزنی احتمال
و بهصورت زیر معرفی شد.
r
p r ,k r
*
k 0
که در آن
r
p r ,k r
*
r
r 1 1
k 0
) r r k ( 1
! r k
2
k
k
! ) k ! ( r k
r k
r k
)p r , k ( 1
*
در گشتاورهای خطی نیز مثل گشتاورهای عادی
زیر است .
نسبتهایی معرفی میشود ،که بهصور
r 3
59
r
2
r
2
1
دنباله گشتاورهای خطی
1 اندازه مکان ،اندازه مقیاس یا پراکندگی است.
LCv گشتاورهای باالتر خطی را نشان میدهد که نباید با
نماد ضریب تغییرات اشتباه شود.
3 اندازه چولگی LCsو 4اندازه کشیدگی LCkبیان
میکند.
نسبت گشتاورهای خطی نمونه را با tو trنشان داده
میشود که با قرار دادن برآورد rبدست میآید.
اگر r3باشد آنگاه قدر مطلق rکمتر از واحد است.
60
دیاگرام نسبت گشتاورهای خطی
61
شرح دیاگرام نسبت گشتاورهای خطی
نمودار نسبتهای گشتاور خطی .توزیعهای دو و سه پارامتری به ترتیب
با نقاط و خطوط نشان داده شدهاند .نشانههای اختصاری به کار رفته
برای توزیعهای مختلف عبارتند از:
Eبرای توزیع نمایی،
Gبرای توزیع گامبل،
Lبرای توزیع لجستیک،
Nبرای توزیع نرمال،
Uبرای توزیع یکنواخت،
GPAبرای پرتو تعمیم یافته،
GEVبرای مقادیر حدی تعمیم یافته،
GLOبرای لجستیک تعمیم یافته،
LN3برای لوگ نرمال سه پارامتری،
PE3برای پیرسون نوع .3
ناحیهی هاشور خورده شامل مقادیر ممکن است که از رابطهی زیر به
دست میآیند.
62
63
برآورد پارامترهای توزیع PGDبرحسب توزیعهای خطی
پارامترهای این توزیع بر حسب گشتاورهای خطی به
صورت زیر است 1 Γ 1 / .
1
Γ (1 ) /
3
)
6(1 2
10 1 3
1 2
2 1 2
2 1 3
1 2
5 1 4
3
4
در اینجا (.)نشان دهندهی تابع گاما است.
64
t x 1e t dt
0
Γ x
بستههای نرمافزاری Rراجع به گشتاورهای خطی و تحلیل منطقهای
بسته lmom
این بسته شامل توابع مربوط به گشتاورهای خطی است .محاسبهی
گشتاورهای خطی توزیعها ،برآورد پارامترها ،ترسیم نمودار
نسبتهای گشتاور خطی و ،...از جمله امکاناتی است که این بسته
برای کاربران فراهم میکند.
بسته lmomco
65
نویسندهJ. R. M. Hosking :
نویسندهWilliam H. Asquith :
این بسته جهت اجرای گشتاورهای خطی ،شامل برآورد
گشتاورهای خطی ،برآورد گشتاورهای وزنی احتمال ،برآورد
پارامترهای توزیعهای آشنا و نه چندان آشنای متعدد ،و برآورد
گشتاورهای خطی این توزیعها از روی پارامترهای آنها نوشته
شده است.
دنباله بستههای نرمافزاری Rراجع به گشتاور...
بسته RFA
نویسندهMathieu Ribatet :
این بسته شامل چندین تابع جهت اجرای یک تحلیل فراوانی
منطقهای است .چند نمونه داده نیز در این بسته ارائه شده است.
بسته nsRFA
66
نویسندهAlberto Viglione :
این بسته شامل مجموعهای از ابزارهای آماری برای کاربردهای
عملی روشهای تحلیل فراوانی منطقهای در هیدرولوژی است .این
بسته بیشتر با تمرکز روی روش مقدار (سیالب) نمایه،
هیدرولوگها را در جهت منطقهبندی مقدار نمایه ،تشکیل مناطق
همگن با منحنیهای رشد مشابه و برازش توابع توزیع بر
منحنیهای رشد منطقهای تجربی یاری میدهد.
دنباله بستههای نرمافزاری Rراجع به گشتاور...
بسته lomomRFA
نویسندهJ. R. M. Hosking :
این بسته دارای توابعی برای تحلیل فراوانی منطقهای منطبق
با روش Hoskingو )1997( Wallisدر « Regional Frequency
»Analysis: An Approach Based on L-momentsاست.
67
برآورد پارامترهای توزیع GPD
تابع درستنمایی GPDبه صورت زیر است.
) L(,)=f(y1,…,yk|,
تابع لگاریتم درستنمایی به صورت زیر است.
) l ( , ) log L ( ,
k
) k log (1 1 / ) log(1 y i /
i 1
اگر =0باشد .آنگاه
i
k
y
i 1
68
1
l ( ) k log
برآورد
بطور کلی برآوردها به دو دسته تقسیم میشوند.
برآورد نقطهای
برآورد فاصلهای
برآورد نقطهای همانطور که از اسمش پیداست ،نشان از
یک نقطه دارد .مانند:
X ; m
S ; s
در این نوع برآورد ،اثر تغییر نمونهها در میزان برآورد
مشهود نیست.
این نوع برآورد دارای احتمال متناظر نیست.
69
دانشگاه صنعت آب و برق
برآورد فاصلهای
برای این برآورد یک نوع فاصله در نظر گرفته میشود.
Lower # < population parameter < Upper #
به عنوان مثال:
P(Lower # < µ< Upper # )=1-α
احتمال متناظر این فاصله برابر با 1-است.
usually 90%, 95%, or 99%
)( = 10%), ( = 5%), ( = 1%
70
دانشگاه صنعت آب و برق
تکمیل و بسط بیشتر روش MLE
همانطور که مالحظه شد در حالتی که با پارامتر مواجه
بودیم و از منحنی ) (curveدرستنمایی صحبت شد.
اکنون اگر تعداد پارامترها بیش از یکی باشد راجع به
سطح ) (surfaceدرستنمایی صحبت میشود.
برای تعدد پارامترها الزم است که تعاریف و قضیههایی
را مرور نمود.
71
تعریف ماتریس واریانس-کوواریانس
اگر بردار X=(X1,…Xn)Tدر نظر باشد .ماتریس
مذکور به شرح زیر است.
1, n
n , n
1,1
i,j
j ,i
n ,1
که در آن i,i=var(Xi)و i,j=cov(Xi,Xj)که ij
است.
اگر متغیر فوق از تابع احتمال چند متغیره نرمال پیروی
کند بردار میانگین به صورت =(1,…,n)Tو
است
ماتریس واریانس-کوواریانس پارامتر
X
M.V
آنN (
دوم ,
)
T
72
تابع نرمال چند متغیره
تابع احتمال توام نرمال به صورت زیر تعریف ميشود.
1
f (x )
exp ( x ) ( x ) x
n
1
T
1
2
1/ 2
k /2
) (2
X
که در آن | |دترمینان ماتریس واریانس-کوواریانس
است و هر تابع حاشیهای آن خود تابع نرمال میباشد.
ماتریس اطالعات مورد انتظار expected
information matrixکه با ) IE(نشان میدهند
میزان مقدار مورد انتظار سطح لگاریتم درستنمایی را
نشان میدهد.
73
(Fisher) ماتریس اطالعات مورد انتظار
I E ( )
e
1,1
( )
e 1,d ( )
e d ,d ( )
e i , j ( )
e j , i ( )
e d ,1 ( )
2
e i , j ( ) E
l ( )
i j
n
l ( ) log L ( )
log f
i
( x i ; )
i 1
74
ماتریس اطالعات مشاهده )(Fisher
این ماتریس برای تخمین جمالت ماتریس IEبهکار
میرود زیرا =oعموما مجهول است .بنابر تعریف
ماتریس مذکور عبارت است از:
l ( )
) l (
2
1 d
) l (
) l (
2
i j
2
d
2
1
2
) l (
2
2
j i
) l (
2
d i
I O ( )
برای تابع چند متغیره ،معکوس ماتریس اطالعات معادل
ماتریس-کوواریانس است.
75
فاصله اطمینان
قضیه :اگر x1,…xnمقادیر محقق شده از یک توزیع
درستنمایی و
(دارای چند پارامتر) باشد و ) l(.تابع
ˆ
برآورد به روش MLEباشد ،آنگاه برای nبزرگ
داریم:
) ) ˆ M V N ( , I (
1
E
d
با استفاده از این قضیه میتوان فاصله اطمینان را برای
) o=(1,…,nبدست آورد .اگر معکوس ) IE(oبا
بزرگ داریم:
ˆ
برایNn(
i,jنشان دهیم ،آنگاه) ,
i ,i
i
i
اگر برای 1-αفاصله اطمینان در نظر
76
داریم:
گرفته
شودˆi ،
z
i ,i
2
دنباله فاصله اطمینان
ˆ
برآوردگر به روش
قضیه :برای اندازه بزرگ nاگر
MLEپارامتر ) o=(1,…,dبا ماتریس واریانس-
کوواریانس Vباشد ،آنگاه اگر ) =g(یک تابع
ˆ o=g(داریم:
N ( , o
اسکالر باشد MLE ،از )V)
که
V V
با
o
T
d
,
,
1
این قضیه به روش دلتا شهرت دارد.
77
R کد مربوط در
: برای روش دلتا داریم
library(POT)
data(ardieres)
ardieres <- clust(ardieres, 4, 10 / 365, clust.max =
TRUE)
fitted <- fitgpd(ardieres[,"obs"], 5, 'mle')
rbind(gpd.fishape(fitted),gpd.fiscale(fitted))
: نتیجه میشود که
conf.inf.shape conf.sup.shape
[1,] -0.01491947
0.5846731
[2,]
1.79820282
3.8749015
78
روش پروفیل درستنمایی
اگر پارامتر =(1,…,n)باشد ˆو برآورد آن به
روش MLEباشد و فرض کنید که ) l(x;لگاریتم تابع
درستنمایی باشد ،آنگاه
2
2 l ( x ; ˆ ) 2 l ( x ; ) n n
1
بنابراین اگر یک پارامتر داشته باشیم تابع حدی
خواهد بود.
دو ویژگی را میتوان برای پروفیل درستنمایی در نظر
گرفت.
2
غالبا روش دلتا فاصله اطمینان بزرگتری را ایجاد میکند.
فاصله اطمینان نامتقارن است (عدم قطعیت روی کران باال
بیشتر است).
79
R کدهای مربوط به
set.seed(6)
x <- rgpd(100, 0, 5, 0.3)
mle <- fitgpd(x, 0)
ic.delta <- gpd.firl(mle, p = 0.8, conf = 0.9)
ic.prof <- gpd.pfrl(mle, p = 0.8, conf = 0.9,
range = c(5, 15))
rbind(ic.delta, ic.prof)
. نتیجه به صورت زیر است
conf.inf conf.sup
ic.delta 7.326755 10.59160
ic.prof 7.474747 10.70707
80
تعریف تابع برازش توزیع پرتو عمومی fitgpd
)fitgpd(data, threshold, est = "mle", corr...
:data داده های مورد نظر است.
:threshold آستانه مورد نظر
:est روش برآورد پارامترها
: این آرگومان منطقی ،ماتریس همبستگی را نشان
میدهد.
81
تعریف تابع فاصله اطمینان confint
confint(object, parm, level = 0.95, ...,
)range, prob, prof = TRUE
:object شی است که قرار فاصله اطمینان آن تعیین
گردد.
:parm پارامتری که قرار است فاصله اطمینان برای آن
محاسبه گردد.
:level سطح فاصله اطمینان
:range دامنه است که حد باال و پایین را نشان میدهد.
:prob احتمال عدم تجاوز Non-Exceedance
:prof پروفیل درستنمایی که نامتقارن است .اگر مقدار
82آن FALSEباشد ،روش دلتا اعمال میگردد.
توابع فاصله اطمینان Fisher
83
روش فیشر برای تعیین فاصله اطمینان پارامتر shape
) gpd.fishape(object, conf
:objectشی است که قرار فاصله اطمینان آن تعیین گردد.
:confسطح اطمینان
روش فیشر برای تعیین فاصله اطمینان پارامتر scale
) gpd.fiscale(fitted, conf
:objectشی است که قرار فاصله اطمینان آن تعیین گردد.
:confسطح اطمینان
روش فیشر برای تعیین فاصله اطمینان سطح دروه بازگشت )return level (rl
) gpd.firl(fitted, prob, conf
:objectشی است که قرار فاصله اطمینان آن تعیین گردد.
:احتمال عدم تجاوز Non-Exceedance
:confسطح اطمینان
توابع فاصله اطمینان profile likelihood
روش پروفیل درستنمایی برای تعیین فاصله اطمینان پارامتر
shape
gpd.pfshape(object, range, xlab, ylab, conf ,
)nrang, vert.lines = TRUE, ...
:object شی است که قرار فاصله اطمینان آن تعیین گردد.
:range دامنه است که حد باال و پایین را نشان میدهد.
:ylab,xlab نام محور طولها و نام محور عرضها را
مشخص میکند.
:conf سطح اطمینان
:nrange دامنه است که حد باال و پایین را نشان میدهد.
:vert.lines ترسیم کننده خطوط عمودی
84
دنباله توابع فاصله اطمینان profile likelihood
روش پروفیل درستنمایی برای تعیین فاصله اطمینان پارامتر
scale
gpd.pfscale(object, range, xlab, ylab, conf ,
)nrang, vert.lines = TRUE, ...
:object شی است که قرار فاصله اطمینان آن تعیین گردد.
:range دامنه است که حد باال و پایین را نشان میدهد.
:ylab,xlab نام محور طولها و نام محور عرضها را
مشخص میکند.
:conf سطح اطمینان
:nrange دامنه است که حد باال و پایین را نشان میدهد.
:vert.lines ترسیم کننده خطوط عمودی
85
دنباله توابع فاصله اطمینان profile likelihood
روش پروفیل درستنمایی برای تعیین فاصله اطمینان پارامتر
return level
gpd.pfsrl(object, range, xlab, ylab, conf ,
)nrang, vert.lines = TRUE, ...
:object شی است که قرار فاصله اطمینان آن تعیین گردد.
:range دامنه است که حد باال و پایین را نشان میدهد.
:ylab,xlab نام محور طولها و نام محور عرضها را
مشخص میکند.
:conf سطح اطمینان
:nrange دامنه است که حد باال و پایین را نشان میدهد.
:vert.lines ترسیم کننده خطوط عمودی
86
profile likelihood وfisher کدهای مربوط به برآورد پارامترها و فاصله اطمینان آنها به روش
library(POT)
x <- rgpd(500, 0, 1, 0.25)
mle <- fitgpd(x, 0, corr=T)
mle$param[1:2]
mle$std.err
mle$exceed
mle$corr
par(mfrow=c(1,2))
profile_shape <- confint(mle, parm = "shape")
profile_scale <- confint(mle, parm = "scale")
rbind(profile_shape,profile_scale)
fisher_shape <- gpd.fishape(mle)
fisher_scale <- gpd.fiscale(mle)
rbind(fisher_shape,fisher_scale)
87
88
مناسب بودن مدل برازش یافته
با سه نمودار این مطلب را مورد بررسی قرار می دهیم.
نمودار احتمال-احتمال )(PP-Plot
نمودار چندک-چندک )(QQ-Plot
منحنی چگالی احتمال
)library(POT
)data(ardieres
ardieres <- clust(ardieres, 4, 10 / 365, clust.max
)= TRUE
)'fitted <- fitgpd(ardieres[, "obs"], 6, 'mle
)plot(fitted,which=1
برای بدست آوردن نمودار دوم و سوم باید پارامتر which=2و
which=3گردد .باالخره نمودار چهارم which=4
89
نمودار احتمال-احتمال )(PP-Plot
برای ترسیم نمودار به روش زیر عمل شده است.
) Ptheo ( x i ) F ( x i
rank ( x i ) 0.5
n
Pobs ( x i )
که در آن nتعداد نمونه ها است.
اکنون در یک دستگاه مختصات نقاط زیر را رسم
میکنیم.
( Pobs ( x i ), Ptheo ( x i )) : i 1, , n
برای برازش مناسب باید نقاط حول نیمساز ربع اول
90قرار گیرند.
PP-Plot نمودار
91
نمودار چندک-چندک )(QQ-Plot
برای ترسیم نمودار به روش زیر عمل شده است.
) x theo F ( p i
rank ( x i ) 0.5
n
pi
که در آن nتعداد نمونه ها است.
اکنون در یک دستگاه مختصات نقاط زیر را رسم
میکنیم ( x i , x th eo ( i )) : i 1, , n .
برای برازش مناسب باید نقاط حول نیمساز ربع اول
92قرار گیرند.
QQ-Plot نمودار
93
ardieres <- clust(ardieres, 4, 10 / 365,
clust.max = TRUE)
fitted <- fitgpd(ardieres[, "obs"], 6,
'mle')
plot(fitted,which=4)
94
نمودار بازگشتها مورد نظر
95