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

四類預測人口方法的對比及Logistic人口生長模型的改進

2022-05-31 11:40:04王沛林
保定學院學報 2022年3期
關鍵詞:生長效果模型

王沛林

(北京理工大學 數學與統計學院,北京 100081)

在當今科技經濟飛速發展的社會中,人口問題一直是世界各個國家和各個領域所關注的焦點問題,它是一個社會發展的最基礎也是最關鍵的問題.為了維持社會生活的基本平衡與穩定,把人口數量保持在一個比較理想化的狀態,就要對國家和地區的人口總數進行合理預測,以便制定未來人口政策,使社會有序發展.在早期的研究中,進行人口定量預測的基礎方法就是依據國家統計局的人口結構和規模的資料和數據,以當前人口為出發點,并且對未來人口的變化趨勢作出合乎常理的假設,運用科學的數學模型進行擬合演算,進而預測出未來人口總數、性別比例、結構組成等一些重要人口要素.

社會的進步需要人口穩定可持續發展,基于中國人口數量現狀和上述已知方法,本文通過比較馬爾薩斯人口模型、多項式擬合模型、ARIMA模型和Logistic人口生長模型的預測效果,著眼于未來人口總數的估計分析,關注人口發展的總體趨勢,并將現有Logistic人口生長模型從微分方程和模型參數擬合優化角度進行改進,得到了更加準確的預測方程,以便更好地預測今后一段時間中國的總人口數,具有一定的預見性.

1 中國人口發展預測模型的建立

1.1 數據的選取

擬合數據的選取十分重要.經過資料調查,我們知道從20世紀70年代,在全國范圍內開始大力實施計劃生育政策[1],所以在預測全國總人口時,我們選取1970年為基礎年份,1970年的總人口數量為最初人數.進而考慮到數據選取的原則,數據選擇得越多,即樣本越多,能準確預測擬合人口數量的可能性越大,結果就越精確,誤差就越小,所以選擇的年份不能少于25年,否則將很難發現人口數量變動的趨勢.再有,我國自2016年1月1日起,正式開始實行“全面二孩”政策[2],這項政策的實施在短期內不會對我國總人口數的變化趨勢產生較大的影響.因此,在本文中先不考慮近幾年二胎方針的影響.最后通過查閱中國國家統計局網站,獲得1970—2014年每年中國人口總數,以這45年的數據作為擬合數據,保證了一定的準確性.

1.2 模型假設

1)現有的基本國策“計劃生育”保持不變.

2)不考慮戰爭等不可預見因素.

3)由于環境等因素,人口不能無限制地增長.

1.3 馬爾薩斯人口模型及其應用

1.3.1 模型的理論背景

人口學家馬爾薩斯在生物增長定律中提到:人口的變化率與人口總數成正比[3],且這個比率不隨時間變化,為常量.由此,我們設定人口總數為p(t),人口變化率為l,若已經知道在t0時刻的人口數量值p0,則就可以得到馬爾薩斯人口模型:

這就是非常著名的分析人口數量問題的基礎微分方程模型,通過求解方程,我們可以得到它的另外一種表達形式:

其中,t0表示最開始的時間(即初始年份),p0表示t0時刻的人口總數,t表示一般意義上的時間,p(t)表示t時刻的人口總數.

1.3.2 模型的建立與實現

將1970—2014年每一年的人口總數作為基礎數據,代入模型中.為了進行數據擬合,首先要考慮模型的求解實現問題.

參數估計:l、p0可以用已知數據,運用線性最小二乘法進行估計.(2)式兩邊取對數,得到:

作對數變換后,再以1970—2014年的人口數據擬合(3)式,經計算得到:l=0.011 2,p0=88 521.43.

最后,得到馬爾薩斯人口模型:

其中t0=1970.

得到模型的表達式以后,我們進而對2015—2019年的人口進行預測.

1.3.3 模型的誤差檢驗分析

評價一個模型的好壞以及擬合效果,首先可以通過實際值和擬合值作出的圖像來直觀地對比分析.

如圖1所示,在馬爾薩斯人口模型預測圖中,這45年來實際人口和預測人口的總體增長趨勢大體是一致的,但是到后期預測人口的增長率要明顯大于實際人口的增長率.

圖1 馬爾薩斯人口模型預測1970—2014年總人口

