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

小樣本條件下可變作用代表值的貝葉斯推斷方法

2015-01-07 07:59:38姚繼濤王旭東
關(guān)鍵詞:方法

姚繼濤, 王旭東

(西安建筑科技大學(xué)土木工程學(xué)院,陜西西安710055)

小樣本條件下可變作用代表值的貝葉斯推斷方法

姚繼濤, 王旭東

(西安建筑科技大學(xué)土木工程學(xué)院,陜西西安710055)

為克服工程實(shí)際中測(cè)試數(shù)據(jù)不充足的條件下,采用經(jīng)典統(tǒng)計(jì)學(xué)方法推斷可變作用的代表值時(shí),由統(tǒng)計(jì)不確定性導(dǎo)致推斷結(jié)果偏于冒進(jìn)的缺點(diǎn),基于極小值I型分布分位值的線性回歸推斷方法,提出了小樣本條件下可變作用代表值的線性回歸推斷方法,并以此作為檢驗(yàn)其他推斷方法精度的基準(zhǔn)方法;根據(jù)貝葉斯理論,利用分布參數(shù)的Jeffreys無信息先驗(yàn)分布,提出可變作用標(biāo)準(zhǔn)值和頻遇值的貝葉斯推斷方法.應(yīng)用研究結(jié)果表明:貝葉斯推斷方法較線性回歸推斷方法簡便且應(yīng)用范圍廣,在標(biāo)準(zhǔn)差已知的情況,可給出更優(yōu)的推斷結(jié)果;無參數(shù)信息時(shí),對(duì)于標(biāo)準(zhǔn)值和頻遇值的保證率不低于0.90的情況,貝葉斯推斷方法均具有較好的精度;相對(duì)于其他可信水平,可信水平0.75時(shí)推斷結(jié)果更接近真值,因此,建議取可信水平為0.75.

荷載;代表值;貝葉斯方法;統(tǒng)計(jì)推斷;小樣本容量

可變作用的代表值包括標(biāo)準(zhǔn)值、頻遇值和準(zhǔn)永久值等[1],對(duì)它們的推斷是建立結(jié)構(gòu)設(shè)計(jì)與評(píng)定方法的基礎(chǔ).當(dāng)樣本容量(測(cè)試數(shù)據(jù)數(shù)量)充足時(shí),一般可采用矩法、極大似然法等經(jīng)典統(tǒng)計(jì)學(xué)的推斷方法,但實(shí)際中的測(cè)試數(shù)據(jù)不總是充足的,如仍采用目前經(jīng)典統(tǒng)計(jì)學(xué)的方法,往往導(dǎo)致代表值的推斷結(jié)果因統(tǒng)計(jì)不確定性的影響而偏于冒進(jìn),這時(shí)更合理的選擇是采用小樣本的推斷方法.

可變作用的任意時(shí)點(diǎn)值通常服從極大值Ⅰ型分布,其代表值一般可表達(dá)為該分布的某一分位值[2-3].由于極大值Ⅰ型分布與極小值Ⅰ型分布同屬于極值參數(shù)分布族,可相互轉(zhuǎn)換[4],因此,對(duì)可變作用代表值的推斷可借用目前極小值Ⅰ型分布分位值的線性回歸推斷方法[5].該法可在不同的置信度下考慮統(tǒng)計(jì)不確定性的影響,用于小樣本的場(chǎng)合,在機(jī)械、電子等領(lǐng)域產(chǎn)品壽命的推斷中有著廣泛的應(yīng)用[6-7].但是,采用該法推斷時(shí)需查表確定大量系數(shù)的數(shù)值,應(yīng)用不便,且目前的數(shù)值表并不完全滿足可變作用代表值推斷的需要[8],而建立新的數(shù)值表又需進(jìn)行繁瑣的數(shù)值計(jì)算或蒙特卡洛模擬,存在較大的困難.因此,這種方法不能很好地解決小樣本條件下可變作用代表值的推斷問題.

