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

ARIMA模型在成都市成華區(qū)狂犬病暴露監(jiān)測(cè)數(shù)據(jù)分析中的應(yīng)用

2017-01-09 13:43:20強(qiáng)
關(guān)鍵詞:趨勢(shì)模型

楊 靜 張 強(qiáng)

ARIMA模型在成都市成華區(qū)狂犬病暴露監(jiān)測(cè)數(shù)據(jù)分析中的應(yīng)用

楊 靜1,2張 強(qiáng)1△

目的了解成都市成華區(qū)狂犬病暴露監(jiān)測(cè)數(shù)據(jù)的基本情況和特征;利用狂犬病暴露數(shù)據(jù)建立ARIMA模型,對(duì)2016年狂犬病暴露數(shù)據(jù)進(jìn)行預(yù)測(cè),為人用狂犬病疫苗、免疫球蛋白需求計(jì)劃,犬傷處置報(bào)銷(xiāo)費(fèi)用預(yù)算的制定提供參考依據(jù)。方法應(yīng)用SPSS19.0對(duì)成都市成華區(qū)2009-2014年逐月狂犬病暴露數(shù)據(jù)構(gòu)建ARIMA模型,以2015年狂犬病暴露數(shù)據(jù)為驗(yàn)證樣本,驗(yàn)證模型的預(yù)測(cè)效果,并預(yù)測(cè)2016年逐月狂犬病暴露數(shù)據(jù)。結(jié)果最優(yōu)模型ARIMA(1,0,0)(1,1,0)12能較好地?cái)M合既往時(shí)間段的狂犬病暴露數(shù)據(jù)序列,擬合值與實(shí)際值基本保持一致的曲線趨勢(shì)。2015年1月-12月檢驗(yàn)樣本的預(yù)測(cè)結(jié)果顯示,暴露實(shí)際值均在預(yù)測(cè)可信區(qū)間內(nèi),均方根誤差為28.79,平均絕對(duì)誤差為82.36,平均絕對(duì)誤差百分比為9.08%。結(jié)論ARIMA(1,0,0)(1,1,0)12模型能較好地?cái)M合成都市成華區(qū)狂犬病暴露數(shù)變動(dòng)趨勢(shì),適用于狂犬病暴露人數(shù)的預(yù)測(cè)。

狂犬病 ARIMA模型 預(yù)測(cè)

狂犬病是由狂犬病毒引起的一種侵犯中樞神經(jīng)系統(tǒng)為主的急性人獸共患傳染病。王梅、周航等人[1]的研究顯示,狂犬病的發(fā)病具有明顯的季節(jié)性,發(fā)病高峰在夏、秋季。開(kāi)展狂犬病暴露監(jiān)測(cè)數(shù)據(jù)的分析應(yīng)用,是合理制定疫苗需求計(jì)劃、犬傷處置費(fèi)用報(bào)銷(xiāo)政策等狂犬病防控措施的有力保障。ARIMA模型是時(shí)間序列分析中較為成熟和應(yīng)用較為廣泛的方法之一,通過(guò)差分對(duì)序列線性趨勢(shì)、周期性等確定性信息進(jìn)行提取,具有充分利用歷史觀測(cè)值,短期預(yù)測(cè)效果較好的優(yōu)勢(shì)。本文利用成都市成華區(qū)2009-2015年狂犬病暴露監(jiān)測(cè)月報(bào)表數(shù)據(jù)構(gòu)建ARIMA模型,預(yù)測(cè)2016年逐月狂犬病暴露數(shù),為該區(qū)2016年人用狂犬病疫苗、免疫球蛋白的需求計(jì)劃和犬傷處置報(bào)銷(xiāo)費(fèi)用預(yù)算的制定提供參考依據(jù)。

資料與方法

1.資料來(lái)源

數(shù)據(jù)來(lái)源于2009-2015年成華區(qū)各家犬傷門(mén)診的《成都市成華區(qū)犬傷門(mén)診月報(bào)表》。

2.ARIMA模型建模過(guò)程