通過計算誤差來進一步分析模型的優劣.在統計學中,我們慣用yi來代表第i組數據的真實值,來代表第i組數據的預測值,用n來表示樣本數據的組數,所以有誤差指標為:相對誤差(RE),殘差平方和(resnorm):,標準誤差(RMSE):

相對誤差越小,模型可信度越高;殘差平方和越小,模型擬合效果越好;標準誤差越小,模型精密度越高;Rnew為非線性回歸方程的擬合優度指標,其值越靠近1,模型擬合程度越高[4].

知悉這幾個誤差指標后,用Matlab計算出由馬爾薩斯人口模型擬合的這45年數據的相對誤差最大值為0.066 6,殘差平方和為4.941 7,標準誤差為0.331 4,Rnew為0.971 3.

在物理學中,一般相對誤差不超過5%算是比較精確[5],但此模型最大相對誤差0.066 6要高于0.05,并且Rnew的值離1還有一定距離,由此看來馬爾薩斯人口模型的精確度不是很高.

1.3.4 模型的預測效果及總體評價

將前文利用馬爾薩斯人口模型預測出的2015—2019年中國總人口數與國家統計局網站上的真實數據進行對比,以此來判斷預測效果.

由表1可以計算出,這5年的相對誤差最大值為0.094 6,遠大于0.05,誤差較大,預測效果并不是很理想.

表1 中國2015—2019年實際總人口數與預測總人口數 萬人

進一步繪制2015—2030年的人口預測圖,從圖2中可以看出預測人口的年增長率過高,人口數量呈幾何狀增長,不太符合實際.

圖2 馬爾薩斯人口模型預測2015—2030年總人口

從馬爾薩斯人口模型的擬合效果和預測結果來看,計算出來的誤差較大,模型并不是很精確,其預測的人口數量呈指數增長,增長率過高,故該模型的現代利用價值有限.

1.4 多項式擬合模型及其應用

1.4.1 模型的理論背景

從圖1中我們可以觀察出,人口數量的變化過程線上各點斜率都不太相同,這就表明各個階段的人口增長速率是不一樣的,這時用單純的一條直線來描繪人口的發展顯然是不合適的,所以要用更貼近實際的曲線來擬合人口數量的變動,這就產生了多項式擬合模型.

1.4.2 模型的建立與實現

為了保證多項式擬合模型的合理性和準確性,并且兼顧誤差最小的原則,還要使模型具有可解釋性,故本文選取了三次多項式對1970—2014年的總人口數據進行擬合.

將時間t作為解釋變量,中國總人口數量p(t)作為預測變量,建立多項式確定它們之間的關系.得出如下結果:

1.4.3 模型的誤差檢驗分析

通過圖3擬合圖像可以看出,在1970—2014年中,擬合的曲線很貼近實際的數據,三次多項式擬合模型的短期預測效果很好,沒有太大的誤差.

圖3 三次多項式擬合模型預測1970—2014年總人口

三次多項式的相對誤差最大值為0.0133,殘差平方和為0.122 3,標準誤差為0.052 1,Rnew為0.995 5.由誤差可以看出,多項式擬合模型的精確度較馬爾薩斯人口模型有很大程度的提高.

1.4.4 模型的預測效果及總體評價

利用三次多項式擬合模型的表達式,預測出2015—2019年中國人口總數分別為136 090萬人、136 194萬人、136 218萬人、136 161萬人、136 021萬人,將這5年來的預測人口同實際人口相比較,我們得出最大的相對誤差值為0.028 5,沒有超過0.05,且要遠低于馬爾薩斯人口模型的誤差值,預測效果比較好.

從三次多項式擬合模型的曲線貼合效果和部分預測結果來看,各種誤差都比較小,模型比較合理,但是2015—2019年預測出來的人口數量總體呈現先增長后下降的趨勢,人口數量下降在我國現階段幾乎是不可能的,所以筆者認為多項式擬合模型只適合短時段1~3年內的人口預測,長期預測則會出現大的偏差,效果不會很好.

1.5 ARIMA模型及其應用

1.5.1 模型的理論背景

由前面的多項式擬合模型可以看出,人口數量是隨著時間不斷變化的.但是人口與時間的關系又不能簡單地用多項式來準確的描述,為了尋找更合理的人口預測模型,我們想到了經濟學中的ARIMA模型.

預測一種事物的變化時,用它的過去推斷它的未來,即用時間序列的往日數據表示事物隨時間變化的種種規律,并將兩者之間的模式應用到未來,進而對未來的數據作出預測.ARIMA模型就是一種常用的時間序列模型.