小樣本條件下,貝葉斯法也是一種可選擇的推斷方法[9],除了考慮經(jīng)典統(tǒng)計(jì)學(xué)中的總體信息和樣本信息,它還考慮了關(guān)于總體的先驗(yàn)信息[10].如果對(duì)總體的分布參數(shù)采用無信息的先驗(yàn)分布,降低主觀因素的影響,則貝葉斯法同樣可在小樣本條件下給出合理的推斷結(jié)果,且較線性回歸推斷方法簡便,也便于在更廣的范圍內(nèi)使用[11-12].但是,目前的貝葉斯推斷方法主要適用于總體服從二項(xiàng)分布、泊松分布、正態(tài)分布、指數(shù)分布等的場(chǎng)合,對(duì)于服從極大值Ⅰ型分布的總體,目前尚無合適的推斷方法[13].因此,小樣本條件下對(duì)可變作用代表值的推斷需采用新的貝葉斯方法.

本文將首先建立和評(píng)價(jià)可變作用代表值的線性回歸推斷方法,并以此作為檢驗(yàn)其他推斷方法精度的基準(zhǔn)方法;然后,重點(diǎn)研究和建立可變作用代表值的貝葉斯推斷方法,并通過對(duì)比分析評(píng)價(jià)其優(yōu)劣.由于可變作用的準(zhǔn)永久值一般為可變作用任意時(shí)點(diǎn)分布的均值[1],受統(tǒng)計(jì)不確定性的影響較小,可采用經(jīng)典統(tǒng)計(jì)學(xué)的方法推斷,因此,論文將主要研究可變作用標(biāo)準(zhǔn)值和頻遇值的推斷問題.

1 線性回歸推斷方法

一般假定可變作用的任意時(shí)點(diǎn)值X服從極大值Ⅰ型分布[2],其概率密度函數(shù)為

式中:μ、α為分布參數(shù),-∞<μ<∞,0<α<∞.

可變作用的標(biāo)準(zhǔn)值或頻遇值可表示為X的下側(cè)p分位值xp[1],它滿足

式中:p為標(biāo)準(zhǔn)值或頻遇值的保證率,通常取[0,1]內(nèi)較大的數(shù)值;

為借用極小值Ⅰ型分布分位值的線性回歸推斷方法,設(shè)X的n個(gè)次序統(tǒng)計(jì)量(由小到大排列的樣本)為X(1),X(2),…,X(n),相應(yīng)的測(cè)試值為x(1),x(2),…,x(n),并令

則Y服從參數(shù)為-μ、α的極小值Ⅰ型分布,其次序統(tǒng)計(jì)量和上側(cè)p分位值分別為

根據(jù)極小值Ⅰ型分布分位值的線性回歸推斷方法[14],在置信度C下yp的估計(jì)量為

DI(n,n,j)、CI(n,n,j)為僅與樣本容量n、樣本序位j有關(guān)的系數(shù);

vp,C為統(tǒng)計(jì)量的下側(cè)C分位值.

DI(n,n,j)、CI(n,n,j)、vp,C均可查表確定[8].由于

故置信度C下可變作用標(biāo)準(zhǔn)值或頻遇值xp的估計(jì)量為

這時(shí)根據(jù)測(cè)試值x(1),x(2),…,x(n),便可最終確定xp的估計(jì)值.

該線性回歸推斷方法不僅考慮了樣本容量,還考慮了樣本序位,更充分地利用了樣本提供的信息,并可在不同的置信度下考慮統(tǒng)計(jì)不確定性的影響,適用于小樣本場(chǎng)合.但是,該法在應(yīng)用中需查表確定DI(n,n,j)、CI(n,n,j)、vp,C的數(shù)值,查閱過程非常繁瑣,且目前vp,C數(shù)值表中僅提供了p=0.90,0.95,0.99和n≤25時(shí)的數(shù)值[8],并不完全滿足可變作用標(biāo)準(zhǔn)值和頻遇值推斷的需要,而建立新的數(shù)值表又需進(jìn)行大量的數(shù)值模擬計(jì)算,存在較大的困難.另外,目前僅建立了無參數(shù)信息時(shí)的線性回歸推斷方法,對(duì)于標(biāo)準(zhǔn)差已知的場(chǎng)合,尚無相應(yīng)的方法[14].

2 貝葉斯推斷方法

2.1 無參數(shù)信息時(shí)的情況

設(shè)可變作用任意時(shí)點(diǎn)值X的n個(gè)樣本測(cè)試值為x1,x2,…,xn,樣本均值和樣本標(biāo)準(zhǔn)差分別為xˉ和s,則它們發(fā)生的聯(lián)合概率密度函數(shù)為