(1)時(shí)間變量的定義與序列平穩(wěn)化 將2009年1月-2015年12月成華區(qū)犬傷月報(bào)表數(shù)據(jù)序列的時(shí)間單位定義為年份、季度、月份型。通過(guò)時(shí)間序列圖觀察序列的平穩(wěn)性,對(duì)不平穩(wěn)的序列進(jìn)行數(shù)據(jù)轉(zhuǎn)化、差分處理,使其成為零均數(shù)的平穩(wěn)時(shí)間序列,達(dá)到以下要求[2]:均數(shù)不隨時(shí)間變化;方差不隨時(shí)間變化;自相關(guān)系數(shù)與所在的時(shí)間點(diǎn)無(wú)關(guān),僅與時(shí)間間隔距離有關(guān)。

(2)模型識(shí)別和定階 自回歸移動(dòng)平均模型ARIMA(p,d,q),可以寫(xiě)成Φ(B)wt=θ(B)ut,其中wt為zt的d階差分;ut為隨機(jī)剩余項(xiàng)(又稱(chēng)為白噪聲),p是自回歸的階數(shù),d指差分的次數(shù),q是移動(dòng)平均的階數(shù)。包含季節(jié)趨勢(shì)的時(shí)間序列可通過(guò)ARIMA(p,d,q)(P,D,Q)s過(guò)程來(lái)擬合,其中P,D,Q,s分別是季節(jié)性自回歸階數(shù)、季節(jié)差分次數(shù)、季節(jié)性移動(dòng)平均階數(shù)、季節(jié)周期。定階即利用自相關(guān)圖(ACF)、偏自相關(guān)圖(PACF)和互相關(guān)圖(CCF)確定p、d、q等參數(shù)的過(guò)程,首先根據(jù)ACF、PACF圖形的截尾或拖尾情況進(jìn)行模型的初步擬合,再根據(jù)擬合的結(jié)果進(jìn)行相應(yīng)調(diào)整、不斷修正,并結(jié)合評(píng)價(jià)指標(biāo)由低階向高階選擇模型的辦法確定模型的p、d和q。

(3)參數(shù)的估計(jì) 依據(jù)BIC確定模型的階數(shù),采用最小二乘法估計(jì)出p個(gè)自回歸參數(shù)φ1、φ2…φp,q個(gè)移動(dòng)平均參數(shù)θ1、θ2…θq,P個(gè)季節(jié)自回歸參數(shù)Ф1、Ф2…ФP,以及Q個(gè)季節(jié)移動(dòng)平均參數(shù)

(4)模型的檢驗(yàn)與優(yōu)化 根據(jù)平穩(wěn)的R方、正態(tài)化的BIC準(zhǔn)則判斷模型的擬合優(yōu)度(BIC值相對(duì)較小的模型較好)。由Ljung-BoxQ檢驗(yàn)結(jié)果對(duì)模型殘差序列進(jìn)行白噪聲檢驗(yàn),判斷所建模型的適合性。在所有通過(guò)檢驗(yàn)的模型中,確定標(biāo)準(zhǔn)化的BIC值較小,模型較為簡(jiǎn)潔的為最優(yōu)模型[4]。

(5)模型的驗(yàn)證和預(yù)測(cè) 以2015年1-12月犬傷暴露數(shù)據(jù)為驗(yàn)證樣本,采用均方根誤差(RMSE)、平均絕對(duì)誤差(MAE)、平均絕對(duì)百分比誤差(MAPE)作為評(píng)價(jià)指標(biāo),通過(guò)比較預(yù)測(cè)值和真實(shí)值來(lái)評(píng)價(jià)模型的預(yù)測(cè)精度。選用較優(yōu)模型對(duì)2016年狂犬病逐月暴露數(shù)進(jìn)行預(yù)測(cè)。

3.統(tǒng)計(jì)軟件處理

利用SPSS 19.0統(tǒng)計(jì)軟件建立ARIMA模型,P<0.05表示差異有統(tǒng)計(jì)學(xué)意義。

結(jié) 果

1.成華區(qū)狂犬病暴露監(jiān)測(cè)數(shù)據(jù)基本情況

成都市成華區(qū)2009-2015年狂犬病暴露逐月監(jiān)測(cè)數(shù)呈現(xiàn)明顯的季節(jié)性波動(dòng),每年從1月開(kāi)始呈逐月上升趨勢(shì),至7、8月份達(dá)到峰值,之后逐月下降。2009-2015年的狂犬病暴露數(shù)最低為8441例,最高為10285例,年平均數(shù)為9023例,年暴露監(jiān)測(cè)數(shù)基本圍繞該平均值波動(dòng)。

2.序列的平穩(wěn)化