ARIMA(p,d,q)模型含有 3個參數,這里面 p表示自回歸(AR)階數,q表示移動平均(MA)階數,d 表示模型的差分階數[6].該模型有 3 種基本類別[7]:AR(p)模型、MA(q)模型、ARIMA(p,d,q)模型.

1.5.2 模型的建立與實現

利用1970—2014年的中國總人口數據,建立ARIMA(p,d,q)模型.模型的建立步驟如下:

1)載入數據,進行一次單位根檢驗,判定數據序列的平穩性.本文得到p值為0.99,大于檢驗的臨界值0.05,所以該序列不平穩.

2)對不平穩序列進行差分處理,差分的次數為d,使其成為平穩序列.經過了三次差分操作后,一次單位根檢驗的p值小于了0.05,序列化為平穩序列,進而確定d=3.

3)對已經平穩的序列進行白噪聲檢驗,判斷其是否為純粹的隨機序列,由于對白噪聲進行下一步的處理和預測沒有任何價值,故使用Ljung-Box方法[8]進行檢測.p值為0.008 808,小于標準值0.05[8],所以三次差分后的平穩序列不是隨機序列,可以展開下一步工作.

4)確定模型的階數.畫出序列的自相關圖(ACF)和偏自相關圖(PACF),見圖4.據觀察可以發現ACF 1階截尾,PACF具有拖尾性.綜上,初定模型為ARIMA(0,3,1).接下來用R語言軟件中的自動定階功能(auto.arima)來對上面的定階結果進行準確性檢驗,結果也顯示為ARIMA(0,3,1)模型.

圖4 自相關圖與偏自相關圖

5)確定模型的系數.得出該模型為xt=εt+0.769εt-1.

6)模型的顯著性檢驗.使用Ljung-Box方法[8]檢驗殘差:p 值為 0.287 8,大于 0.05,說明殘差是隨機序列,白噪聲檢驗通過.參數的顯著性檢驗:用R軟件測出的系數-0.769 0除

以它的標準誤差0.122 7,商的絕對值為1.96,大于T檢驗統計量在5%水平的臨界值,斷定系數顯著不等于零[9],檢驗通過.

7)用 ARIMA(0,3,1)模型預測 2015—2019年的中國人口總數.

1.5.3 模型的誤差檢驗分析

從圖5可以看出,ARIMA模型作出的預測數據與實際數據比較吻合,相差不大,效果良好.

圖5 ARIMA模型預測1970—2014年總人口

相對誤差最大值為0.003 3,殘差平方和為0.006 3,標準誤差為0.011 8,Rnew為0.999 0.從誤差值來看,模型的擬合水平非常高,較前面2種模型而言,有著更高的可信度.

1.5.4 模型的預測效果及總體評價

利用ARIMA模型預測2015—2019年的中國人口總數,將預測值與實際值作比較,得出相對誤差的最大值為0.002 8,遠小于前2種模型的誤差值,預測效果有了質的飛躍.

進而利用此模型預測出了2015—2030年中國的總人口數,并繪成圖6.看出未來10年我國人口將呈緩慢增長態勢,到2030年不會破15億大關.

圖6 ARIMA模型預測2015—2030年總人口

根據綜上各種論證來看,ARIMA模型的預測效果遠優于前2種模型,給出了較為合理的人口隨時間變化的趨勢,精確度比較高,有較大的參考價值.

1.6 Logistic人口生長模型及其應用

1.6.1 模型的理論背景

在前面討論過的馬爾薩斯人口模型中有致命的缺點——人口變化率l為常量.它只假設出了人類發展最理想的情況,并得出人口數量將呈指數上漲的極端結論.認識到馬爾薩斯人口模型的不足之后,有學者在其基礎上進行了改進,形成了Logistic人口生長模型.模型的中心觀點為:人口數量并不能無窮盡地上漲,它會受到各種因素的約束,且隨著人口的不斷增加,這種約束力會不斷變強,最終人數會到達一個極值.

1.6.2 模型的建立與實現

從上述的理論基礎中可知,如果想獲取方程(10)的解,就要估計出pm和l這2個參數的值,求解方法有2種.

1.6.2.1 三點等間距法

利用昔日的人口數據,粗略計算出l和pm的值.再采用t0、t1、t2這3個年份的人口總數量p0、p1、p2,且滿足條件t2-t1=t1-t0=a,代入方程式(10)進行計算,可得出:

取1970年作為開端年份t0,取1992年為中間年份t1,取2014年為最末年份t2,此時a=22,將相應數值代入公式(11),獲取l和pm的具體值;將計算出的2個參數值代入方程(10)中,得到模型,模擬圖見圖7.

圖7 Logistic模型三點等間距法預測1970—2014年總人口

1.6.2.2 lsqcurvefit最小二乘法

在Matlab軟件中有一個lsqcurvefit命令,依據最小二乘原理即使模擬數據和實際數據之間的殘差平方和最小,也可擬合出最優化的非線性曲線.

選取 t0=1970,p0=82 992代入式(10),然后將式(10)作為未知的需要擬合的非線性曲線,l和pm為要求解的參數.給出兩參數的初值解,1970—2014年的總人口為實際數據,用lsqcurvefit函數進行擬合,pm為153 147 萬人,l為 0.045 4,有模型,最后畫出對比圖,見圖8.

圖8 Logistic模型最小二乘法預測1970—2014年總人口

1.6.3 模型的誤差檢驗分析

從圖7、圖8中可以剖析出,三點等間距法和最小二乘法的模擬結果都不錯.

三點等間距法:相對誤差最大值為0.018 2,殘差平方和為0.305 0,標準誤差為0.082 3,Rnew為0.992 9.最小二乘法:相對誤差最大值為0.013 5,殘差平方和為0.247 7,標準誤差為 0.074 2,Rnew為 0.993 6.很明顯,最小二乘法還是要略優于三點等間距法,這可能是因為lsqcurvefit函數能夠找到全局最優解,而三點等間距法只是拿出3年的數據來求解方程(10),有局限性.

1.6.4 模型的預測效果及總體評價

綜上,我們選用最優的lsqcurvefit最小二乘法預測出了2015—2019年的總人口數,與實際人口比,最大相對誤差值為0.003 9.進而本研究用Logistic人口生長模型預測出了2015—2030年的總人口,見圖9.

圖9 Logistic人口生長模型預測2015—2030年總人口

從圖9可以發現,在未來10年我國人口增長有逐漸放緩的趨勢,總體較平穩.

總體來說,Logistic人口生長模型比馬爾薩斯人口模型考慮的影響人口變動的因素更多更充分,預測結果也更加確切、合理.但是,從誤差的角度來看,該模型還不如ARIMA模型精確,有改進優化的空間.

2 4種預測方法的比較分析

通過表2我們清楚地認識到,無論是從誤差最小化角度,還是從模型的擬合效果和預測結果方面,Logistic人口生長模型和ARIMA模型都要優于其他兩種模型,并且ARIMA模型的精確性還要比Logistic人口生長模型突出,但是ARIMA模型的局限性在于它只適合短時間內的預測[11],而Logistic人口生長模型作為中長期人口預測的一種典型方法,它的精確度還有待提高.為了更好地進行人口中長期預測,下文將探討Logistic人口生長模型的改良方法.

表2 4種預測模型的比較分析

3 Logistic人口生長模型的改進

通過分析看出,Logistic人口生長模型在預測人口時,其自身的合理性和精確度還有進一步提升的空間,為此我們從2種不同角度探索對此模型修正的方法.

3.1 方程角度改進及其應用

3.1.1 改進的方法依據

從我國的當前人口形勢出發,由于計劃生育政策已實行了多年,我國過去人口猛烈增長的勢頭已經得到了很好的控制,久而久之又出現了生育率低、人口老齡化的現象.如今我國的人口增長率已經非常低,如方程(7)的假設——人口變化率l和人口總數p(t)是最簡單的線性關系,已經不再適用于當下.現在我們從l與p(t)的關系入手,對l(p)函數式進行改進.

作出假設:人口上漲已經快到了負荷的程度,上漲率慢慢變小最終將趨向零.

根據假設,l(p)理應表達為e-p的函數,設:

其中l(p0)是初始年份為t0、總人數為p0時的人口增長率,b、k為待求系數.

于是,修正以后的模型為:

3.1.2 模型的建立與實現

依舊使用1970—2014年的中國人口總數作為研究數據,這樣上下文形成對照,方便比較擬合效果.先來求解方程(12)中的參數值b、k.

首先,對方程(12)的兩邊作對數變換得到:

最后,得到從方程角度進行修改后的模型:

其中起始年份t0為1970年.求出了方程(15)的數值解,并繪制擬合圖像,見圖10.

圖10 方程角度改進的Logistic人口生長模型預測1970—2014年總人口

3.1.3 模型的誤差檢驗分析

根據圖10可看出,從方程角度改進后的模型繪制的擬合曲線與現實中的數據比較貼近,在開始的5~6年和末尾的5~6年預測效果尤其好,在中間年份或多或少出現了偏差.依據程序運行結果,相對誤差最大值是0.026 3,殘差平方和是0.796 3,標準誤差是0.134 5,Rnew是0.988 4.從誤差值的大小來看,此方法修正后的Logistic人口生長模型并沒有較大改善.

3.1.4 模型的預測效果及總體評價

與前文一樣,我們用改良的模型預測出了2015—2019年的中國人口總數,與實際人口比較,這5年的最大相對誤差值為0.004 0,這與沒有修正之前的最大相對誤差值0.003 9差距非常小,幾乎沒有改善.之后,用改進后的模型預測了2015—2030年的中國總人口,見圖11.

圖11 方程角度改進的Logistic人口生長模型預測2015—2030年總人口

將圖11與圖9對比發現,改進前預測2015—2030年的中國總人口大約是從13.8億上漲到14.5億,改進后預測這16年的人口數大約是從13.8億上漲到14.7億多,首尾年都相差不大,說明我國未來人口的變化有很大可能按照這種趨勢發展,并在2030年達到14.6億左右的人口數.

無論是從誤差指標的角度來看,還是從預測效果來看,修正了主方程的Logistic人口生長模型并沒有很大程度上提高模型的效率和精確性,與沒修改之前的預測大致相同,借鑒意義有限.

3.2 參數角度的模擬仿真改進及其應用

3.2.1 改進的方法依據

本文又從另一個角度——使參數值更加精確化,來對模型的修正問題進行探究.

在Logistic人口生長模型中,存在關系函數l(p)=l-sp,l>0,s>0,此式中含有2個待求參數l和s.我們知道,這2個參數主要是通過專家預估[10]或者統計得來,而沒有從函數式本身入手去求其值.所以,修正模型中l(p)的決定方式,根據1970—2014年逐年的總人口增長率與每年人口總數的關系,運用最小二乘線性擬合的方法來明確參數l、s的切實值,即通過真實數據的模擬仿真來確定參數值,而不是無端的估計統計,這樣預測水平可能會有所提高.

3.2.2 模型的建立與實現

應用國家統計局的數據計算出1970—2014年中國總人口的逐年變化率,并制成散點圖.

由圖12可以觀察出,總人口的增長率自1987年開始,幾乎呈直線型下降,因此可以采用直線擬合的方式來得出函數l(p)=l-sp,l>0,s>0的具體表達式.

圖12 1970—2014年中國總人口的增長率變化

于是,1987—2014年逐年的總人口增長率作為l(p),1987—2014年每年的總人口數為p,運用線性最小二乘法來擬合它們之間的關系式l(p)=l-sp,最終得到具體表達式為:

將1987年當作初始年份t0,1987年所對應的總人口數為p0,2個數值與式(16)一起代入方程(6),得到方程:

而后再求解微分方程(17)得到了預測總人口數p(t)的解析式:

實際總人口與預測總人口相比較,見圖13.

圖13 參數角度改進的Logistic人口生長模型預測1987—2014年總人口

3.2.3 模型的誤差檢驗分析

查看擬合圖13,我們發現改進后的模型作出的預測數據與真實數據很接近,總體的增長趨勢也相當一致.

通過計算,得到了該模型下相對誤差最大值是0.006 9,殘差平方和是0.086 0,標準誤差是 0.055 4,Rnew是 0.995 6.改進后模型計算出的誤差指標值要比原模型計算出的誤差值小得多,擬合優度Rnew也更加接近1,結果良好.

3.2.4 模型的預測效果及總體評價

同樣,我們用改良后的模型預測出了2015—2019年的中國總人口數,分別為137 480萬人、137 901萬人、138 298萬人、138 671萬人、139 023萬人,預測值與實際2015—2019年總人口數的最大相對誤差值是0.007 0,與之前未改進時的預測效果相近.接著我們又預測了2015—2030年的中國總人口數,見圖14.

圖14 參數角度改進的Logistic人口生長模型預測2015—2030年總人口

