999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

醫(yī)療統(tǒng)計(jì)周期性預(yù)測問題的比較研究

2016-12-09 07:51:38李望晨于貞杰王在翔張利平
統(tǒng)計(jì)與決策 2016年19期
關(guān)鍵詞:模型

李望晨,于貞杰,王在翔,張利平

(濰坊醫(yī)學(xué)院“健康山東”重大社會風(fēng)險(xiǎn)預(yù)測與治理協(xié)同創(chuàng)新中心,山東濰坊261053)

醫(yī)療統(tǒng)計(jì)周期性預(yù)測問題的比較研究

李望晨,于貞杰,王在翔,張利平

(濰坊醫(yī)學(xué)院“健康山東”重大社會風(fēng)險(xiǎn)預(yù)測與治理協(xié)同創(chuàng)新中心,山東濰坊261053)

文章立足醫(yī)療統(tǒng)計(jì)領(lǐng)域周期性預(yù)測問題,以門診人次為例進(jìn)行建模方法比較,為類似研究提供方法參考。采用方式為:方式1:門診人次季度數(shù)據(jù)由X11法時(shí)序分解后提取季節(jié)指數(shù)消除季節(jié)波動,用ARIMA法和拋物線擬合外推。方式2:周期差分提取季節(jié)信息,二階差分提取長期趨勢,直接用季節(jié)ARIMA法建模。季節(jié)波動和短期相關(guān)間交互影響不清,試用簡單季節(jié)或乘積季節(jié)模型。得出結(jié)論:方式1對周期時(shí)序資料分解后,由組合方法依次直觀反映;方式2直接用季節(jié)模型擬合預(yù)測,乘積季節(jié)模型性能并非有替代性,簡單季節(jié)模型更優(yōu)。算例為載體用幾類方法比較研究,對于醫(yī)院統(tǒng)計(jì)領(lǐng)域周期性預(yù)測問題有借鑒意義。

時(shí)間序列;組合模型;X11法;ARIMA;乘積季節(jié)模型

0 引言

醫(yī)療統(tǒng)計(jì)工作領(lǐng)域經(jīng)常見時(shí)間序列預(yù)測問題,如傳染病或慢性病發(fā)病數(shù)、發(fā)病率、門診或住院人次、醫(yī)療營業(yè)收入、衛(wèi)生費(fèi)用支出等。這些事物受到社會、經(jīng)濟(jì)、環(huán)境等復(fù)雜因素的影響,對其解讀分析也是衛(wèi)生政策制定者、公共衛(wèi)生干預(yù)者、衛(wèi)生行政管理者關(guān)心的問題。

由于影響因素復(fù)雜性,傳統(tǒng)回歸模型結(jié)構(gòu)模式確定、參數(shù)選擇困難,指標(biāo)難以量化。時(shí)序資料本身隨時(shí)間延伸表現(xiàn)規(guī)律性,適合于時(shí)間序列分析。時(shí)間序列方法在數(shù)學(xué)、統(tǒng)計(jì)學(xué)、計(jì)量經(jīng)濟(jì)和信息科學(xué)推動下,理論體系和應(yīng)用實(shí)現(xiàn)較為成熟[1,2],近年來在醫(yī)療衛(wèi)生領(lǐng)域得到推廣應(yīng)用,如ARIMA法、曲線擬合法、灰色系統(tǒng)法、數(shù)據(jù)挖據(jù)技術(shù),其中ARIMA法有代表性,對于隨機(jī)時(shí)序資料擬合較好,對于趨勢性、周期性特點(diǎn)資料也有很強(qiáng)適應(yīng)性。醫(yī)療統(tǒng)計(jì)領(lǐng)域常見周期性數(shù)據(jù)資料[3](如門診人次、收入、費(fèi)用),多采集為月度或季度資料,觀察期長、數(shù)據(jù)豐富、周期明顯、長期平緩,這類事物規(guī)律性好于傳染病疫情時(shí)序資料。對此有必要引入適當(dāng)方法建立模型,借助算例實(shí)現(xiàn)為類似研究提供方法借鑒。