觀察原始序列自相關(guān)圖發(fā)現(xiàn)序列有周期性變化規(guī)律,周期為12個(gè)月。原始偏自回歸函數(shù)在k=1后呈現(xiàn)余弦衰減波形。因此,需對(duì)原始序列進(jìn)行一次季節(jié)性差分,觀察差分后的自相關(guān)和偏相關(guān)分析圖,可見(jiàn),其自相關(guān)函數(shù)在k=1后呈現(xiàn)衰減趨勢(shì),僅當(dāng)k=1和k=12時(shí),自回歸系數(shù)明顯突破了可信區(qū)間界值,其偏自回歸函數(shù)在k=1后呈現(xiàn)逐漸衰減至零的趨勢(shì),并落入可信區(qū)間,此時(shí)的時(shí)間序列已基本趨于平穩(wěn),見(jiàn)圖1。

圖1 成都市成華區(qū)狂犬病暴露監(jiān)測(cè)數(shù)一階季節(jié)差分自相關(guān)和偏相關(guān)分布圖

3.模型的識(shí)別和定階

由于成都市成華區(qū)2009年1月-2014年12月狂犬病暴露數(shù)序列存在明顯的季節(jié)性趨勢(shì),季節(jié)性周期為12個(gè)月,故選用ARIMA(p,d,q)(P,D,Q)S過(guò)程來(lái)擬合建模。對(duì)序列進(jìn)行了1次季節(jié)性差分,因此確定d=0,D=1。根據(jù)一階季節(jié)差分的ACF圖,k=1后函數(shù)呈現(xiàn)衰減趨勢(shì),只有k=1、k=2和k=12時(shí),自回歸系數(shù)突破了可信區(qū)間界值,可以選擇q=1或2,Q=1。根據(jù)一階季節(jié)差分的PACF圖,k=1后函數(shù)呈現(xiàn)衰減趨勢(shì),可以選擇p=1,P=1。因此可以選出5組模型,ARIMA(1,0,1)(1,1,1)12、ARIMA(1,0,0)(1,1,1)12、ARIMA(1,0,0)(1,1,0)12、ARIMA(2,0,1)(1,1,1)12、ARIMA(2,0,0)(1,1,0)12。

4.模型參數(shù)估計(jì)

模型ARIMA(1,0,0)(1,1,0)12的標(biāo)準(zhǔn)化BIC值=9.304,在擬合的所有模型中最小;模型擬合效果度量Ljung-Box Q差異無(wú)統(tǒng)計(jì)學(xué)意義(Q=22.076,P=0.141),模型的殘差為白噪聲;且模型中的參數(shù)檢驗(yàn)均有意義,見(jiàn)表1,說(shuō)明所擬合模型是有效的。

表1 成都市成華區(qū)狂犬病暴露監(jiān)測(cè)數(shù)的ARIMA模型參數(shù)估計(jì)

5.預(yù)測(cè)效果分析

應(yīng)用模型ARIMA(1,0,0)(1,1,0)12對(duì)2009年1月-2015年12月的狂犬病暴露數(shù)進(jìn)行回代預(yù)測(cè),結(jié)果顯示擬合值與實(shí)際值基本保持一致的曲線趨勢(shì),且實(shí)際值均在預(yù)測(cè)可信區(qū)間內(nèi),預(yù)測(cè)均方根誤差為28.79,平均絕對(duì)誤差為82.36,平均絕對(duì)誤差百分比為9.08%,見(jiàn)圖2和表2。同時(shí)預(yù)測(cè)2016年狂犬病暴露數(shù)顯示,2016年各月狂犬病暴露數(shù)的趨勢(shì)繼續(xù)跟歷年數(shù)據(jù)趨勢(shì)一致,暴露高峰將出現(xiàn)在5~8月,見(jiàn)圖2和表3。

圖2 成都市成華區(qū)狂犬病暴露監(jiān)測(cè)數(shù)ARIMA(1,0,0)(1,1,0)12預(yù)測(cè)模型擬合圖

表2 ARIMA(1,0,0)(1,1,0)12模型預(yù)測(cè)成華區(qū)2015年狂犬病暴露數(shù)的驗(yàn)證結(jié)果

表3 ARIMA(1,0,0)(1,1,0)12模型對(duì)成華區(qū)2016年狂犬病暴露數(shù)的預(yù)測(cè)結(jié)果

