亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        最小廣義特征值在多元方差分析中的應(yīng)用探討

        2018-05-22 13:17:28江忠偉
        統(tǒng)計(jì)與決策 2018年9期
        關(guān)鍵詞:平方和正態(tài)分布組內(nèi)

        江忠偉

        (中國(guó)人民銀行南通市中心支行,江蘇 南通 226007)

        0 引言

        多元方差分析是一元方差分析的推廣,在選擇檢驗(yàn)統(tǒng)計(jì)量方面,通常的做法是:考慮到組內(nèi)差異是由隨機(jī)誤差造成的,組間差異可能是由隨機(jī)誤差和系統(tǒng)誤差共同引起的,與一元方差分析的基本思想相同。在一元方差分析中,若各個(gè)總體之間沒有顯著差異,則組間離差平方和與組內(nèi)離差平方和近似相等。可以證明組間離差平和與組內(nèi)離差平方和的比值服從F分布,給定顯著性水平后,就可以算出臨界值即得出拒絕域。與一元方差分析不同的是:多元統(tǒng)計(jì)分析需要將一元方差分析中的組間離差平方和、組內(nèi)離差平方和推廣為組間離差陣以及組內(nèi)離差陣。然后基于組間離差陣與組內(nèi)離差陣的比值構(gòu)建檢驗(yàn)統(tǒng)計(jì)量,可以證明該統(tǒng)計(jì)量為wilks統(tǒng)計(jì)量,給定顯著性水平后,就可以算出臨界值即得出拒絕域[1]。另外還有一些其他的檢驗(yàn)統(tǒng)計(jì)量,例如Hotelling跡檢驗(yàn)統(tǒng)計(jì)量[2,3]、Pil?lai-Bartlett準(zhǔn)則檢驗(yàn)統(tǒng)計(jì)量(Pillai-Bartlett criterion)[4,5]Roy最大特征值檢驗(yàn)統(tǒng)計(jì)量(Roy’s Largest Root)[6],具體表達(dá)形式見表1。

        表1 四種檢驗(yàn)統(tǒng)計(jì)量匯總

        通過推導(dǎo)證明,四個(gè)檢驗(yàn)統(tǒng)計(jì)量經(jīng)過適當(dāng)?shù)淖冃尉D(zhuǎn)化成服從F分布的檢驗(yàn)統(tǒng)計(jì)量[7]。在進(jìn)行多元方差分析時(shí)選擇哪個(gè)檢驗(yàn)統(tǒng)計(jì)量,是一個(gè)很有實(shí)際意義的問題。Stevens[7]對(duì)上述四個(gè)檢驗(yàn)統(tǒng)計(jì)量拒絕原假設(shè)能力進(jìn)行了比較,結(jié)果表明:在相同條件下,Roy最大特征值檢驗(yàn)統(tǒng)計(jì)量拒絕能力最強(qiáng)。Olson[8]對(duì)上述四種檢驗(yàn)統(tǒng)計(jì)量的檢驗(yàn)穩(wěn)健性進(jìn)行了比較,結(jié)果表明:通常,Pillai-Bartlett準(zhǔn)則檢驗(yàn)統(tǒng)計(jì)量的穩(wěn)健性好。

        綜上所述,四個(gè)檢驗(yàn)統(tǒng)計(jì)量經(jīng)過適當(dāng)?shù)淖冃尉梢孕纬梢粋€(gè)以F分布為漸近分布的隨機(jī)變量,據(jù)此可以在給定的顯著性水平下,設(shè)置一個(gè)小概率事件:當(dāng)原假設(shè)成立時(shí),檢驗(yàn)統(tǒng)計(jì)量的取值落入構(gòu)建的小概率事件中,則拒絕原假設(shè)。例如,利用wilks檢驗(yàn)統(tǒng)計(jì)量進(jìn)行檢驗(yàn)的思路為:首先利用似然比原則導(dǎo)出服從wilks分布的檢驗(yàn)統(tǒng)計(jì)量;由于對(duì)wilks檢驗(yàn)統(tǒng)計(jì)量不夠熟悉,通常將wilks檢驗(yàn)統(tǒng)計(jì)量轉(zhuǎn)換成F檢驗(yàn)統(tǒng)計(jì)量;最后結(jié)合一個(gè)給定的顯著性水平,就確定了拒絕域,即檢驗(yàn)法則。其三個(gè)檢驗(yàn)統(tǒng)計(jì)量也是按照這種思路:先利用樣本資料導(dǎo)出一個(gè)統(tǒng)計(jì)量,再將該檢驗(yàn)統(tǒng)計(jì)量轉(zhuǎn)換成F檢驗(yàn)統(tǒng)計(jì)量,最后結(jié)合一個(gè)給定的顯著性水平確定拒絕域。有一個(gè)很自然的想法是:能否先對(duì)樣本資料進(jìn)行變換,然后再根據(jù)變換后的樣本資料構(gòu)建F檢驗(yàn)統(tǒng)計(jì)量進(jìn)行方差分析?

        1 基本思路

        多元方差分析的主要任務(wù)是檢驗(yàn)因子的不同處理(類型變量)對(duì)不同處理下得到的樣本觀測(cè)值(數(shù)值變量)有無顯著影響,即分類自變量對(duì)數(shù)值因變量有無顯著影響。該模型可以表述為:設(shè)分類自變量有K個(gè)處理,可以將每個(gè)處理看成一個(gè)總體,則有總體:

        從這K個(gè)總體抽取如下樣本:

        其中是相互獨(dú)立的。

        檢驗(yàn):

        H0:至少有一組i≠j,使得μi≠μj,H1:μ1=…=μK,可以對(duì)m個(gè)總體中的所有樣品做同一變換即選擇一個(gè)p維行向量與所有的樣品進(jìn)行線性組合,顯然:若H0:至少有一組i≠j,使得μi≠μj成立,則選取任意一個(gè)p維行向量,必有H0:至少有一組i≠j,使得≠成立;反之也是如此。

        另一方面,由于服從p維多元正態(tài)分布的向量的分量的線性組合仍然服從正態(tài)分布,所以變換之后的樣品數(shù)據(jù)仍然服從正態(tài)分布。據(jù)此可以構(gòu)建F檢驗(yàn)統(tǒng)計(jì)量進(jìn)行一元方差分析。但F檢驗(yàn)統(tǒng)計(jì)量的取值是無法確定的,雖然樣本觀測(cè)值是已知的,但p維行向量是未知的。如何求出?假設(shè)檢驗(yàn)的目的是尋找證據(jù)支持本文的觀點(diǎn)。通常的做法是設(shè)置兩個(gè)對(duì)立事件,然后尋找一個(gè)特例拒絕與本文觀點(diǎn)對(duì)立的觀點(diǎn),這樣可以從一定置信水平上認(rèn)為本文觀點(diǎn)是正確的。因?yàn)榫芙^一個(gè)觀點(diǎn)只需要找到一個(gè)特例就行了,而接受一個(gè)觀點(diǎn)需要考慮所有的情況(通常是做不到的),因此只需尋找特例來拒絕原假設(shè)。利用矩陣的譜分解以及向量的線性表出等知識(shí),可以解出上述F檢驗(yàn)統(tǒng)計(jì)量的最小值以及相對(duì)應(yīng)l′的具體形式。如何利用這個(gè)極端值?一般的,對(duì)于假設(shè)檢驗(yàn)中的原假設(shè)H0,可以認(rèn)為H0是根據(jù)實(shí)際問題提出來的,往往是從過去經(jīng)驗(yàn)中總結(jié)出來的,沒有充分理由不能拒絕它。所以在多元方差分析中,當(dāng)原假設(shè)為:H0:至少有一組i≠j,使得μi≠μj,若原假設(shè)為真,即各個(gè)總體的均值向量有顯著差異,此時(shí)各水平的系統(tǒng)誤差不為零,此時(shí)F檢驗(yàn)統(tǒng)計(jì)量(為組間離差平方和與組內(nèi)離差平方和的比值)會(huì)很大。但若由樣本計(jì)算出的F檢驗(yàn)統(tǒng)計(jì)量的值小到可以將其看成一個(gè)小概率事件,則可以認(rèn)為原假設(shè)是不正確的,此時(shí)有較大把握拒絕原假設(shè)H0,接受備擇假設(shè)H1。

        2 依據(jù)樣本資料直接構(gòu)造F檢驗(yàn)統(tǒng)計(jì)量

        設(shè)分類自變量有K個(gè)處理,可以將每個(gè)處理看成一個(gè)子總體,則有總體:

        從這K個(gè)子總體抽取如下樣本:

        其中是 相 互 獨(dú) 立的。按照上文的內(nèi)容,選擇一個(gè)p維向量l′與所有樣品相乘,得出線性組合后的樣本:

        樣本數(shù)據(jù)經(jīng)過線性組合后均變成了一維數(shù)據(jù),由上文可知,檢驗(yàn)H0:至少有一組i≠j,使得μi≠μj與檢驗(yàn)H0:至少有一組i≠j,使得≠是等價(jià)的。這樣就將多元方差分析轉(zhuǎn)換為一元方差分析??梢詷?gòu)造F檢驗(yàn)統(tǒng)計(jì)量進(jìn)行一元方差分析。這里存在兩個(gè)問題:第一個(gè)問題是該樣本數(shù)據(jù)經(jīng)歷線性組合之后是否仍然服從正態(tài)分布;第二個(gè)問題是變換后的樣本數(shù)據(jù)的組間離差平方和與組內(nèi)離差平方和是否仍然獨(dú)立。接下來分別論證這兩個(gè)問題。

        2.1 樣本數(shù)據(jù)線性組合后正態(tài)性證明

        在一元正態(tài)分布中,若Z~N(0 ,1) ,則X=μ+σ Z~N(μ,σ2)。類似的在多元正態(tài)分布中,可以類似的定義多元正態(tài)分布。設(shè)相互獨(dú)立且有相同的分布N(0 ,1),μ為p維常數(shù)向量,A為p階常數(shù)矩陣,則稱:x=μ+的分布為多元正態(tài)分布,記作

        可以利用上述定義證明樣本數(shù)據(jù)進(jìn)行線性組合后仍然服從正態(tài)分布。具體過程如下:

        協(xié)差陣∑可以分解為:∑=

        則可以寫成μj+

        故得證。

        2.2 線性組合后的數(shù)據(jù)組間離差與組內(nèi)離差平方和獨(dú)立性證明

        由上知樣本數(shù)據(jù)進(jìn)行線性組合后仍然服從正態(tài)分布,可以計(jì)算出變換后的樣本數(shù)據(jù)的總離差平方和SST、組間離差平方SSB和組內(nèi)離差平方和SSE,經(jīng)過適當(dāng)變形之后總離差平方和SST、組間離差平方SSB和組內(nèi)離差平方和SSE均服從卡方分布,若組間離差平方SSB和組內(nèi)離差平方和SSE相互獨(dú)立,則可以構(gòu)造出F檢驗(yàn)統(tǒng)計(jì)量進(jìn)行方差分析。下面證明組間離差平方SSB和組內(nèi)離差平方和SSE相互獨(dú)立。

        變換后樣本數(shù)據(jù)的總離差平方和SST、組間離差平方SSB和組內(nèi)離差平方和SSE為:

        p維行向量l′為一個(gè)常數(shù)向量,要證明組間離差平方SSB和組內(nèi)離差平方和SSE之間相互獨(dú)立,即證明組間離差陣B和組內(nèi)離差陣E相互獨(dú)立。隨機(jī)矩陣的獨(dú)立性可以利用的科克朗(Cochran)定理來證明:設(shè)X~Nn×p(M,In?Σ ),C和D為n階對(duì)稱矩陣,X′CX與X′DX獨(dú)立,當(dāng)且僅當(dāng)CD=0。另外,若A是投影陣則I-A也是投影陣并且有A(I-A)=0成立。利用科克朗(Cochran)定理以及投影陣的性質(zhì),可以很方便地證明組間離差平方SSB和組內(nèi)離差平方和SSE之間相互獨(dú)立。具體證明過程如下:

        資料陣Y~Nn×p(M,In?Σ ),其中M的各行是各個(gè)子總體的均值向量的轉(zhuǎn)置按照各個(gè)子總體的觀測(cè)次數(shù)重復(fù)排列而成。

        可以將總離差陣改寫成:

        其中:

        可以驗(yàn)證:

        故C為投影陣且rank(C)=n-1;

        類似的有:

        也可以將組內(nèi)離差陣E寫成:E=Y′C*Y

        其中,C*=diag(C2,…,CK)

        顯然C*也是投影陣并且rank(C*)=rank(C1)+rank(C2)+…+rank(CK)=n-K;

        組間離差陣B可以改寫成:

        其中,

        顯然有,C**=C**′, (C**)2=C**,故C**是投影陣并且rank(C**)=trC**=trC+trC*=K-1。

        C,C*,C**均為投影陣,并且有C=C*+C**,所以C*C**=0,由科克朗(Cochran)定理知組間離差陣B和組內(nèi)離差陣E是相互獨(dú)立的,故組間離差平方SSB和組內(nèi)離差平方和SSE之間相互獨(dú)立。

        綜上所述,本文可以構(gòu)造出F檢驗(yàn)統(tǒng)計(jì)量:

        3 構(gòu)建檢驗(yàn)法則

        可以將原假設(shè)和備擇假設(shè)設(shè)為:

        H0:H0:存在μi≠μj,i≠j;H1:μ1=…=μK

        由上文知,可以將原假設(shè)和備擇假設(shè)改寫成:

        H0:存在l′μi≠l′μj,i≠j;H1:l′μ1= … =l′μK

        并且這兩組原假設(shè)和備擇假設(shè)的檢驗(yàn)結(jié)果是等價(jià)的。檢驗(yàn)統(tǒng)計(jì)量為:

        是一個(gè)已知分布的統(tǒng)計(jì)量,只需要給出顯著性水平α就可以確定拒絕域的臨界值Fα即得出檢驗(yàn)法則。

        該F檢驗(yàn)統(tǒng)計(jì)量與一般的F統(tǒng)計(jì)量有所不同,其中的p維行向量l′事先并不知道,所以無法計(jì)算出檢驗(yàn)統(tǒng)計(jì)量的具體數(shù)值。但考慮到檢驗(yàn)的初衷:拒絕與本文觀點(diǎn)對(duì)立的觀點(diǎn),從而證明本文的觀點(diǎn)是正確的。故只需要找到一個(gè)特例說明與本文觀點(diǎn)對(duì)立的觀點(diǎn)是錯(cuò)誤的。原假設(shè)H0:存在l′μi≠l′μj,i≠j成立時(shí),即系統(tǒng)誤差不為零。所以組間離差平方和與組內(nèi)離差平方和應(yīng)該相差很大。若將樣本觀測(cè)值帶入檢驗(yàn)統(tǒng)計(jì)量F,計(jì)算得出的結(jié)果很小,小到可以看成是一個(gè)小概率事件,則我們有充分的理由拒絕原假設(shè)。所以上述的假設(shè)檢驗(yàn)問題就轉(zhuǎn)化為已知樣本數(shù)據(jù)的條件下求解F檢驗(yàn)統(tǒng)計(jì)量的最小值,再與臨界值Fα(下分為數(shù))做出比較。F檢驗(yàn)統(tǒng)計(jì)量的最小值的計(jì)算過程如下:

        組內(nèi)離差陣組間離差陣顯然E、B為正定矩陣并且是對(duì)稱矩陣,檢驗(yàn)統(tǒng)計(jì)量F可以改寫成:

        其中是p階對(duì)稱矩陣,故其特征值是實(shí)數(shù);又因?yàn)闉檎ň仃?,故其特征值全部大于零?/p>

        由矩陣的譜分解知:

        其中λ1≥λ2≥…≥λp為B相對(duì)于E的廣義特征值,β1,β2,…,βp為B相對(duì)于E的廣義特征值λ1≥λ2≥…≥λp所對(duì)應(yīng)的標(biāo)準(zhǔn)化特征向量。β2,…,βp為一組線性無關(guān)的p維向量,對(duì)β2,…,βp做適當(dāng)變換后,可以將其看成p維向量空間中的一組標(biāo)準(zhǔn)正交基,該正交基仍然記作β2,…,βp。

        由向量的線性表出知:

        其中a2,…,ap為常數(shù)。

        將式(2)、式(3)帶入式(1)得:

        當(dāng)l=βp時(shí),等號(hào)成立。

        綜上所述,檢驗(yàn)法則為:當(dāng)時(shí),有充分理由拒絕原假設(shè),接受備擇假設(shè);當(dāng)時(shí),不拒絕原假設(shè)。

        4 利用投影思想進(jìn)行多元方差分析的優(yōu)點(diǎn)

        傳統(tǒng)的構(gòu)造檢驗(yàn)統(tǒng)計(jì)量的步驟為:先構(gòu)造出一個(gè)統(tǒng)計(jì)量,該統(tǒng)計(jì)量的分布是不為我們所熟悉的,為此一般的做法是將該統(tǒng)計(jì)量做適當(dāng)變換使得變換后的統(tǒng)計(jì)量的分布漸近服從一個(gè)我們熟悉的分布即F分布。這一過程通常計(jì)算量較大,并且理論性強(qiáng)不易理解。若直接從投影后的樣本資料出發(fā)構(gòu)建F分布。首先,從推導(dǎo)過程中可以發(fā)現(xiàn),所使用的都是基本的統(tǒng)計(jì)知識(shí)以及一些線性代數(shù)知識(shí),推導(dǎo)過程也十分簡(jiǎn)單,可以方便大家理解以及運(yùn)用該分析方法;其次,隨著計(jì)算機(jī)的普及以及儲(chǔ)存技術(shù)的發(fā)展,所研究的數(shù)據(jù)往往是海量、高維的數(shù)據(jù),這是挖掘數(shù)據(jù)中有價(jià)值信息的一個(gè)障礙,利用投影思想可以將高維度數(shù)據(jù)變換成低維度,這種思想的應(yīng)用無疑帶來了巨大的便利。

        5 模擬

        為了證實(shí)方法的正確性,分兩步進(jìn)行模擬。

        第一步利用R軟件產(chǎn)生9個(gè)子總體,每個(gè)子總體有20個(gè)樣品,這9個(gè)子總體的均值向量和協(xié)方差陣相同,所有樣品均為5維向量(見表2),其中均值向量和協(xié)方差陣是隨機(jī)選取的,在此基礎(chǔ)上進(jìn)行多元方差分析,驗(yàn)證檢驗(yàn)結(jié)果是否能夠拒絕原假設(shè)。

        表2 第一個(gè)子總體前十個(gè)樣品的5維向量

        利用計(jì)算出的組間離差陣相對(duì)于組內(nèi)離差陣最小廣義特征值為λp=0.0129,故F檢驗(yàn)統(tǒng)計(jì)量的取值為F=,該分位點(diǎn)對(duì)應(yīng)的p=0.02703,非常接近0,因此有充分理由拒絕原假設(shè)。

        第二步繼續(xù)利用R軟件產(chǎn)生9個(gè)子總體,每個(gè)子總體有20個(gè)樣品,與第一步不同的是,這9個(gè)子總體的均值向量不相同,所有樣品均為5維向量(見表3),其中均值向量和協(xié)差陣是隨機(jī)選取的,在此基礎(chǔ)上進(jìn)行多元方差分析,驗(yàn)證檢驗(yàn)結(jié)果是否為不能拒絕原假設(shè)。

        表3 第一個(gè)子總體前十個(gè)樣品的5維向量

        利用計(jì)算出的組間離差陣相對(duì)于組內(nèi)離差陣最小廣義特征值為λp=0.0228,故F檢驗(yàn)統(tǒng)計(jì)量的取值為F=,該分位點(diǎn)對(duì)應(yīng)的p=0.136,不是一個(gè)非常小的數(shù)值,因此沒有充分理由拒絕原假設(shè)。

        6 結(jié)論

        本文首先利用投影思想構(gòu)建的F檢驗(yàn)統(tǒng)計(jì)量在模擬試驗(yàn)中,當(dāng)各個(gè)子總體均值向量之間不存在差異時(shí),檢驗(yàn)結(jié)果拒絕原假設(shè),接受備擇假設(shè);當(dāng)各個(gè)子總體均值向量之間確實(shí)存在差異,F(xiàn)檢驗(yàn)統(tǒng)計(jì)量的取值不能夠拒絕原假設(shè),故可以達(dá)到多元方差分析的目的。在假設(shè)檢驗(yàn)過程中,當(dāng)沒有充分理由拒絕原假設(shè)時(shí),這時(shí)很多人便認(rèn)為原假設(shè)是正確的。贊同這個(gè)觀點(diǎn)的人并沒有考慮原假設(shè)錯(cuò)誤但檢驗(yàn)統(tǒng)計(jì)量取值沒有落入拒絕域中的概率的大?。醇{偽的概率),若原假設(shè)錯(cuò)誤時(shí)建議統(tǒng)計(jì)量取值沒有落入拒絕域中的概率很大,這時(shí)認(rèn)為原假設(shè)是正確的顯然是不可信的。此時(shí)可以認(rèn)為檢驗(yàn)工作并沒有取得實(shí)質(zhì)進(jìn)展。如何有效克服這個(gè)問題有待更進(jìn)一步的探討。

        參考文獻(xiàn):

        [1] Finch H.Comparison of the Performance of Nonparametric and Para?metric MANOVA Test Statistics When Assumptions Are Violated[J].Methodology,2005,1(1).

        [2] Kapstad H,Hanestad B R,Langeland N,et al.Cutpoints for Mild,Moderate and Severe Pain in Patients With Osteoarthritis of the Hip or Knee Ready for Joint Replacement Surgery[J].BMC Musculoskele?tal Disorders,2008,9(1).

        [3] Hatlen M A,Arora K,Vacic V,et al.Integrative Genetic Analysis of Mouse and Human AML Identifies Cooperating Disease Alleles[J].The Journal of Experimental Medicine,2016,213(1).

        [4] Ullah I,Jones B.Regularised Manova for High-Dimensional Data[J].Australian&New Zealand Journal of Statistics,2015,57(3).

        [5] Chiani M.Distribution of the Largest Root of a Matrix for Roy’s Test in Multivariate Analysis of Variance[J].Journal of Multivariate Analy?sis,2016,(143).

        [6] Haase R F,Ellis M V.Multivariate Analysis of Variance[J].Journal of Counseling Psychology,1987,34(4).

        [7] Stevens J P.Power of the Multivariate Analysis of Variance Tests[J].Psychological Bulletin,1980,88(3).

        [8] Olson C L.On Choosing a Test Statistic in Multivariate Analysis of Variance[J].Psychological Bulletin,1976,83(4).

        猜你喜歡
        平方和正態(tài)分布組內(nèi)
        用心說題 提高效率 培養(yǎng)能力
        費(fèi)馬—?dú)W拉兩平方和定理
        利用平方和方法證明不等式賽題
        基于對(duì)數(shù)正態(tài)分布的出行時(shí)長(zhǎng)可靠性計(jì)算
        勾股定理的擴(kuò)展
        正態(tài)分布及其應(yīng)用
        關(guān)于四奇數(shù)平方和問題
        正態(tài)分布題型剖析
        χ2分布、t 分布、F 分布與正態(tài)分布間的關(guān)系
        合作學(xué)習(xí)組內(nèi)交流討論時(shí)間的遵循原則
        日韩女优av一区二区| 91精品欧美综合在线观看| av网址大全在线播放| 婚外情长久的相处之道| 国产精品久久久久久久| 中文字幕亚洲无线码| 亚洲Av无码专区尤物| 蜜乳一区二区三区亚洲国产| 色欲人妻aaaaaaa无码| 韩国精品一区二区三区无码视频| 国产一区二区三区精品久久呦| 久久久噜噜噜噜久久熟女m| 亚洲日韩成人无码| 少妇无码一区二区三区| 亚洲色www无码| 日本一级片一区二区三区| 无码国产色欲xxxx视频| chinese国产乱在线观看| 久久视频在线视频精品| 日本一区二区三区人妻| 无码吃奶揉捏奶头高潮视频| 国产成人国产在线观看| 亚洲综合中文一区二区| 99国产精品99久久久久久| 男人边吻奶边挵进去视频| 国产在线天堂av| 一本色道久久亚洲精品| 久久久久亚洲av成人无码| 日韩高清毛片| 伊人狼人影院在线视频| 国产亚洲欧洲aⅴ综合一区| 色噜噜狠狠色综合成人网| 国产亚洲无码1024| 亚洲av成人av三上悠亚| 精品亚洲一区二区三区在线观看 | h动漫尤物视频| 人妖一区二区三区在线| 精品久久久久成人码免费动漫| 国产精品自产拍在线观看免费| 久久久噜噜噜久久熟女| 中国人妻与老外黑人|