第一種方式,鑒于月度或季度周期變化時(shí)序資料周期性、趨勢性較明顯,可用時(shí)間序列分解法提取信息,X-11過程法可用于提取季節(jié)指數(shù)[4],趨勢因素可用曲線擬合法或ARIMA法。時(shí)序分解后從兩種視角建立擬合模型,外推預(yù)測及合并信息后再預(yù)測未來情況。第二種方式,傳染病預(yù)測為當(dāng)前熱點(diǎn)問題,乘積ARIMA法對傳染病月度數(shù)據(jù)表現(xiàn)應(yīng)用代表性。對于季節(jié)性變化資料,首先用步長差分方法提取季節(jié)信息,低階差分提取長期趨勢信息,季節(jié)方法與ARIMA模型結(jié)合建模型,必要時(shí)比較簡單季節(jié)模型與乘積季節(jié)模型研究效果。

1 對象資料

以醫(yī)療統(tǒng)計(jì)領(lǐng)域時(shí)序資料為立足點(diǎn)建立模型和進(jìn)行實(shí)證研究,方法學(xué)視角下對于門診或住院人次、醫(yī)療收入等周期性指標(biāo)有普適借鑒意義。某二級甲等醫(yī)院2004—2014年按季度統(tǒng)計(jì)門診人次數(shù)據(jù)為例,以2004—2013年建立擬合模型,對2014年進(jìn)行預(yù)測驗(yàn)證。門診人次數(shù)據(jù)隨季度變化呈周期波動,年度間有平緩增長特點(diǎn),見表1和圖1所示。

表1 某醫(yī)院門診人次數(shù)據(jù)

2 建模方式一

原始序列有明顯趨勢變化特點(diǎn),對2004—2013年時(shí)序分解,提取季節(jié)指數(shù)后剩余序列有趨勢增長及隨機(jī)變化特點(diǎn),用曲線擬合法或ARIMA法擬合趨勢規(guī)律并外推預(yù)測未來[5],與季節(jié)指數(shù)相乘后還原為2014年情況。組合方法在SAS軟件環(huán)境下計(jì)算實(shí)現(xiàn)方便。

2.1季節(jié)指數(shù)提取

由于門診人次數(shù)據(jù)有明顯季節(jié)變化特點(diǎn),可經(jīng)時(shí)間序列分解提取長期趨勢、周期性季節(jié)指數(shù)、不規(guī)則變化部分等。X-11法由美國國情調(diào)查局建立,由于需要用到11次移動平均而得名,可據(jù)此對門診人次數(shù)據(jù)提取季節(jié)指數(shù): 92.882%、102.616%、99.544%、104.963%,消除季節(jié)因素后長期趨勢見圖2、剩余隨機(jī)序列見圖3所示。

圖1 原始數(shù)據(jù)

圖2 長期趨勢

圖3 不規(guī)則波動

2.2ARIMA法

ARIMA表達(dá)式Φ(B)▽dxt=Θ(B)εt,其中φi為自回歸系數(shù),θj為移動平均系數(shù)。Φ(B)=(1-φ1B-…-φpBp)為自回歸多項(xiàng)式,Θ(B)=(1-θ1B-…-θqBq)為移動平均多項(xiàng)式, Bkxt=xt-k為k步延遲算子,▽d=(1-B)d為d階差分算子。

根據(jù)數(shù)據(jù)資料確定模型結(jié)構(gòu)、識別參數(shù)、檢驗(yàn)?zāi)P惋@著性,否則要不斷調(diào)整優(yōu)化模型。趨勢變化序列要經(jīng)低階差分后轉(zhuǎn)化為平穩(wěn)序列。由AIC、SBC最小準(zhǔn)則選擇最優(yōu)模型,由最小二乘法識別參數(shù)。由殘差純隨機(jī)性檢驗(yàn)后,由模型進(jìn)行短期外推預(yù)測。

正試期第2、第4周采用全糞收集方式采集糞樣,每組用蛇皮袋固定于隨機(jī)選取的4只試驗(yàn)羊尾后,每日更換兩次并稱重,將飼料樣與糞樣帶回實(shí)驗(yàn)室進(jìn)行常規(guī)分析計(jì)算養(yǎng)分表觀消化率。正試期第2、第4周每組隨機(jī)選取4只試驗(yàn)羊,晨飼后3 h利用負(fù)壓原理經(jīng)口腔采集瘤胃液各30 ml,四層紗布過濾后,采用上海雷磁pHS-3C型精密pH計(jì)立即測定pH值,之后將濾液分裝在3個(gè)10 ml滅菌離心管中,置于-20℃下保存待測NH3-N濃度和MCP含量。