可以看出,按此種模型進行推測,我國的人口增長率減小的趨勢將加快,人口增長將放慢.

利用真實數據對未知函數進行模擬仿真,進而求解出更加有實際意義的參數值,這種改進方法經過實踐證明具有不錯的效果,使模型的擬合優度有了較大提高,從誤差值來看也比原始模型更具有說服力,精確性大幅度提高.

3.3 模型的改進效果剖析

無論是從模型存在的偏差還是從模型的擬合優度指標來看,從參數角度進行的模擬仿真改進遠遠要比從方程入手對模型的改進效果理想得多.改進后的模型,參數值被賦予了更加具體的內涵而且也精確了許多,同時,與原有的Logistic人口生長模型相比,從參數角度改進后的模型精確度提高了,預測水平也升高了.

本文通過比較馬爾薩斯人口模型、多項式擬合模型、ARIMA模型、Logistic人口生長模型這4種人口預測方法的優劣,認為ARIMA模型和Logistic人口生長模型的預測效果比較理想,并用比較合理的Logistic人口生長模型對未來的人口進行了預測,進一步針對Logistic人口生長模型存在的不足提出了2種改進辦法,使模型的精確度進一步提高.最后預測,到2030年,我國的總人口將在14.2億~14.6億波動.

猜你喜歡
生長效果模型
一半模型
按摩效果確有理論依據
碗蓮生長記
小讀者(2021年2期)2021-03-29 05:03:48
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
生長在哪里的啟示
華人時刊(2019年13期)2019-11-17 14:59:54
迅速制造慢門虛化效果
生長
文苑(2018年22期)2018-11-19 02:54:14
抓住“瞬間性”效果
中華詩詞(2018年11期)2018-03-26 06:41:34
3D打印中的模型分割與打包
主站蜘蛛池模板: 色综合久久无码网| 国产福利一区在线| 亚洲高清无在码在线无弹窗| 青青草欧美| 波多野一区| 国产高清无码第一十页在线观看| 5555国产在线观看| 国产成人a在线观看视频| 国产精品夜夜嗨视频免费视频 | 国产香蕉97碰碰视频VA碰碰看| 国产乱子伦视频在线播放| 久久熟女AV| 一级高清毛片免费a级高清毛片| 国产国拍精品视频免费看| 色婷婷亚洲综合五月| 中文字幕无码电影| 免费毛片视频| 亚洲成人动漫在线| 亚洲欧美在线精品一区二区| 992Tv视频国产精品| 日韩欧美国产另类| 午夜视频免费一区二区在线看| 囯产av无码片毛片一级| 国产亚洲成AⅤ人片在线观看| 国产精品粉嫩| 亚洲精品在线影院| 在线看片免费人成视久网下载| 日韩国产欧美精品在线| 久久99热66这里只有精品一| 中国精品自拍| 国产区精品高清在线观看| 一区二区自拍| 国产精品手机视频| 日本道综合一本久久久88| 欧美国产精品不卡在线观看| 国产成人高清精品免费5388| 国产精品hd在线播放| 欧美国产在线看| 久久国语对白| 伊人大杳蕉中文无码| 无码内射在线| 国产又爽又黄无遮挡免费观看| 亚洲最新地址| 国产无码在线调教| 久久人妻xunleige无码| 99久久免费精品特色大片| 国产精品天干天干在线观看| 欧美在线一级片| 国产麻豆永久视频| 欧美成人看片一区二区三区| 精品福利网| 九九热在线视频| 亚洲综合在线网| 国产精品3p视频| 97国产在线观看| 久久夜色撩人精品国产| 一级毛片不卡片免费观看| 一区二区三区毛片无码| 国产精品3p视频| 在线亚洲精品福利网址导航| 免费可以看的无遮挡av无码| 欧美国产日本高清不卡| 亚洲三级a| 污视频日本| 91精品aⅴ无码中文字字幕蜜桃| 久久毛片基地| 国产精品七七在线播放| 亚洲精品综合一二三区在线| 在线观看91香蕉国产免费| www.狠狠| 97在线观看视频免费| 日韩毛片视频| 国产乱人伦AV在线A| 真实国产精品vr专区| 亚洲日韩久久综合中文字幕| 婷婷久久综合九色综合88| 国产杨幂丝袜av在线播放| 国产原创第一页在线观看| 欧美黄网在线| 亚洲永久免费网站| 国产乱子伦精品视频| 无遮挡国产高潮视频免费观看|