取未知參數(shù)μ、α的先驗(yàn)分布為Jeffreys無信息先驗(yàn)分布[13],即

則μ、α的聯(lián)合后驗(yàn)分布為

按式(3)進(jìn)行變量代換后,可得xp、α的聯(lián)合后驗(yàn)分布為

為簡化問題,將式(15)中的指數(shù)函數(shù)近似表達(dá)為泰勒級(jí)數(shù)中的線性項(xiàng),即

式中:歐拉常數(shù)CE≈0.577 2,展開點(diǎn)為(xi-xp+kα-CEα)/α=0.

由于X的均值

在展開點(diǎn)處有xi-μX=0,將式(16)代入式(15)并求和可獲得相對(duì)精確的結(jié)果.這時(shí)xp、α的聯(lián)合后驗(yàn)分布可整理為

將式(17)中的最后一項(xiàng)按泰勒級(jí)數(shù)展開,即

并令

則U,V的聯(lián)合分布為

式(21)中最后一個(gè)分式為自由度n+m的χ2分布的概率密度函數(shù)[4].對(duì)v積分后,可得U的邊緣分布

故U服從自由度n-1、參數(shù)λ的非中心t分布[4],

按區(qū)間估計(jì)法中的上限估計(jì)值,最終可得可變作用標(biāo)準(zhǔn)值或頻遇值的估計(jì)值,即

式中:t(n-1,λ,1-C)為自由度n-1、參數(shù)λ的非中心t分布的上側(cè)1-C分位值,C為可信水平;

由于目前t(n-1,λ,1-C)的數(shù)值表中所考慮的參數(shù)λ不能滿足標(biāo)準(zhǔn)值和頻遇值推斷的需要,可采用下列方法近似計(jì)算t(n-1,λ,1-C)的值[15],

式中:z1-C為標(biāo)準(zhǔn)正態(tài)分布的上側(cè)1-C分位值.

2.2 標(biāo)準(zhǔn)差σX已知的情況

可變作用任意時(shí)點(diǎn)值X的標(biāo)準(zhǔn)差σX已知時(shí),分布參數(shù)同樣取未知參數(shù)μ的先驗(yàn)分布為Jeffreys無信息先驗(yàn)分布[13],即

則按貝葉斯公式變量代換后,可得分位值xp的后驗(yàn)分布為

則U的分布為

故U服從參數(shù)為n的標(biāo)準(zhǔn)伽馬分布Ga(n,1)[4].這時(shí)按區(qū)間估計(jì)法中的上限估計(jì)值,可得可變作用標(biāo)準(zhǔn)值或頻遇值的估計(jì)值,即

式中:γ(n,1,C)為標(biāo)準(zhǔn)伽馬分布Ga(n,1)的下側(cè)C分位值,C為可信水平;k2=γ(n,1,C)/n.

這時(shí)的貝葉斯推斷方法在建立過程中并未采用近似手段.

3 對(duì)比分析

可變作用標(biāo)準(zhǔn)值和頻遇值的保證率p一般均不低于0.90[1-2],如對(duì)風(fēng)、雪荷載,p=0.98(標(biāo)準(zhǔn)值)和0.90(頻遇值)[16],因此,僅針對(duì)p≥0.90的情況對(duì)比分析貝葉斯推斷方法的精度.設(shè)可變作用任意時(shí)點(diǎn)值X的10個(gè)次序統(tǒng)計(jì)量的測(cè)試值為x(1),x(2),…,x(10)(見表1),其統(tǒng)計(jì)結(jié)果為

為便于對(duì)比分析,設(shè)xp的保證率p為目前vp,C數(shù)值表中的0.90、0.95和0.99,且σX已知時(shí)取σX=s.

表2列舉了不同可信水平C和保證率p下xp的推斷結(jié)果,推斷過程中DI(10,10,j)、CI(10,10,j)的數(shù)值及相關(guān)計(jì)算結(jié)果見表1.

表1 樣本測(cè)試值和DI(10,10,j)、CI(10,10,j)的值Tab.1 Test values of sample and values of DI(10,10,j),CI(10,10,j) kN/m

表2 xp的推斷結(jié)果Tab.2 The inference result of xp