門診人次序列消除季節(jié)指數(shù),長期趨勢序列用ARIMA法擬合。趨勢序列經(jīng)二階差分消除趨勢。在p,q≤6中自動尋優(yōu)定階,BIC(1,0)=11.757最小,認(rèn)為AR(1)模型最合適,這也符合簡化原則。條件最小二乘法識別參數(shù),說明參數(shù)顯著(P=0.0001),(1+0.57965B)(1-B)2xt=εt為表達(dá)式。經(jīng)模型顯著性檢驗(yàn),殘差為純隨機(jī)序列,延遲階數(shù)6,12,18,24時(shí),P=0.7355,0.7495,0.8347,0.0672>0.05,模型對序列擬合顯著有效。對2014年各季度趨勢外推預(yù)測,得到70202,71347,72614,73810,由季節(jié)指數(shù)合并計(jì)算2014年門診人次65205,73213,72283,77473,相對誤差絕對值為0.59%,1.60%,1.21%,1.81%。

2.3曲線擬合法

門診人次序列消除季節(jié)因素后有平穩(wěn)長期趨勢特點(diǎn),可用曲線擬合法建立模型,如二次或三次拋物線、指數(shù)曲線、Logistic曲線等。經(jīng)試選,本例適合三次拋物線反映趨勢,對近期數(shù)據(jù)變化規(guī)律擬合更好,表達(dá)式為xt=52710+457×t-23.344×t2+0.555×t3。決定系數(shù)為R2=0.982,F(xiàn)檢驗(yàn)P<0.001,模型擬合不錯(cuò),對2014年各季度趨勢外推預(yù)測,得到70424,71808,73286,74860,由季節(jié)指數(shù)合并計(jì)算2014年門診人次65411,73687,72952,78576,相對誤差絕對值為0.91%,2.26%,2.15%,3.26%。

門診人次月度數(shù)據(jù)消除季節(jié)波動后,ARIMA法和趨勢外推法均可用于建立模型。ARIMA法為隨機(jī)時(shí)序方法,將歷史序列和殘差序列納入線性模型。不易對事物規(guī)律直觀解釋,優(yōu)點(diǎn)是擬合性能和短期預(yù)測精度高。趨勢外推曲線能解釋事物隨時(shí)間變化規(guī)律,須有曲線變化特點(diǎn),序列波動明顯存在時(shí)不太適合。經(jīng)計(jì)算比較,ARIMA效果優(yōu)于多項(xiàng)式曲線。趨勢擬合與外推效果比較見圖4—圖5,組合模型對原始序列擬合及外推效果見圖6所示。

圖4 ARIMA模型

圖5 三次拋物線模型

圖6 三種組合模型擬合外推效果比較

3 建模方式二

考慮門診人次數(shù)據(jù)隨季度呈周期波動變化、隨年度長期平緩趨勢增長特點(diǎn),借鑒當(dāng)前傳染病預(yù)測文獻(xiàn)中常見季節(jié)ARIMA模型[6~7],用以建立擬合模型及外推預(yù)測。

3.1簡單季節(jié)AR IMA模型

由于門診人次序列不平穩(wěn),用周期步長差分消除季節(jié)波動,用低階差分消除平緩增長趨勢,將數(shù)據(jù)轉(zhuǎn)化為平穩(wěn)序列。模型結(jié)構(gòu)如下:

門診人次序列隨季節(jié)周期波動、隨年度有長期趨勢,季節(jié)變化用4步差分,經(jīng)一階差分或二階差分消除趨勢特點(diǎn),過度差分會損失信息,差分序列若平穩(wěn)化則可以建立模型。經(jīng)計(jì)算分析,以2階差分消除趨勢、以4步差分消除季節(jié)波動,經(jīng)白噪聲檢驗(yàn)得知,延遲6階時(shí),P=0.1452>0.05,該序列可以建模分析。觀察ACF、PACF是否落入兩倍標(biāo)準(zhǔn)差內(nèi),調(diào)試階數(shù)、確定結(jié)構(gòu)和識別參數(shù)。模型結(jié)構(gòu)以簡單為原則,發(fā)現(xiàn)ACF一階以后、PACF一階及四階時(shí)落入兩倍標(biāo)準(zhǔn)差范圍并繼續(xù)拖尾,考慮低階疏系數(shù)模型,模型ARIMA((4),1,1)表達(dá)式如下:

(1+0.38669B4)(1-B4)(1-B)2xt=(1+0.90475)εt。經(jīng)延遲6,12,18,24階時(shí),P=0.4509,0.8271,0.9703,0.5971,說明殘差為白噪聲序列。對2014年門診人次各季度外推預(yù)測,得到65127,72420,71452,76067;計(jì)算相對誤差絕對值為0.47%,0.50%,0.05%,0.03%,預(yù)測效果非常好。

3.2乘積季節(jié)AR IMA模型

如果季節(jié)波動與長期趨勢間有復(fù)雜交互影響,要用乘積季節(jié)ARIMA模型。該方法在傳染病發(fā)病率、發(fā)病數(shù)預(yù)測中是當(dāng)前熱點(diǎn),可能是傳染病月度或季度資料特點(diǎn)決定的。對于本例門診人次預(yù)測來說,簡單季節(jié)模型經(jīng)計(jì)算預(yù)測效果不錯(cuò),可認(rèn)為季節(jié)波動與短期相關(guān)間并無復(fù)雜關(guān)系,以下仍用乘積季節(jié)模型驗(yàn)證和進(jìn)行比較。

以ARMA(p,q)提取短期相關(guān)性,以ARMA(P,Q)提取季節(jié)相關(guān)性。假設(shè)二者存在交互乘積關(guān)系,構(gòu)造為乘積模型: ARIMA(p,d,q)×(P,D,Q)S:Φ(B)ΦS(B)▽D▽dxt=Θ(B)ΘS(B)εt;

其中Θ(B)=(1-θ1B-…-θqBq),

前面得知長期趨勢用2階差分,季度周期用四步差分。調(diào)試階數(shù)及識別參數(shù),模型結(jié)構(gòu)ARIMA(1,2,0)×(1,1,0)4,表達(dá)式為(1-0.18884B)(1-0.32985B4)(1-B4)(1-B)2xt。=εt殘差序列延遲6,12,18,24時(shí),P=0.1178,0.5641,0.6726,0.1108>0.05,認(rèn)為模型顯著。對2014年門診人次各季度外推預(yù)測,得到65061,72547,71978,76881,計(jì)算相對誤差絕對值0.37%, 0.68%,0.78%,1.04%。然后根據(jù)ACF和PACF特點(diǎn),繼續(xù)選擇可能模型結(jié)構(gòu)如ARIMA(1,2,1)×(1,1,0)4,表達(dá)式為(1+0.40425B)(1-0.25856B)(1-B4)(1-B)2xt=(1+0.9655B)εt。延遲6,12,18,24時(shí),殘差序列P=0.1878,0.6384,0.7670,0.2152>0.05,認(rèn)為模型顯著。對2014年門診人次各季度外推預(yù)測,得到64988,72311,71860,76605,計(jì)算相對誤差絕對值0.25%, 0.35%,0.62%,0.67%。

乘積季節(jié)ARIMA模型預(yù)測精度稍微改善,又反復(fù)試取了多種模型結(jié)構(gòu),但相對簡單季節(jié)ARIMA模型并未明顯改善預(yù)測精度,認(rèn)為季節(jié)波動和短期相關(guān)無交互關(guān)系,比較見圖7所示。

圖7 兩種季節(jié)模型擬合外推效果

以上從兩種方式、四種模型對門診人次季度數(shù)據(jù)進(jìn)行預(yù)測。第一種方式:X11過程法與ARIMA、多項(xiàng)式曲線設(shè)計(jì)組合模型。第二種方式:假設(shè)季節(jié)波動和短期相關(guān)性是否有交互關(guān)系,構(gòu)造簡單季節(jié)ARIMA模型、乘積季節(jié)ARIMA模型。2004—2013年季度數(shù)據(jù)用于確定模型結(jié)構(gòu)、識別參數(shù),對2014年情況預(yù)測情況比較見表2所示。

表2 模型預(yù)測精度比較