討 論

ARIMA模型是一種精確度較高的短期預(yù)測(cè)方法,通過(guò)季節(jié)性差分和非季節(jié)性差分削弱序列趨勢(shì)性及季節(jié)周期性的干擾,并結(jié)合模型參數(shù)對(duì)時(shí)間序列進(jìn)行擬合和預(yù)測(cè)[5]。在利用ARIMA模型對(duì)時(shí)間序列進(jìn)行預(yù)測(cè)時(shí),為保證模型的預(yù)測(cè)精度至少需要50個(gè)以上的歷史統(tǒng)計(jì)數(shù)據(jù)[6]。本次研究利用過(guò)去6年的狂犬病暴露逐月監(jiān)測(cè)數(shù)據(jù),建立ARIMA(1,0,0)(1,1,0)12模型對(duì)成都市成華區(qū)狂犬病暴露數(shù)據(jù)進(jìn)行預(yù)測(cè)。在實(shí)際應(yīng)用中,需不斷用新的數(shù)據(jù)對(duì)已建模型進(jìn)行修正,提高模型預(yù)測(cè)的精度[7]。吳家兵[8]等人提出如果網(wǎng)絡(luò)模型預(yù)測(cè)對(duì)象的慣性趨勢(shì)發(fā)生了較大的變化(如采取了新的防控措施),則需要收集新的數(shù)據(jù)對(duì)模型進(jìn)行修正或重新擬合。

目前,ARIMA模型已廣泛應(yīng)用于傳染病發(fā)病的預(yù)測(cè)[9-12]。在模型的擬合過(guò)程中,首先對(duì)原始時(shí)間序列進(jìn)行觀察,如果未達(dá)平穩(wěn)化要求,則進(jìn)行差分或(和)季節(jié)差分,使其達(dá)到平穩(wěn)化的要求,確定D或(和)d。隨后通過(guò)對(duì)ACF圖和PACF圖的觀察識(shí)別,對(duì)自回歸模型和移動(dòng)平均模型的p、q進(jìn)行定階,產(chǎn)生幾個(gè)試用模型。依據(jù)BIC值診斷模型的擬合優(yōu)度,并根據(jù)簡(jiǎn)潔、殘差不相關(guān)的原則篩選出最優(yōu)模型。本研究最終確定的最優(yōu)模型為ARIMA(1,0,0)(1,1,0)12。該模型對(duì)成都市成華區(qū)2009-2015年狂犬病逐月暴露監(jiān)測(cè)數(shù)據(jù)實(shí)際值進(jìn)行了較好的擬合,回代預(yù)測(cè)2015年1月-12月狂犬病暴露數(shù)預(yù)測(cè)值與實(shí)際值的平均絕對(duì)誤差百分比為9.08%,顯示預(yù)測(cè)數(shù)據(jù)與實(shí)際數(shù)據(jù)吻合程度較高,提示利用ARIMA(1,0,0)(1,1,0)12模型能對(duì)狂犬病暴露數(shù)進(jìn)行較好的預(yù)測(cè)。

全國(guó)狂犬病年暴露人數(shù)逾4000萬(wàn)[13]。狂犬病暴露的監(jiān)測(cè)工作對(duì)狂犬病防控效果的評(píng)估、防控策略的調(diào)整及疫情趨勢(shì)的預(yù)測(cè)分析均具有重要意義[14]。本研究結(jié)果顯示:成都市成華區(qū)2009年-2015年的狂犬病暴露數(shù)圍繞9023例的年平均值上下波動(dòng),說(shuō)明成都市成華區(qū)近年積極開(kāi)展狂犬病暴露監(jiān)測(cè)及暴露后的規(guī)范化處置、犬傷處置費(fèi)用限額報(bào)銷(xiāo)及健康教育等綜合防制措施的成效顯著。ARIMA(1,0,0)(1,1,0)12模型的擬合預(yù)測(cè)曲線顯示:每年狂犬病逐月暴露數(shù)呈季節(jié)性的單峰分布,7、8月份為高峰值月,這與王梅[1]等人研究結(jié)果一致;對(duì)2016年的預(yù)測(cè)結(jié)果顯示:2016年逐月狂犬病暴露數(shù)與歷年數(shù)據(jù)趨勢(shì)一致,暴露高峰將出現(xiàn)在5-8月。