由表2中的結(jié)果可知:矩法的推斷結(jié)果最低,偏于冒進(jìn),且C越高,與其它方法的差異越大,這主要是因?yàn)槲闯浞挚紤]小樣本條件下統(tǒng)計(jì)不確定性的影響;無參數(shù)信息時(shí),采用近似展開式的貝葉斯法的推斷結(jié)果精度較高,且p、C越高,相對(duì)誤差越小,并偏于保守,適用于任意保證率不低于0.90的情況;σX已知時(shí),貝葉斯法的推斷結(jié)果因獲知更多的參數(shù)信息而優(yōu)于無參數(shù)信息時(shí)的結(jié)果,包括線性回歸推斷方法的結(jié)果,且p、C越高,優(yōu)勢(shì)越明顯,因此,當(dāng)σX未知時(shí),可對(duì)σX取一個(gè)偏大的值,按σX已知時(shí)的方法推斷.研究中還對(duì)比分析了樣本容量分別為5和20時(shí),不同p、C值的推斷結(jié)果,可得到同樣的結(jié)論.

無參數(shù)信息時(shí)貝葉斯推斷方法的誤差主要來源于式(16)所采用的近似方法.按該式可將μ、α的聯(lián)合后驗(yàn)分布整理為

式(32)實(shí)際上為正態(tài)分布N(μ′,σ′2)的分布參數(shù)μ′、σ′的聯(lián)合后驗(yàn)分布,其中

這相當(dāng)于以正態(tài)分布N(μ′,σ′2)替代了原先的極大值Ⅰ型分布.由于極大值Ⅰ型分布、正態(tài)分布在概率密度函數(shù)曲線的右端更為相似,因此,在推斷保證率p較高的分位值xp時(shí)可獲得相對(duì)精確的結(jié)果.

由表2中的結(jié)果還可知,C≥0.90時(shí),貝葉斯法的推斷結(jié)果要明顯高于矩法的推斷結(jié)果,且對(duì)C的變化敏感;相對(duì)而言,C=0.60,0.75時(shí)的推斷結(jié)果較適宜,且兩者的數(shù)值差別不大,但按C=0.75推斷時(shí)可更充分地考慮統(tǒng)計(jì)不確定性的影響,相對(duì)誤差也更小,因此,相對(duì)其他可信水平,按C=0.75推斷可獲得更合適的結(jié)果,建議取

為便于應(yīng)用,表3列舉了可信水平C=0.75時(shí)的k1、k2值.

表3 k1、k2數(shù)值表(C=0.75)Tab.3 Numerical tables of k1and k2(C=0.75)

4 結(jié) 論

(1)測(cè)試數(shù)據(jù)不足時(shí),對(duì)可變作用準(zhǔn)永久值的推斷可采用經(jīng)典統(tǒng)計(jì)學(xué)的方法,但推斷標(biāo)準(zhǔn)值和頻遇值時(shí),經(jīng)典統(tǒng)計(jì)學(xué)方法的推斷結(jié)果偏于冒進(jìn),它未充分考慮小樣本條件下統(tǒng)計(jì)不確定性的影響.

(2)小樣本條件下,可借用極小值Ⅰ型分布分位值的線性回歸推斷方法,推斷可變作用的標(biāo)準(zhǔn)值和頻遇值,但需查表確定大量系數(shù)的數(shù)值,應(yīng)用不便,且目前的數(shù)值表不能完全滿足推斷的需要.

(3)用本文中的貝葉斯推斷方法,同樣可在小樣本條件下推斷可變作用的標(biāo)準(zhǔn)值和頻遇值,但較線性回歸推斷方法簡便,適用范圍更廣且在標(biāo)準(zhǔn)差已知的情況下,可給出更精確的推斷結(jié)果.

(4)無參數(shù)信息時(shí),對(duì)于保證率不低于0.90的情況,文中貝葉斯推斷方法得到的可變作用標(biāo)準(zhǔn)值和頻遇值具有較高的精度.

(5)相對(duì)于其他可信水平,貝葉斯推斷方法在可信水平為0.75時(shí)的推斷結(jié)果要更為合適,建議取可信水平為0.75.

[1] 中國建筑科學(xué)研究院.GB50153—2008工程結(jié)構(gòu)可靠性設(shè)計(jì)統(tǒng)一標(biāo)準(zhǔn)[S].北京:中國建筑工業(yè)出版社,2009.