4 討論

門診、收入等醫(yī)療統(tǒng)計(jì)工作領(lǐng)域多數(shù)指標(biāo)為例,月度、季度數(shù)據(jù)隨時(shí)間變化往往有規(guī)律性。影響因素復(fù)雜而無法確定模型結(jié)構(gòu)或識別參數(shù),可以嘗試由時(shí)序數(shù)據(jù)由時(shí)間序列分析方法建模,時(shí)序數(shù)據(jù)可分解為長期趨勢、季節(jié)波動和隨機(jī)無序部分,由數(shù)據(jù)資料計(jì)算季節(jié)指數(shù),消除周期波動后分析長期趨勢,也可直接應(yīng)用季節(jié)ARIMA模型。

一種方式是,鑒于季節(jié)指數(shù)較明顯,以X11法作為時(shí)間序列分解法提取信息,原始序列不再有季節(jié)波動特點(diǎn),考慮這類序列有平滑且類似曲線變化特點(diǎn),以某類曲線擬合時(shí)序規(guī)律,也可由ARIMA法直接建模。經(jīng)外推預(yù)測后將預(yù)測值與季節(jié)指數(shù)合成未來預(yù)測值。另一種方式是,序列季節(jié)波動與長期趨勢不給予分解提取,而是直接建立季節(jié)ARIMA模型,若季節(jié)波動與短期相關(guān)交互不明顯則用簡單季節(jié)模型,若交互影響存在則用乘積季節(jié)模型。

以某醫(yī)院門診人次預(yù)測為例,季度數(shù)據(jù)周期性變化有規(guī)律,歷史數(shù)據(jù)中季節(jié)波動且平穩(wěn)增長,觀測期長、數(shù)據(jù)連貫、規(guī)律性強(qiáng),明顯有季節(jié)波動和平穩(wěn)延續(xù)特點(diǎn),后期增幅明顯,長期較平緩且無轉(zhuǎn)折因素影響。采用以上幾種模型演示實(shí)施流程,為醫(yī)療統(tǒng)計(jì)工作提供方法借鑒。兩種方式又分別給出兩種算法,比較模型原理特點(diǎn)、擬合性能和預(yù)測精度。以2004—2013年數(shù)據(jù)擬合模型,預(yù)測2014年情況并與真實(shí)值計(jì)算相對誤差,X11-ARIMA法1%左右,X11-三次多項(xiàng)式2%左右,簡單季節(jié)ARIMA0.5%以下,乘積季節(jié)ARIMA0.5%左右。

第一種方式為組合建模方法,適用于季節(jié)指數(shù)較明確情況,消除季節(jié)波動以后的長期趨勢規(guī)律,ARIMA法好于曲線擬合法。第二種方式直接用季節(jié)ARIMA模型,經(jīng)算例驗(yàn)證表明季節(jié)ARIMA模型擬合性能、預(yù)測精度均好于第一種方式。預(yù)測領(lǐng)域多數(shù)文獻(xiàn)表明,月度或季度數(shù)據(jù)基本上用乘積季節(jié)ARIMA模型,這種思路未必千篇一律,經(jīng)本例驗(yàn)證,門診人次數(shù)據(jù)中季節(jié)波動和短期相關(guān)性無交互影響,簡單季節(jié)ARIMA模型甚至優(yōu)于乘積季節(jié)ARIMA模型,具體資料情況下可以進(jìn)行試選優(yōu)化,以便改善模型對資料擬合及預(yù)測適應(yīng)性。

[1]王燕.應(yīng)用時(shí)間序列分析[M].北京:中國人民大學(xué)出版社,2012.

[2]徐國祥.統(tǒng)計(jì)預(yù)測與決策[M].上海:上海財(cái)經(jīng)大學(xué)出版社,2012.

[3]耿娟.ARIMA模型在醫(yī)院門診量預(yù)測中的應(yīng)用[J].中國衛(wèi)生統(tǒng)計(jì), 2014,31(8).

[4]申銅倩,劉文東,胡建利等.X11-ARIMA過程在痢疾疫情預(yù)測中的應(yīng)用研究[J].中國衛(wèi)生統(tǒng)計(jì),2014,31(3).