結(jié)合本研究的結(jié)果,在下一步工作中應(yīng)繼續(xù)加強(qiáng)狂犬病暴露數(shù)據(jù)監(jiān)測(cè),不斷收集和使用新的數(shù)據(jù)修正預(yù)測(cè)模型,以提高預(yù)測(cè)精度,更好地提供參考依據(jù);根據(jù)預(yù)測(cè)值和趨勢(shì),科學(xué)制定人用狂犬病疫苗的需求計(jì)劃,做好經(jīng)費(fèi)預(yù)算,促進(jìn)犬傷處置費(fèi)用報(bào)銷(xiāo)政策的順利實(shí)施;在暴露高峰期加大健康教育工作力度。

[1]王梅,周航,殷文武,等.中國(guó)2005-2011年人狂犬病不同地區(qū)季節(jié)分布特征研究.中華流行病學(xué)雜志,2012,33(11):1151-1154.

[2]張文彤,董偉.SPSS統(tǒng)計(jì)分析高級(jí)教程.第2版.高等教育出版社,2013:395-398.

[3]陳斌,周伴群,焦亮,等.ARIMA模型在狂犬病暴露監(jiān)測(cè)中的應(yīng)用.中國(guó)預(yù)防醫(yī)學(xué)雜志,2011,12(5):427-430.

[4]潘浩,鄭楊,吳寰宇,等.ARIMA模型預(yù)測(cè)上海市手足口病發(fā)病趨勢(shì).預(yù)防醫(yī)學(xué)情報(bào)雜志,2011,27(6):408-411.

[5]丁磊,丁淑軍,張萌,等.應(yīng)用時(shí)間序列分析研究秋冬型恙蟲(chóng)病時(shí)間分布特征及趨勢(shì).中華流行病學(xué)雜志,2012,33(7):698-701.

[6]劉重程,李宏通,唐雅清,等.ARIMA模型在細(xì)菌性痢疾預(yù)測(cè)中的應(yīng)用.中國(guó)預(yù)防醫(yī)學(xué)雜志,2011,12(10):842-844.

[7]張?jiān)剑鮿匐y,劉媛,等.應(yīng)用ARIMA模型對(duì)呼吸系統(tǒng)疾病月住院量及住院費(fèi)用的預(yù)測(cè).中國(guó)衛(wèi)生統(tǒng)計(jì),2015,32(2):197-200.

[8]吳家兵,葉臨湘,尤爾科.ARIMA模型在傳染病發(fā)病率預(yù)測(cè)中的應(yīng)用.?dāng)?shù)理醫(yī)藥雜志,2007,20(1):90-92.

[9]李驪,錢(qián)俊,楊軍,等.三種模型對(duì)廣東省副傷寒逐月發(fā)病數(shù)預(yù)測(cè)的比較.中國(guó)衛(wèi)生統(tǒng)計(jì),2014,31(2):197-201.

[10]焦亮,阮峰,黃利群,等.基于ARIMA的流感癥狀預(yù)測(cè)模型.實(shí)用預(yù)防醫(yī)學(xué),2010,17(8):1482-1486.

[11]吳偉,郭軍巧,安淑一,等.應(yīng)用ARIMA-GRNN模型對(duì)腎綜合征出血熱發(fā)病率時(shí)間序列數(shù)據(jù)的預(yù)測(cè)研究.中國(guó)衛(wèi)生統(tǒng)計(jì),2015,32(2):211-213.

[12]陳偉,陳正利,李少芳,等.ARIMA模型在河南省梅毒月發(fā)病率預(yù)測(cè)中的應(yīng)用.中國(guó)衛(wèi)生統(tǒng)計(jì),2013,30(4):604-606.

[13]Cui PY,Hang Z,HuiW.Analysis on Factors Related to Rabies Epidemic in China from 2007-2011.Virologica Sinica,2012,27(2):132-143.

[14]周興余,劉學(xué)成,張佳珂.2010年四川省狂犬病監(jiān)測(cè).預(yù)防醫(yī)學(xué)情報(bào)雜志,2012,28(1):34-37.

(責(zé)任編輯:劉 壯)

ARIMA Model in Prediction of the Rabies Surveillance in Chenghua District of Chengdu City

Yang Jing,Zhang Qiang
(Department of Epidemiology and Health Statistics,School of West China Public Health,Sichuan University(610041),Chengdu)