[2] 中國建筑科學(xué)研究院.GBJ68—84建筑結(jié)構(gòu)設(shè)計(jì)統(tǒng)一標(biāo)準(zhǔn)[S].北京:中國建筑工業(yè)出版社,1984.

[3] 馮云芬,貢金鑫,王建超.樓面活荷載、風(fēng)荷載的頻遇值和準(zhǔn)永久值的確定[J].工業(yè)建筑,2012,42(7):74-78.FENGYunfen,GONGJinxin,WANGJianchao.Determination of frequent value and quasi-permanent value of floor live load and wind load[J].Industrial Construction,2012,42(7):74-78.

[4] 茆詩松.統(tǒng)計(jì)手冊(cè)[M].北京:科學(xué)出版社,2003:212-223.

[5] 周源泉.質(zhì)量可靠性增長與評(píng)定方法[M].北京:北京航空航天大學(xué)出版社,1997:151-174.

[6] 黃炎生,鄧浩,李濤.考慮驗(yàn)證荷載的既有結(jié)構(gòu)可靠度計(jì)算[J].華南理工大學(xué)學(xué)報(bào):自然科學(xué)版,2008,36(10):12-15.HUANG Yansheng,DENG Hao,LI Tao.Reliability calculation of existing structure under proof load[J].JournalofSouthChinaUniversityofTechnology:Natural Science Edition,2008,36(10):12-15.

[7] AGGARWALA R.Progressive interval censoring:some mathematicalresultsapplicationtoinference[J].Communications in Statistics:Theory and Methods,2001,30:1921-1935.

[8] 第四機(jī)械工業(yè)部標(biāo)準(zhǔn)化研究所.可靠性試驗(yàn)用表[M].北京:國防工業(yè)出版社,1979:1-87.

[9] CHOPIN N,LELIèVRE T,STOLTZ G.Free energy methods for Bayesian inference:efficient exploration of univariateGaussianmixtureposteriors[J].Stat.Comput,2012,22:897-916.

[10] TSIONAS E G.Bayesian inference for multivariate gamma distributions[J].Statistics and Computing,2004,14(3):223-233.

[11] PENG Xiuyun,YAN Zaizai.Bayesian estimation for generalizedexponentialdistributionbasedon progressivetype-Ⅰintervalcensoring[J].Acta Mathematicae Applicatae Sinica,2013,29(2):391-402.

[12] LIN YJ,LIOYL.Bayesianinferenceunder progressive type-Ⅰinterval censoring[J].Journal of Applied Statistics,2012,39(1):1811-1824.

[13] 茆詩松.貝葉斯統(tǒng)計(jì)[M].北京:中國統(tǒng)計(jì)出版社,1999:75-79.

[14] 戴樹森,費(fèi)鶴良.可靠性試驗(yàn)及其統(tǒng)計(jì)分析(上冊(cè))[M].北京:國防工業(yè)出版社出版,1983:617-628.

[15] 姚繼濤.基于不確定性推理的既有結(jié)構(gòu)可靠性評(píng)定[M].北京:科學(xué)出版社,2011:50-57.

[16] 中國建筑科學(xué)研究院.GB50009—2001建筑結(jié)構(gòu)荷載規(guī)范[S].北京:中國建筑工業(yè)出版社,2001.

(中文編輯:秦萍玲 英文編輯:蘭俊思)

Bayesian Methods for Inferring Representative Values of Variable Actions in Small Sample Situations

YAO Jitao, WANG Xudong
(School of Civil Engineering,Xi'an University of Architecture and Technology,Xi'an 710055,China)

There is a need sometimes in engineering to infer the representative values of variable actions in cases when the test data is not sufficient,but the classical statistics methods do not take into account of the succedent influences of statistical uncertainty,and the inferred results are always on the aggressive side.In order to overcome the above shortcomings,applying the current linear regression estimation for inferring the fractiles of type I minimum distribution,the linear regression estimation of representative values of variable actions is proposed as a reference method for the precision inspection of other inference methods.According to Bayesian theory,applying Jeffreys non-informative prior distribution,a Bayesian method for inferring characteristic and frequent values of variable actions is put forward.The results show that the Bayesian inference method is more convenient than the liner regression estimation,and easy to use in a broader context.The method can give more advantageous results when standard deviation is known.If no parameter information is available,the Bayesian inference method has better precision when the guarantee rates of characteristic and frequent values are no less than 0.90.Compared with other confidence degrees,the inferred results are more close to the true value when the confidence degree equals 0.75;therefore,it is recommended to take a confidence degree of 0.75.