[5]李望晨.基于增長特征法與ARIMA的人均衛(wèi)生事業(yè)費(fèi)趨勢預(yù)測比較研究[J].中國衛(wèi)生統(tǒng)計(jì),2014,31(3).

[6]陸波,閔思韜,閔紅星等.應(yīng)用ARIMA模型預(yù)測麻疹發(fā)病率的可行性研究[J].中國衛(wèi)生統(tǒng)計(jì),2014,32(1).

[7]張愛紅,周培,申銅倩等.乘積季節(jié)ARIMA模型在食源性疾病預(yù)測中的應(yīng)用[J].中國衛(wèi)生統(tǒng)計(jì),2014,31(3).

(責(zé)任編輯/浩天)

C913.4

A

1002-6487(2016)19-0078-03

教育部人文社科基金資助項(xiàng)目(15YJCZH087;14YJAZH101;14YJA630098);山東自然科學(xué)基金資助項(xiàng)目(ZR2015HL101);山東統(tǒng)計(jì)科研項(xiàng)目(KT16230;KT16231)

李望晨(1980—),男,山東濰坊人,碩士,副教授,研究方向:衛(wèi)生管理統(tǒng)計(jì)。于貞杰(1971—),女,山東濰坊人,博士,教授,研究方向:衛(wèi)生事業(yè)管理。王在翔(1963—),男,山東濰坊人,碩士,教授,研究方向:衛(wèi)生管理統(tǒng)計(jì)。(通訊作者)張利平(1980—),男,山東濰坊人,博士研究生,講師,研究方向:職業(yè)衛(wèi)生評價(jià)。

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機(jī)模型
提煉模型 突破難點(diǎn)
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 97在线国产视频| 亚洲精品无码AⅤ片青青在线观看| 国产黄网站在线观看| 亚洲国模精品一区| 亚洲天堂视频在线观看| 中国一级特黄大片在线观看| 精品福利视频导航| jizz国产在线| 久久人搡人人玩人妻精品| 无码不卡的中文字幕视频| 无码国产伊人| 精品视频在线观看你懂的一区| 青青青国产精品国产精品美女| 亚洲aaa视频| AV老司机AV天堂| 狠狠色婷婷丁香综合久久韩国| 欧美视频二区| 免费无码又爽又黄又刺激网站| 成人免费黄色小视频| 国产精品理论片| 无码aaa视频| 色综合久久88色综合天天提莫| 免费看a毛片| 免费无码网站| 久久视精品| 色婷婷综合在线| 精品午夜国产福利观看| 免费啪啪网址| 色天堂无毒不卡| 成人噜噜噜视频在线观看| 激情综合婷婷丁香五月尤物| 精品国产自在在线在线观看| 99re经典视频在线| 午夜影院a级片| 国产亚洲精品无码专| 在线精品视频成人网| 国产成人8x视频一区二区| 中文无码毛片又爽又刺激| 日韩毛片在线播放| 谁有在线观看日韩亚洲最新视频| 免费jizz在线播放| 干中文字幕| 高清无码一本到东京热 | 欧美日韩一区二区在线播放 | 日本a∨在线观看| 狠狠色噜噜狠狠狠狠奇米777| 71pao成人国产永久免费视频| av无码久久精品| 高清无码手机在线观看| 亚洲色图欧美| 日韩一区精品视频一区二区| 欧亚日韩Av| 大学生久久香蕉国产线观看| 欧美精品色视频| 欧美午夜性视频| 国产91精品久久| 色综合久久综合网| 都市激情亚洲综合久久| 在线一级毛片| 国产女人在线观看| 99久久精品视香蕉蕉| 视频在线观看一区二区| 精品夜恋影院亚洲欧洲| 99久久精品免费观看国产| 久久精品电影| 欧美日韩一区二区三| 国产欧美日韩va另类在线播放| 久久人午夜亚洲精品无码区| 无码AV高清毛片中国一级毛片| 精品国产黑色丝袜高跟鞋| 亚洲精品无码在线播放网站| 亚洲天堂区| 亚洲男人的天堂在线| 伊在人亚洲香蕉精品播放| 国产导航在线| 欧美一级大片在线观看| 青草视频免费在线观看| 中文字幕精品一区二区三区视频| www.狠狠| 国产簧片免费在线播放| 波多野一区| 色综合五月婷婷|