ObjectiveTo understand the basic characteristics of rabies surveillance data in Chenghua district of Chengdu.To establish ARIMA model based on rabies surveillance data collected from recent years,and to forecast the rabies exposed data from January to December in 2016 in Chenghua district,with the purpose of providing evidence for the standardization construction of dog injury clinic and demand planning of rabies vaccine and rabies immunoglobulin.MethodsSPSS 19.0 was used to establish ARIMA model based on the monthly rabies surveillance data from 2009 to 2014,and case numbers of rabies exposed from January to December in 2015 were used as sample to examine the model accuracy.The optimal ARIMA model was used to predict the monthly rabies exposed numbers in 2016.ResultsThe optimal ARIMA model was ARIMA(1,0,0)(1,1,0)12,which could fit very well with the rabies exposure series in the past period of time.Case numbers of rabies exposed from January to December in 2015 were used as sample to exam ine the model accuracy,the results showed all actual values fell in the 95% confidence intervals of expected values,the mean square error was 28.79,the mean absolute error was 82.36,the mean absolute percentage error was9.08%.ConclusionARIMA(1,0,0)(1,1,0)12could simulate the trend of rabies exposure in the Chenghua district of Chengdu,and can be applied for forecasting the case number of rabies exposed.

Rabies;Surveillance;ARIMA model;Prediction

1.四川大學(xué)華西公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系(610041)

2.成都市成華區(qū)疾病預(yù)防控制中心

△通信作者:張強(qiáng),E-mail:qiangzhang@scu.edu.cn

猜你喜歡
趨勢(shì)模型
一半模型
趨勢(shì)
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
初秋唇妝趨勢(shì)
Coco薇(2017年9期)2017-09-07 21:23:49
3D打印中的模型分割與打包
SPINEXPO?2017春夏流行趨勢(shì)
“去編”大趨勢(shì)
FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
趨勢(shì)
主站蜘蛛池模板: 日本三级黄在线观看| 久久一级电影| 国产黄色免费看| 国产丰满大乳无码免费播放| 天堂网亚洲系列亚洲系列| 精品综合久久久久久97超人| av色爱 天堂网| 超级碰免费视频91| 青草视频在线观看国产| 亚洲侵犯无码网址在线观看| 国产精品入口麻豆| 黄色免费在线网址| 在线播放国产99re| 亚洲69视频| 伊人欧美在线| 无码精品福利一区二区三区| 国产在线小视频| 在线欧美日韩国产| 日本高清在线看免费观看| 伊人久久大香线蕉综合影视| 伊人狠狠丁香婷婷综合色| 欧美久久网| 尤物精品视频一区二区三区| 国产综合亚洲欧洲区精品无码| 国产精品亚洲专区一区| 国产日韩欧美在线播放| 欧美啪啪一区| 亚洲无线观看| 中文字幕无码中文字幕有码在线| 国产成人高清亚洲一区久久| 人禽伦免费交视频网页播放| 国产一区二区影院| 欧美一级在线看| 亚洲性影院| 99精品热视频这里只有精品7| 久久a级片| 性欧美久久| 毛片一级在线| 国产网站黄| 国产精品99一区不卡| 成人欧美日韩| 九九热免费在线视频| 久草中文网| 欧美成人aⅴ| 国产流白浆视频| 91美女视频在线| 亚洲成A人V欧美综合| 91精品综合| 欧美a√在线| 久久久精品久久久久三级| 久久久精品国产SM调教网站| 午夜限制老子影院888| 亚洲成人在线免费| 免费在线色| 国产女人在线| 亚洲性一区| 国产99视频精品免费视频7| 免费看黄片一区二区三区| 亚洲无码精品在线播放| 婷婷中文在线| 精品一区二区无码av| 欧美五月婷婷| 免费在线看黄网址| 国产不卡网| 色综合a怡红院怡红院首页| 丰满人妻一区二区三区视频| 日韩区欧美区| 亚洲成人在线网| 精品国产三级在线观看| 国产区91| 欧美特黄一级大黄录像| 538精品在线观看| 久久综合九色综合97网| 国模沟沟一区二区三区| 亚洲中文字幕久久无码精品A| 久久精品中文字幕免费| 国产爽歪歪免费视频在线观看| 特级精品毛片免费观看| 色精品视频| 91探花国产综合在线精品| 狼友视频一区二区三区| 国内99精品激情视频精品|