load;characteristic value;Bayesian approach;statistical inference;small sample

TU312.1

:A

0258-2724(2014)06-0995-07

10.3969/j.issn.0258-2724.2014.06.010

2013-09-29

國家自然科學(xué)基金資助項(xiàng)目(50678143,51278401)

姚繼濤(1965-),男,教授,博士,研究方向?yàn)榻Y(jié)構(gòu)可靠性理論,電話:029-82202359,E-mail:yaojitao@163.com

姚繼濤,王旭東.小樣本條件下可變作用代表值的貝葉斯推斷方法[J].西南交通大學(xué)學(xué)報(bào),2014,49(6):995-1001.

猜你喜歡
方法
中醫(yī)特有的急救方法
中老年保健(2021年9期)2021-08-24 03:52:04
高中數(shù)學(xué)教學(xué)改革的方法
化學(xué)反應(yīng)多變幻 “虛擬”方法幫大忙
變快的方法
兒童繪本(2020年5期)2020-04-07 17:46:30
學(xué)習(xí)方法
用對(duì)方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
最有效的簡單方法
山東青年(2016年1期)2016-02-28 14:25:23
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 91小视频版在线观看www| 在线免费观看AV| 中文字幕天无码久久精品视频免费| 一级不卡毛片| 久久一色本道亚洲| 亚洲精品黄| 亚洲不卡无码av中文字幕| 久久久久九九精品影院 | 国产免费黄| 国产在线98福利播放视频免费| 波多野结衣一二三| 国产日韩av在线播放| 亚洲欧美日韩精品专区| 国产精品专区第一页在线观看| 欧美视频免费一区二区三区| 久久黄色免费电影| 尤物成AV人片在线观看| 亚洲中文字幕久久无码精品A| 影音先锋丝袜制服| 亚洲一区毛片| 色综合成人| 亚洲精品第五页| 五月六月伊人狠狠丁香网| 免费在线不卡视频| 亚洲天堂免费| 手机成人午夜在线视频| 久久久久人妻精品一区三寸蜜桃| 最新亚洲人成无码网站欣赏网 | 国产成人高清在线精品| 亚卅精品无码久久毛片乌克兰| 免费看美女毛片| 日本成人不卡视频| 久久久久久尹人网香蕉 | 欧洲极品无码一区二区三区| 伊人久久精品亚洲午夜| 91精品国产一区自在线拍| 免费一级毛片在线播放傲雪网| 婷婷丁香色| 国产精品久久久久久影院| 国产在线欧美| 午夜在线不卡| 日韩欧美在线观看| 在线免费看片a| 精品无码一区二区在线观看| 国产麻豆另类AV| 亚洲经典在线中文字幕| 日韩少妇激情一区二区| 国内精品小视频福利网址| 午夜欧美理论2019理论| 1024你懂的国产精品| 在线观看热码亚洲av每日更新| 精品无码国产自产野外拍在线| 中文字幕在线看视频一区二区三区| 99ri精品视频在线观看播放| 久久久精品国产SM调教网站| 午夜福利视频一区| 国产在线观看一区二区三区| 久久精品aⅴ无码中文字幕| 亚洲精品成人片在线观看| 超碰91免费人妻| 成人无码一区二区三区视频在线观看| 激情视频综合网| 国产成+人+综合+亚洲欧美| 婷婷激情亚洲| 亚洲中文久久精品无玛| 婷婷色狠狠干| 中国一级特黄大片在线观看| 国产色爱av资源综合区| 中文无码毛片又爽又刺激| 2024av在线无码中文最新| 国产欧美另类| 国产丰满大乳无码免费播放| 欧美三级自拍| 在线va视频| 国产欧美精品一区二区| 国产亚洲精品97在线观看| 亚洲AV无码一区二区三区牲色| 欧美翘臀一区二区三区| 福利在线一区| 老司机午夜精品视频你懂的| 粗大猛烈进出高潮视频无码| 日韩免费毛片|