日韩av黄I国产麻豆传媒I国产91av视频在线观看I日韩一区二区三区在线看I美女国产在线I麻豆视频国产在线观看I成人黄色短片

歡迎訪問 生活随笔!

生活随笔

當前位置: 首頁 > 编程资源 > 编程问答 >内容正文

编程问答

MATLAB统计与回归

發(fā)布時間:2023/12/10 编程问答 40 豆豆
生活随笔 收集整理的這篇文章主要介紹了 MATLAB统计与回归 小編覺得挺不錯的,現(xiàn)在分享給大家,幫大家做個參考.

?

11.1 前言

統(tǒng)計的技巧與資料分析常常形影不離。一般統(tǒng)計使用加法、累加法、平均值,中間值等等,由於處理的對象是矩陣資料,故其基本統(tǒng)計之技巧已經(jīng)廣為應用,其觀念也會在正常之運作中出現(xiàn)。統(tǒng)計學中比較特殊應用者為機率、亂數(shù)、常態(tài)分配等,而配合應用者為其相關(guān)之圖表。

在MATLAB中,有一個統(tǒng)計學工具箱,內(nèi)藏各種統(tǒng)計學上需要應用的指令,可以執(zhí)行上述與統(tǒng)計學有關(guān)之內(nèi)容。這些相關(guān)的指令大部份以M-檔案組成,所以可利用type 這個功能檢視其內(nèi)容。甚至可以更改其檔案名稱與內(nèi)容,增加自己需要的功能,使其成為新的指令。此外,有些指令尚搭配繪圖介面,因而可以在繪圖模式下,進行資料與圖之適配,形成具體的方程式或?qū)嶒炇?#xff0c;以供未來研究者使用。

統(tǒng)計工具箱中提供約二百餘個指令檔案,其中對機率分配方面則提供廿餘種機率型態(tài),每種均有其相關(guān)的函數(shù),諸如:

  • 機率密度函數(shù)(pdf)
  • 累積分佈函數(shù)(cdf)
  • 累積分佈函數(shù)之反函數(shù)
  • 亂數(shù)產(chǎn)生器
  • 均數(shù)與變異數(shù)

?

?

11.2.1平均值MEAN


統(tǒng)計學上對資料處理常用趨中的處理,求取均值或中間值等,均會取中的特徵。求取一個矩陣或向量之平均值時可用指令MEAN,其格式如下:

? M=mean(A,dim)

若A為向量,其結(jié)果M為單一值,亦即向量中各元素之平均;若A為矩陣,則結(jié)果M為一列向量,其中元素為各行之平均值。dim為方向性參數(shù),其預設值為1,表示結(jié)果係行向平均,故M為列向量;若dim=2,則係沿列向平均,結(jié)果M為行向量。例如:

A = [1 2 3; 3 3 6; 4 6 8; 4 7 7]
A =
??? 1???? 2???? 3
??? 3???? 3???? 6
??? 4???? 6???? 8
??? 4???? 7???? 7

M1=mean(A), M2=mean(A,2)

M1 =
?? 3.0000??? 4.5000??? 6.0000

M2 =
??? 2
??? 4
??? 6
??? 6

?

?

11.2.2中間值median


中間值亦可利用平均值的指令型式求得,其正式指令名稱為median。但其求得之值若非正好中間值,則會以接近中間值之兩值加以平均,其結(jié)果與mean之平均值仍有不同,下面以前述之B矩陣為例,比較median 與mean兩者執(zhí)行後之不同處。

M1=median(B),M2=mean(B)?

M1 =
??? 5.5000??? 5.5000??? 5.0000??? 6.0000
M2 =
??? 5.0000??? 5.5000??? 4.7500??? 6.0000?

其他的平均值則有:
幾何平均(geomean):各元素乘積再開總數(shù)之次方。中間有零值時,其結(jié)果為零。
調(diào)諧平均(harmmean):各元素倒數(shù)和之倒數(shù)乘以總數(shù)。
修剪平均(trimmean):去頭去尾再平均的方式,其頭尾部份為第二參數(shù)(%)之一半比例。

下面的例子為這些平均值的比較:

X=1:2:20
locate=[geomean(X) harmmean(X) mean(X) median(X) trimmean(X,25)]?

X =
???? 1???? 3???? 5???? 7???? 9??? 11??? 13??? 15??? 17??? 19
locate =
??? 7.6139??? 4.6877?? 10.0000?? 10.0000?? 10.0000?

由上述這些值看來,有些的結(jié)果相同,有些則因其分離的情況而有異。中間值及修剪平均值是不管偏異的大小,一般為避免有極端的數(shù)值的影響,採用這種平均法。

?

?

11.2.3最大值與最小值 MIN & MAX


要求一矩陣或向量中元素之最大與最小值時,其指令之型式如下:

? C = max|min(A)
? C = max|min(A,B)
? C = max|min(A,[],dim)
? [C,I] = max(...)

若A為向量,其結(jié)果C為單一值,亦即向量中各元素之最大或最小;若A為矩陣,則結(jié)果C為一列向量,其中元素為各行之最大或最小。dim為方向性參數(shù),其預設值為1,表示結(jié)果係行向取得最大或最小,故C為列向量;若dim=2,則係沿列向操作,結(jié)果M為行向量。注意要dim之參數(shù)時,需加在第三位置。此外,在輸出項中,I表示最大或最小元素之位置,不過此項功能僅求最大值時適用。例如:

A= [7 2 3 4; 2 4 5 6; 4 6 8 5; 6 7 6 1];

[C,Index]=max(A)?

A =
???? 7???? 2???? 3???? 4
???? 2???? 4???? 5???? 6
???? 4???? 6???? 8???? 5
???? 6???? 7???? 6???? 1
C =
???? 7???? 7???? 8???? 6
Index =
???? 1???? 4???? 3???? 2?

在這個指令中,比較特殊的是兩個矩陣A與B之最大或最小比較,其結(jié)果C應為與A或B相同的矩陣,但比較A與B中對應元素間之最大與最小。例如:


B=[1 8 7 6;4 6 2 8;7 5 6 4;8 3 4 6]
[Cab]=min(A,B)?

B =
???? 1???? 8???? 7???? 6
???? 4???? 6???? 2???? 8
???? 7???? 5???? 6???? 4
???? 8???? 3???? 4???? 6
Cab =
???? 1???? 2???? 3???? 4
???? 2???? 4???? 2???? 6
???? 4???? 5???? 6???? 4
???? 6???? 3???? 4???? 1?

?

?

11.2.4變方值VAR


變方值為各樣本品與平均值差之平方和。又稱為變值常態(tài)檢定公式,其指令型式為var(X,1)。其計算之公式如下:

VAR(X)=[(X1-r)2+(X2-r)2+...+(Xn-r)2]/(n-1)

式中,n為其總項目,若n>1,則採用n-1以獲得無偏差之變方值。相關(guān)指令如下:

  Y=var(X)
  Y=var(X,1)
  Y=var(X, w)
  Y=var(X, w, dim)

Y=var(X)為執(zhí)行上式計算之無偏差變方,Y=var(X,1)則採用n為除數(shù)。另w則為加權(quán)值向量,此權(quán)值必須為正值,且長度與X相同。

若X為矩陣時,通常預設為行向計算,但可以利用dim=2參數(shù)改為以列向為計算基礎,其結(jié)果為行向量。var指令會將其元素除以總和,因此權(quán)值總和為一。若w值為零,其結(jié)果如var(X);若為1則如var(X,1)。

範例



>> x=1:9
x =

??? 1???? 2???? 3???? 4???? 5???? 6???? 7???? 8???? 9

>> var(x) %除以(n-1)
ans =

?? 7.5000

>> var(x,1) %除以n

ans =
?? 6.6667

>> w=[1 1 1 1 1 0.5 0.5 0.5 0.5] %設為權(quán)值,向量與x同
w =

?? 1.0000??? 1.0000??? 1.0000??? 1.0000??? 1.0000??? 0.5000??? 0.5000??? 0.5000??? 0.5000

>> var(x,w)
ans =

?? 5.9184

?

?

11.2.5標準差STD


標準誤差為各樣本品與平均值間之常態(tài)差,其值實際上為上述變方var執(zhí)行結(jié)果之開方值,其計算公式如下:

?

11.2.6 共方差COV


共方差為兩向量之觀察值與其平均差之乘積和,其計算之函數(shù)定義如下:


cov(x1,x2)=E[(x1-u1)(x2-u2)]
 其中 ui=Exi

根據(jù)上式,其相關(guān)指令格式為:

C = cov(X)
C = cov(x,y)

在COV之指令,若 X為向量,其回應值應為變方值。若其為矩陣,則各列為觀察值,各行則成為變數(shù),而COV(X)則為共方矩陣。其對角線元素 DIAG(COV(X))即為每行之變方差向量。若將之排序後,即SQRT(DIAG(COV(X))),其結(jié)果為標準差之向量。以下為例:

>> x=[4??? -2???? 1;? 9???? 5???? 7]
C1=cov(x)
M=diag(cov(x))
sort(diag(cov(x)))?

x =
???? 4??? -2???? 1
???? 9???? 5???? 7
C1 =
?? 12.5000?? 17.5000?? 15.0000
?? 17.5000?? 24.5000?? 21.0000
?? 15.0000?? 21.0000?? 18.0000
M =
?? 12.5000
?? 24.5000
?? 18.0000
ans =
?? 12.5000
?? 18.0000
?? 24.5000?


COV(X,Y)則是求取 X與 Y兩等長度之向量之共方差, X與 Y兩向量即使為列向量亦會自動改為行向量,其效果等於COV( [X(:) Y(:)] )。這兩個指令均設法加以常態(tài)化,故母數(shù)除以N-1,以消除偏差。若要維持使用N為母數(shù),則可增加參數(shù)1,即採用 COV(X,1) 或 COV(X,Y,1)指令之型式。

?

?

11.2.7相關(guān)係數(shù) corrcoef


兩個變數(shù)相關(guān)性可由相關(guān)係數(shù)求得。其指令型式如下:

R = corrcoef(X)
R = corrcoef(x,y)
[R,P]=corrcoef(...)
[R,P,RLO,RUP]=corrcoef(...)
[...]=corrcoef(...,'param1',val1,'param2',val2,...)


基本上,R=CORRCOEF(X)在於計算一個R矩陣,其內(nèi)有X陣列行間之相關(guān)係數(shù)。而CORRCOEF(X,Y)則計算 X 與 Y兩行向量之相關(guān)係數(shù),其意義與CORRCOEF([X Y])相同。

假設 C為共方矩陣,且 C = COV(X),則R=CORRCOEF(X)之定義為:

?? R(i,j) = C(i,j)/sqrt{C(i,i)C(j,j)}

輸入項中除XY等資料矩陣外,尚可輸入其他特定變數(shù)與常數(shù)。這些可以用 'PARAM1',VAL1成對表示,其項目包括:

???? 參數(shù)值:
????? 'alpha'??? 顯著水準,預設值為0.05(即95%信任區(qū)間)
????? 'rows'???? 使用 'all' (預設值)表示使用所有列值;
???????????????? 'complete'表示使用沒有含NaN 值之列;
???????????????? 'pairwise'表示計算R(i,j)時使用不含
???????????????? NaN值之 i行或 j行。

輸出值中, P表示檢驗無關(guān)係假設之P值矩陣。每一個P值代表隨機可以觀察得到之最大值域。若 P(i,j)值很小,例如小於 0.05,則R(i,j) 之關(guān)係甚為顯著。

此外,有RLO與RUP代表95%信任水準之下限與上限矩陣,其大小與R相同。

例一



>> x=1:5
x =

???? 1???? 2???? 3???? 4???? 5

>> y=x.^3
y =

???? 1???? 8??? 27??? 64?? 125

>> r=corrcoef(x,y)
r =

??? 1.0000??? 0.9431
??? 0.9431??? 1.0000

答案中r之值愈接近於1,其相關(guān)性愈高。此例中,對角線為自己對自己(即x對x;y對y)故其相關(guān)性為1,其餘x對y或y對x,兩者相關(guān)性一樣,其數(shù)值為0.9431,也相當高。

例二


利用常態(tài)分配亂數(shù)指令randn產(chǎn)生30X4大小之資料,開始時先利用第四行建立與其他行間之關(guān)係,以橫向加總於第四行。其後以corrcoef求相關(guān)係數(shù)r及機率p。就機率而言,p值愈小,表示兩者之相異性更強,其結(jié)果可利用find指令找出小於0.05以下之機率項目。

?????? x = randn(20,4);?????? % uncorrelated data
?????? x(:,4) = sum(x,2)???? % introduce correlation
?????? [r,p] = corrcoef(x)??
?????????????? % compute sample correlation and p-values
?????? [i,j] = find(p<0.05);
?????????????? % find significant correlations [i,j]

x =
?? 0.0828?? -0.5703?? -0.0716?? -0.5850
?? 0.7662?? -1.4986?? -2.4146?? -4.2576
?? 2.2368?? -0.0503?? -0.6943??? 2.2430
?? 0.3269??? 0.5530?? -1.3914?? -0.0113
?? 0.8633??? 0.0835??? 0.3296??? 0.7592
?? 0.6794??? 1.5775??? 0.5985??? 2.2962
?? 0.5548?? -0.3308??? 0.1472?? -0.3822
?? 1.0016??? 0.7952?? -0.1014??? 2.6212
?? 1.2594?? -0.7848?? -2.6350?? -2.4089
?? 0.0442?? -1.2631??? 0.0281?? -1.3408
? -0.3141??? 0.6667?? -0.8763?? -1.7822
?? 0.2267?? -1.3926?? -0.2655?? -1.1188
?? 0.9967?? -1.3006?? -0.3276??? 2.0588
?? 1.2159?? -0.6050?? -1.1582?? -0.2577
? -0.5427?? -1.4886??? 0.5801?? -2.8740
?? 0.9122??? 0.5585??? 0.2398??? 1.9573
? -0.1721?? -0.2774?? -0.3509?? -2.2362
? -0.3360?? -1.2937??? 0.8921?? -0.5890
?? 0.5415?? -0.8884??? 1.5783?? -0.4617
?? 0.9321?? -0.9865?? -1.1082?? -0.4434
r =
?? 1.0000??? 0.1950?? -0.3475??? 0.5143
?? 0.1950??? 1.0000??? 0.0929??? 0.5785
? -0.3475??? 0.0929??? 1.0000??? 0.3822
?? 0.5143??? 0.5785??? 0.3822??? 1.0000
p =
?? 1.0000??? 0.4100??? 0.1333??? 0.0203
?? 0.4100??? 1.0000??? 0.6969??? 0.0075
?? 0.1333??? 0.6969??? 1.0000??? 0.0963
?? 0.0203??? 0.0075??? 0.0963??? 1.0000
ans =
??? 4???? 1
??? 4???? 2
??? 1???? 4
??? 2???? 4

?

?

11.2.8群組函數(shù)grpstats


前面討論到之平均值求法,通常應用於整個陣列之值,若要應用到比較複雜的分組平均問題,則必須使用不同的函數(shù)才能達成。此項指令之格式如下:

means = grpstats(X, group)
[means, sem, counts, name] = grpstats(X, group, whichstats)
grpstats(x, group, alpha)

輸入?yún)?shù)中X為求平均值之對象,X可為多行,其平均結(jié)果也會多行。group則為與X同列長之陣列,可能由多項分組之向量組成,其內(nèi)容可為字串列或細胞陣列之文字,如{G1 G2 G3}。若X中之元素同屬分組中之一項,則其平均值會出現(xiàn)在該項下。

在輸出項中,第一項means為群組平均,sem為組內(nèi)標準差,counts為各組之項數(shù),name則為各組之名稱。上述項目並非一成不變,亦可以在輸入?yún)?shù)whichstats內(nèi)依自己之需要進行設定,這個設定有特定的名稱,其名稱必須使用細胞陣列。項目包括:
  • 'mean' 組平均
  • 'sem' 組標準差
  • 'numel' 各組之數(shù)目
  • 'gname' 各組之名稱
  • 'std' 標準差
  • 'var' 變異值
  • 'meanci' 平均值之95%上下範圍
  • 'predci' 新值預測之95%信任範圍

輸入?yún)?shù)中有alpha,可改變其顯著水準,其預設值為0.05,或為95%之信任水準。

輸出項中,means 即為各分組項之平均值,sem為該分組項之標準差,counts為該分組下之觀察值數(shù)目,而name則為該分組之名稱。

範例一:



 ?? x = 1???? 2???? 3???? 4???? 5???? 6???? 7???? 8???? 9??? 10
 group= 1???? 1???? 1???? 1???? 1???? 2???? 2???? 2???? 2???? 2;

>> grpstats(x,group)

ans =

3
8

上述結(jié)果為分兩組的平均,前五項為一組,後五項為一組。結(jié)果第一組平均為3,第二組為8。
組別間,其項數(shù)並不一定要相同,例如:

範例二:



x =
? 1???? 2???? 3???? 4???? 5???? 6???? 7???? 8???? 9

>> group=[1 1 1 1 2 2 2 2 2]
group =

? 1???? 1???? 1???? 1???? 2???? 2???? 2???? 2???? 2

>> [m,s,c]=grpstats(x,group)

m =

2.5000
7.0000


s =

0.6455
0.7071

c =

? 4
? 5

其輸出之第一項為平均值,第二項為標準差,第三項為各組之項數(shù)。故即使各組之樣本數(shù)不同也可以得到對應組之統(tǒng)計資料。

範例三:


設有200個觀測值分成四小組,每一觀測值分成五項,其平均範圍由1-5。為製造這樣的數(shù)據(jù),下面之例子實際上應用了許多特定的函數(shù):
  • unidrnd(4,100,1) 平均製造一個100X1的陣列,其中之數(shù)值分配為1:4的整數(shù)範圍,以每項分別以1,2,3,4隨機出現(xiàn)。
  • normrnd(true_mean,1) 常態(tài)分配之亂數(shù)函數(shù),其平均值為true_mean,其標準差為1。
  • true_mean((ones(100,1),:) 利用原來設定之(ones(100,1),:)陣列,重覆100次。

執(zhí)行此程式後,由於n為細胞陣列,故全改為字串才能同時顯現(xiàn)其結(jié)果,其結(jié)果如下:


group = unidrnd(4,100,1);%create a 100X1 matrix in random [ 1,2,3,4 ]
true_mean = 1:5;
true_mean = true_mean(ones(100,1),:); %100X5 matrix
x = normrnd(true_mean,1); %randomize
[m, s, c,n] = grpstats(x,group);

[n num2cell(m)]
ans =
'1'??? [0.9584]??? [1.8200]??? [2.8412]??? [4.1669]??? [5.0220]
'2'??? [0.8972]??? [1.8393]??? [2.9423]??? [4.0044]??? [4.9437]
'3'??? [0.9768]??? [2.1093]??? [3.1565]??? [3.9860]??? [5.0585]
'4'??? [1.1164]??? [2.2249]??? [2.8920]??? [4.1323]??? [5.3251]


範例四:


利用matlab所附的carsmall.mat示範檔案,其中參數(shù)項目包括重量(Weight)、年份(Model_Year)等資料,利用該項資料求其年份下之平均車重、預測值、年份名稱及各年份下之數(shù)量。最後並利用errorbar繪出其範圍。

% cargroup.m
load carsmall
[Weight Model_Year]'
[means,p,year,count] = grpstats(Weight,Model_Year,...
{'mean','predci','gname','numel'})
n = length(means);
errorbar((1:n)',means,p(:,2)-means)
set(gca,'xtick',1:n,'xticklabel',year)
title('95% prediction intervals for mean weight by year')


先將上述程式存為cargroup.m檔案,執(zhí)行後應有許多資料,其中僅選擇本題所需要者。其過程如下:

>> cargroup

ans =

Columns 1 through 7

???? 3504??????? 3693??????? 3436??????? 3433??????? 3449??????? 4341??????? 4354
?????? 70????????? 70????????? 70????????? 70????????? 70????????? 70????????? 70

Columns 8 through 14

???? 4312??????? 4425??????? 3850??????? 3090??????? 4142??????? 4034??????? 4166
?????? 70????????? 70????????? 70????????? 70????????? 70????????? 70????????? 70

= == == == == == == == == == == ==
Columns 92 through 98

???? 2835??????? 2665??????? 2370??????? 2950??????? 2790??????? 2130??????? 2295
?????? 82????????? 82????????? 82????????? 82????????? 82????????? 82????????? 82

Columns 99 through 100

???? 2625??????? 2720
?????? 82????????? 82


means =

1.0e+003 *

3.4413
3.0787
2.4535


p =

1.0e+003 *

1.7770??? 5.1056
1.3832??? 4.7742
1.7184??? 3.1887


year =

'70'
'76'
'82'


count =

35
34
31


?

11.2.9表列函數(shù)tabulate



TABLE = tabulate(x)

這個指令可將一向量X之觀測值製成一表格,其第一行為X向量中之相同數(shù)值,第二行為該數(shù)值出現(xiàn)之次數(shù),最後一行為該值出現(xiàn)之百分比。若X值為含有文字串之陣列或細胞陣列,則第一行為陣列內(nèi)之獨一的名稱,其餘兩行則相同。下面為利用rand之隨機函數(shù)取100個值作比較。

tabulate(round(rand(100,1)*6))?
? Value??? Count?? Percent
????? 0??????? 8????? 8.00%
????? 1?????? 16???? 16.00%
????? 2?????? 22???? 22.00%
????? 3?????? 21???? 21.00%
????? 4?????? 10???? 10.00%
????? 5?????? 16???? 16.00%
????? 6??????? 7????? 7.00%?


tabulate(fix(rand(100,1)*6))?
? Value??? Count?? Percent
????? 0?????? 14???? 14.00%
????? 1?????? 14???? 14.00%
????? 2?????? 21???? 21.00%
????? 3?????? 15???? 15.00%
????? 4?????? 22???? 22.00%
????? 5?????? 14???? 14.00%?

由上面之結(jié)果比較可以看出:利用四捨五入法之round函數(shù),其在兩端之機率顯然少很多。
利用前節(jié)之汽車carsmall.mat資料,亦可以tabulate指令作簡單之統(tǒng)計:

load carsmall
>> tabulate(Origin)
??? Value??? Count?? Percent
????? USA?????? 69???? 69.00%
?? France??????? 4????? 4.00%
??? Japan?????? 15???? 15.00%
? Germany??????? 9????? 9.00%
?? Sweden??????? 2????? 2.00%
??? Italy??????? 1????? 1.00%

看起來美國地區(qū)還是美國車多。

上述tablulate指令之左邊若不給予參數(shù),則會自動產(chǎn)生上述之格式,分成三行,即名稱、數(shù)量及百分比。若結(jié)果給予一個參數(shù)時,其內(nèi)容會轉(zhuǎn)為細胞陣列。因此必要時,須利用cell2mat函數(shù)轉(zhuǎn)換為數(shù)值陣列。以上述資料為例,下面的型式會有不同的結(jié)果:

>> a=tabulate(Origin)

a =

??? 'USA'??????? [69]??? [69]
??? 'France'???? [ 4]??? [ 4]
??? 'Japan'????? [15]??? [15]
??? 'Germany'??? [ 9]??? [ 9]
??? 'Sweden'???? [ 2]??? [ 2]
??? 'Italy'????? [ 1]??? [ 1]

顯然其內(nèi)容均是細胞陣列,若要作成餅圖,則必須稍微變換:

>> pie(cell2mat(a(:,2)),a(:,1))


?

11.2.10百分值prctile


在一般大量樣本之情況下,可以利用百分值去確定樣本之合理對應值,由此百分比與對應值之關(guān)係可以瞭解資料之外形、位置以及擴散度。其指令格式如下:

Y = prctile(X, p)

此指令計算X之樣本值中一個大於p%部份之對應值位置,此值並不一定是原有之觀測值,只求其比例位置。輸入?yún)?shù) p 必須落在[0 100]間,可為常數(shù)或向量。若 p = 50% 時,則Y值應對應X之中間值(median)。X之資料可為向量或矩陣,而 p 則可能為一向量或其中之常數(shù)。下表說明可能之之四種狀況:

X????? p?????????????? Y= prctile(X, p)
====?? ========?? ===========================================
向量?? 數(shù)值常數(shù)?? 數(shù)值等於X中第p級百分值
矩陣?? 數(shù)值常數(shù)?? 列向量含X中每一行內(nèi)第p級百分值,其中
????????????????????? y(i)=prctile(X(:,j),p)
向量?? 向量?????? Y向量與p同,其第i項為X中第p(i)級百分值
????????????????????? y(i)=prctile(X,p(i))
矩陣?? 向量?????? Y矩陣中,其第(i,j)項為X中第j行第p(i)級百分值
????????????????????? y(i,j)=prctile(X(:,j),p(i))
====?? ========?? ===========================================

執(zhí)行例:



>> x=randn(100,1);
>> Y=prctile(x,[10:10:90])
Y =

-1.2258 -1.0402 -0.6829 -0.4320 -0.2074? 0.0287? 0.3215? 0.6855? 1.2649

%Y與p之向量相同,與X為列或行向量無關(guān)。

%X2若為矩陣,則p先與X之行向量作百段分值。

X2 =
? 0???? 1???? 2???? 3???? 4
? 1???? 2???? 3???? 4???? 5
? 2???? 3???? 4???? 5???? 6

>>Y2=prctile(X2,50)
Y2 =
? 1???? 2???? 3???? 4???? 5

>>X2=[0:4;1:5; 2:6],Y2=prctile(X2,[20 30 50])

X2 =
? 0???? 1???? 2???? 3???? 4
? 1???? 2???? 3???? 4???? 5
? 2???? 3???? 4???? 5???? 6

>>Y2 =
0.1000??? 1.1000??? 2.1000??? 3.1000??? 4.1000
0.4000??? 1.4000??? 2.4000??? 3.4000??? 4.4000
1.0000??? 2.0000??? 3.0000??? 4.0000??? 5.0000

X2若為矩陣,p為向量時,每個p(i)先與X之行向量作百段分值,並Y依序完成行向量之值。

?

?

11.2.11細分值quantile



細分值與百分值之意義類似,但其區(qū)間為小數(shù),介於[0 1]之間,以配合累積密度函數(shù)之使用,其指令格式如下:

Y = quantile(X, p)
Y = quantile(X, p, dim)

其輸出值Y為X觀測值中傳回值,p為數(shù)值或累積機率值之向量。當X為向量時,Y之大小與p相同,而Y(i)則是第p(i)之細分對應值。當X為矩陣,則Y之第i列為第p(i)對X之行向量之細分值。

其細分方向亦可利用dim設定,但Y在dim指定的方向長度應與p之長度相同。

細分值常應用於累加機率,故其範圍為[0 1]。若X為一具N元素之向量,則QUANTILE依下列方式運算:
  • 1)X值經(jīng)過排序,並細分為(0.5/N), (1.5/n), ..., ((N-0.5)/N) 細分段
  • 2) 線性內(nèi)插法用於計算 (0.5/N) 與 ((N-0.5)/N)間機率之細分段。
  • 3) X中之最大值與最小值作為機率外圍之細段。


X=1:10;y=quantile(X,.50)?
y =
?? 5.5000

y = quantile(X,[.025 .25 .50 .75 .975])
y =
?? 1.0000??? 3.0000??? 5.5000??? 8.0000?? 10.0000

?

堆疊矩陣之使用,前面也曾述及,其相關(guān)語法如下:


B = repmat(A,m,n)
B = repmat(A,[m n])
B = repmat(A,[m n p...])

這是一個處理大矩陣且內(nèi)容有重複時使用之。其功能是以A之內(nèi)容堆疊在一(M x N)的矩陣B中。B矩陣之大小由MXN及A矩陣之內(nèi)容決定。例如:

>>B=repmat( [1 2;3 4],2,3)

B =
?? 1???? 2???? 1???? 2???? 1???? 2
?? 3???? 4???? 3???? 4???? 3???? 4
?? 1???? 2???? 1???? 2???? 1???? 2
?? 3???? 4???? 3???? 4???? 3???? 4

其結(jié)果變?yōu)?X6。A也可以置放文字串,如:

>>C=repmat(' Long live the king!', 2,2)
C =
Long live the king! Long live the king!
Long live the king! Long live the king!

也可置放特定常數(shù):

>> D=repmat(NaN,2,5)

D =
NaN?? NaN?? NaN?? NaN?? NaN
NaN?? NaN?? NaN?? NaN?? NaN

?

?

11.2.13中央慣性矩MOMENT


統(tǒng)計分析之技巧中,通常以觀察值之平均值為參考標,經(jīng)處理後再進行判定。MOMENT的指令就是採用相差距離之次方作為計算之基準,如:

mk = E(x-u)k

式中之m值變成以不同k值產(chǎn)生之中央慣性矩,也是各種差異度之期望值。k值為與均值差之次方,又稱階數(shù)(order),其值必須為正。利用觀察資料依各值對中央均值之值距作階數(shù)之乘方而累加之和,再除以樣本數(shù)。當k=1時稱為第一慣性矩,所得結(jié)果應為零,亦即為均值之位置;k=2時即為第二慣性矩,即為上述討論之變方,只是此時之除數(shù)為n。其指令型式如下:

m = moment(X,order)
舉例:

>> X=randn(4,6)
m1=moment(X,1)
m2=moment(X,2)
m3=moment(X,3)
m4=moment(X,4)

X =
-1.7258??? 0.1387??? 1.1508?? -0.3735?? -1.5731?? -0.2433
0.8132?? -0.8595?? -0.6080?? -0.8320??? 2.0157??? 0.1733
1.4419?? -0.7523??? 0.8062??? 0.2869?? -0.0720??? 0.9232
0.6723??? 1.2296??? 0.2171?? -1.8189??? 2.6289?? -0.1786
m1 =
0???? 0???? 0???? 0???? 0???? 0
m2 =
1.4524??? 0.7053??? 0.4445??? 0.5872??? 2.8011??? 0.2149
m3 =
-1.6611??? 0.3293?? -0.1237?? -0.1293?? -1.1069??? 0.0795
m4 =
4.6600??? 0.8526??? 0.3402??? 0.6391?? 11.1516??? 0.0919

上述之結(jié)果可知,m1均為零,m2與m4因為偶數(shù)次方故均為正數(shù);m3則有正負。

?

?

11.2.14歪度SKEWNESS


歪度為第三中央慣性矩除以標準差之三次方,由此可以測出分佈的偏向,或稱為費雪(Fisher)歪度公式。其計算方式如下:

? s = E(x-μ)33

其中μ為x之平均值,σ為標準差,E表示括號內(nèi)之期望值。y即為所謂之歪度。其指令之型式如下:

s = skewness(X)
s = skewness(X,flag)

輸入?yún)?shù)flag為校正偏差用,flag=1時不作校正,為其預設值;flag=0時則有校正。通常少量之樣本代表一群體時,所得歪度會產(chǎn)生偏差,故視樣本大小需進行校正。因此,修正時,必須執(zhí)行skewness(X,0)指令。

若x為向量,則s 會得到樣本歪度,或為單一值。當X為矩陣時,s歪度將以列向量表示。若其結(jié)果s等於零,表示對稱分佈;s>0則為正偏(即高峰偏左);s<0為負偏(高峰在右)。進行判斷時,若其絕對值小於1.96,則屬常態(tài)分配,大於1.96則為非常態(tài)。例如:

>> X=randn(4,5)
s=skewness(X)

X =
2.0941??? 1.6820?? -0.1586?? -0.5266??? 0.2486
0.0802??? 0.5936??? 0.8709?? -0.6855??? 0.1025
-0.9373??? 0.7902?? -0.1948?? -0.2684?? -0.0410
0.6357??? 0.1053??? 0.0755?? -1.1883?? -2.2476
s =
0.2794??? 0.4979??? 0.9682?? -0.4978?? -1.1201

由歪度結(jié)果得知,這此數(shù)值雖屬常態(tài)分佈,但也有中央偏向問題。第一、二、三行之高峰偏左,第四及第五行則高峰偏右。

?

?

11.2.15峰度KURTOSIS


峰度是採用第四中央慣性矩除以標準差的四次方而得。其公式如下:

? k =E(x-μ)44

此處k值稱為kurtosis,其指令格式如下:

k = kurtosis(X)
k = kurtosis(X,flag)

其中flag之功能與skewness相同,為校正偏差用,flag=1時為其預設值,不作校正;flag=0作校正。kurtosis(X)之意義在於表示峰度之趨勢分佈。在常態(tài)分配中,峰度值為3。大於3時其峰度將高於常態(tài)峰度;小於3時則低於常態(tài)峰度。以此可以作為山峰高度之判斷。
例如:

X=randn(100,5);k=kurtosis(X)

k =
?? 3.1023??? 2.6613??? 2.7877??? 3.5796??? 2.7580

這種常態(tài)分佈狀態(tài)下,當樣本增多時,峰度值將趨近於3。上述五行之資料當中,第一及四行之山峰較高,其餘之峰頂則較低。(峰度值3即表示標準的山形的意識,也比較容易記)

?

在大部份之機率分佈計算中,有些指令尚提供參數(shù)及信任水準的計算方法。為獲得相關(guān)指令,可以使用help stats進行查詢。例如:


>> help stats
? Statistics Toolbox
? Version 5.0 (R14) 05-May-2004

?

Distributions.


?? Parameter estimation.
??? betafit???? - Beta parameter estimation.
??? binofit???? - Binomial parameter estimation.
??? dfittool??? - Distribution fitting tool.
??? evfit?????? - Extreme value parameter estimation.
??? expfit????? - Exponential parameter estimation.
??? gamfit????? - Gamma parameter estimation.
??? lognfit???? - Lognormal parameter estimation.
??? mle???????? - Maximum likelihood estimation (MLE).
??? mlecov????? - Asymptotic covariance matrix of MLE.
??? nbinfit???? - Negative binomial parameter estimation.
??? normfit???? - Normal parameter estimation.
??? poissfit??? - Poisson parameter estimation.
??? raylfit???? - Rayleigh parameter estimation.
??? unifit????? - Uniform parameter estimation.
??? wblfit????? - Weibull parameter estimation.

??

Probability density functions (pdf).


??? betapdf???? - Beta density.
??? binopdf???? - Binomial density.
??? chi2pdf???? - Chi square density.
??? evpdf?????? - Extreme value density.
??? exppdf????? - Exponential density.
??? fpdf??????? - F density.
??? gampdf????? - Gamma density.
??? geopdf????? - Geometric density.
??? hygepdf???? - Hypergeometric density.
??? lognpdf???? - Lognormal density.
??? mvnpdf????? - Multivariate normal density.
??? nbinpdf???? - Negative binomial density.
??? ncfpdf????? - Noncentral F density.
??? nctpdf????? - Noncentral t density.
??? ncx2pdf???? - Noncentral Chi-square density.
??? normpdf???? - Normal (Gaussian) density.
??? pdf???????? - Density function for a specified distribution.
??? poisspdf??? - Poisson density.
??? raylpdf???? - Rayleigh density.
??? tpdf??????? - T density.
??? unidpdf???? - Discrete uniform density.
??? unifpdf???? - Uniform density.
??? wblpdf????? - Weibull density.

??

Cumulative Distribution functions (cdf).


??? betacdf???? - Beta cdf.
??? binocdf???? - Binomial cdf.
??? cdf???????? - Specified cumulative distribution function.
??? chi2cdf???? - Chi square cdf.
??? ecdf??????? - Empirical cdf (Kaplan-Meier estimate).
??? evcdf?????? - Extreme value cumulative distribution function.
??? expcdf????? - Exponential cdf.
??? fcdf??????? - F cdf.
??? gamcdf????? - Gamma cdf.
??? geocdf????? - Geometric cdf.
??? hygecdf???? - Hypergeometric cdf.
??? logncdf???? - Lognormal cdf.
??? nbincdf???? - Negative binomial cdf.
??? ncfcdf????? - Noncentral F cdf.
??? nctcdf????? - Noncentral t cdf.
??? ncx2cdf???? - Noncentral Chi-square cdf.
??? normcdf???? - Normal (Gaussian) cdf.
??? poisscdf??? - Poisson cdf.
??? raylcdf???? - Rayleigh cdf.
??? tcdf??????? - T cdf.
??? unidcdf???? - Discrete uniform cdf.
??? unifcdf???? - Uniform cdf.
??? wblcdf????? - Weibull cdf.

??

Critical Values of Distribution functions.


??? betainv???? - Beta inverse cumulative distribution function.
??? binoinv???? - Binomial inverse cumulative distribution function.
??? chi2inv???? - Chi square inverse cumulative distribution function.
??? evinv?????? - Extreme value inverse cumulative distribution function.
??? expinv????? - Exponential inverse cumulative distribution function.
??? finv??????? - F inverse cumulative distribution function.
??? gaminv????? - Gamma inverse cumulative distribution function.
??? geoinv????? - Geometric inverse cumulative distribution function.
??? hygeinv???? - Hypergeometric inverse cumulative distribution function.
??? icdf??????? - Specified inverse cdf.
??? logninv???? - Lognormal inverse cumulative distribution function.
??? nbininv???? - Negative binomial inverse distribution function.
??? ncfinv????? - Noncentral F inverse cumulative distribution function.
??? nctinv????? - Noncentral t inverse cumulative distribution function.
??? ncx2inv???? - Noncentral Chi-square inverse distribution function.
??? norminv???? - Normal (Gaussian) inverse cumulative distribution function.
??? poissinv??? - Poisson inverse cumulative distribution function.
??? raylinv???? - Rayleigh inverse cumulative distribution function.
??? tinv??????? - T inverse cumulative distribution function.
??? unidinv???? - Discrete uniform inverse cumulative distribution function.
??? unifinv???? - Uniform inverse cumulative distribution function.
??? wblinv????? - Weibull inverse cumulative distribution function.

??

Random Number Generators.


??? betarnd???? - Beta random numbers.
??? binornd???? - Binomial random numbers.
??? chi2rnd???? - Chi square random numbers.
??? evrnd?????? - Extreme value random numbers.
??? exprnd????? - Exponential random numbers.
??? frnd??????? - F random numbers.
??? gamrnd????? - Gamma random numbers.
??? geornd????? - Geometric random numbers.
??? hygernd???? - Hypergeometric random numbers.
??? iwishrnd??? - Inverse Wishart random matrix.
??? lognrnd???? - Lognormal random numbers.
??? mvnrnd????? - Multivariate normal random numbers.
??? mvtrnd????? - Multivariate t random numbers.
??? nbinrnd???? - Negative binomial random numbers.
??? ncfrnd????? - Noncentral F random numbers.
??? nctrnd????? - Noncentral t random numbers.
??? ncx2rnd???? - Noncentral Chi-square random numbers.
??? normrnd???? - Normal (Gaussian) random numbers.
??? poissrnd??? - Poisson random numbers.
??? randg?????? - Gamma random numbers (unit scale).
??? random????? - Random numbers from specified distribution.
??? randsample? - Random sample from finite population.
??? raylrnd???? - Rayleigh random numbers.
??? trnd??????? - T random numbers.
??? unidrnd???? - Discrete uniform random numbers.
??? unifrnd???? - Uniform random numbers.
??? wblrnd????? - Weibull random numbers.
??? wishrnd???? - Wishart random matrix.

??

Statistics.


??? betastat??? - Beta mean and variance.
??? binostat??? - Binomial mean and variance.
??? chi2stat??? - Chi square mean and variance.
??? evstat????? - Extreme value mean and variance.
??? expstat???? - Exponential mean and variance.
??? fstat?????? - F mean and variance.
??? gamstat???? - Gamma mean and variance.
??? geostat???? - Geometric mean and variance.
??? hygestat??? - Hypergeometric mean and variance.
??? lognstat??? - Lognormal mean and variance.
??? nbinstat??? - Negative binomial mean and variance.
??? ncfstat???? - Noncentral F mean and variance.
??? nctstat???? - Noncentral t mean and variance.
??? ncx2stat??? - Noncentral Chi-square mean and variance.
??? normstat??? - Normal (Gaussian) mean and variance.
??? poisstat??? - Poisson mean and variance.
??? raylstat??? - Rayleigh mean and variance.
??? tstat?????? - T mean and variance.
??? unidstat??? - Discrete uniform mean and variance.
??? unifstat??? - Uniform mean and variance.
??? wblstat???? - Weibull mean and variance.

?? Likelihood functions.
??? betalike??? - Negative beta log-likelihood.
??? evlike????? - Negative extreme value log-likelihood.
??? explike???? - Negative exponential log-likelihood.
??? gamlike???? - Negative gamma log-likelihood.
??? lognlike??? - Negative lognormal log-likelihood.
??? nbinlike??? - Negative likelihood for negative binomial distribution.
??? normlike??? - Negative normal likelihood.
??? wbllike???? - Negative Weibull log-likelihood.

?

Descriptive Statistics.


??? bootstrp??? - Bootstrap statistics for any function.
??? corr??????? - Linear or rank correlation coefficient.
??? corrcoef??? - Linear correlation coefficient with confidence intervals.
??? cov???????? - Covariance.
??? crosstab??? - Cross tabulation.
??? geomean???? - Geometric mean.
??? grpstats??? - Summary statistics by group.
??? harmmean??? - Harmonic mean.
??? iqr???????? - Interquartile range.
??? kurtosis??? - Kurtosis.
??? mad???????? - Median Absolute Deviation.
??? mean??????? - Sample average (in MATLAB toolbox).
??? median????? - 50th percentile of a sample.
??? moment????? - Moments of a sample.
??? nanmax????? - Maximum ignoring NaNs.
??? nanmean???? - Mean ignoring NaNs.
??? nanmedian?? - Median ignoring NaNs.
??? nanmin????? - Minimum ignoring NaNs.
??? nanstd????? - Standard deviation ignoring NaNs.
??? nansum????? - Sum ignoring NaNs.
??? nanvar????? - Variance ignoring NaNs.
??? prctile???? - Percentiles.
??? quantile??? - Quantiles.
??? range?????? - Range.
??? skewness??? - Skewness.
??? std???????? - Standard deviation (in MATLAB toolbox).
??? tabulate??? - Frequency table.
??? trimmean??? - Trimmed mean.
??? var???????? - Variance (in MATLAB toolbox).

?

Linear Models.


??? addedvarplot - Created added-variable plot for stepwise regression.
??? anova1????? - One-way analysis of variance.
??? anova2????? - Two-way analysis of variance.
??? anovan????? - n-way analysis of variance.
??? aoctool???? - Interactive tool for analysis of covariance.
??? dummyvar??? - Dummy-variable coding.
??? friedman??? - Friedman's test (nonparametric two-way anova).
??? glmfit????? - Generalized linear model fitting.
??? glmval????? - Evaluate fitted values for generalized linear model.
??? kruskalwallis - Kruskal-Wallis test (nonparametric one-way anova).
??? leverage??? - Regression diagnostic.
??? lscov?????? - Least-squares estimates with known covariance matrix.
??? lsqnonneg?? - Non-negative least-squares.
??? manova1???? - One-way multivariate analysis of variance.
??? manovacluster - Draw clusters of group means for manova1.
??? multcompare - Multiple comparisons of means and other estimates.
??? polyconf??? - Polynomial evaluation and confidence interval estimation.
??? polyfit???? - Least-squares polynomial fitting.
??? polyval???? - Predicted values for polynomial functions.
??? rcoplot???? - Residuals case order plot.
??? regress???? - Multivariate linear regression.
??? regstats??? - Regression diagnostics.
??? ridge?????? - Ridge regression.
??? robustfit?? - Robust regression model fitting.
??? rstool????? - Multidimensional response surface visualization (RSM).
??? stepwise??? - Interactive tool for stepwise regression.
??? stepwisefit - Non-interactive stepwise regression.
??? x2fx??????? - Factor settings matrix (x) to design matrix (fx).

?

Nonlinear Models.


??? nlinfit???? - Nonlinear least-squares data fitting.
??? nlintool??? - Interactive graphical tool for prediction in nonlinear models.
??? nlpredci??? - Confidence intervals for prediction.
??? nlparci???? - Confidence intervals for parameters.

? Design of Experiments (DOE).
??? bbdesign??? - Box-Behnken design.
??? candexch??? - D-optimal design (row exchange algorithm for candidate set).
??? candgen???? - Candidates set for D-optimal design generation.
??? ccdesign??? - Central composite design.
??? cordexch??? - D-optimal design (coordinate exchange algorithm).
??? daugment??? - Augment D-optimal design.
??? dcovary???? - D-optimal design with fixed covariates.
??? ff2n??????? - Two-level full-factorial design.
??? fracfact??? - Two-level fractional factorial design.
??? fullfact??? - Mixed-level full-factorial design.
??? hadamard??? - Hadamard matrices (orthogonal arrays).
??? lhsdesign?? - Latin hypercube sampling design.
??? lhsnorm???? - Latin hypercube multivariate normal sample.
??? rowexch???? - D-optimal design (row exchange algorithm).

?

Statistical Process Control (SPC).


??? capable???? - Capability indices.
??? capaplot??? - Capability plot.
??? ewmaplot??? - Exponentially weighted moving average plot.
??? histfit???? - Histogram with superimposed normal density.
??? normspec??? - Plot normal density between specification limits.
??? schart????? - S chart for monitoring variability.
??? xbarplot??? - Xbar chart for monitoring the mean.

?

Multivariate Statistics.


?? Cluster Analysis.
??? cophenet??? - Cophenetic coefficient.
??? cluster???? - Construct clusters from LINKAGE output.
??? clusterdata - Construct clusters from data.
??? dendrogram? - Generate dendrogram plot.
??? inconsistent - Inconsistent values of a cluster tree.
??? kmeans????? - k-means clustering.
??? linkage???? - Hierarchical cluster information.
??? pdist?????? - Pairwise distance between observations.
??? silhouette? - Silhouette plot of clustered data.
??? squareform? - Square matrix formatted distance.

??

Dimension Reduction Techniques.


??? factoran??? - Factor analysis.
??? pcacov????? - Principal components from covariance matrix.
??? pcares????? - Residuals from principal components.
??? princomp??? - Principal components analysis from raw data.
??? rotatefactors - Rotation of FA or PCA loadings.

??

Plotting.


??? andrewsplot - Andrews plot for multivariate data.
??? biplot????? - Biplot of variable/factor coefficients and scores.
??? glyphplot?? - Plot stars or Chernoff faces for multivariate data.
??? gplotmatrix - Matrix of scatter plots grouped by a common variable.
??? parallelcoords - Parallel coordinates plot for multivariate data.

??

Other Multivariate Methods.


??? barttest??? - Bartlett's test for dimensionality.
??? canoncorr?? - Cannonical correlation analysis.
??? cmdscale??? - Classical multidimensional scaling.
??? classify??? - Linear discriminant analysis.
??? mahal?????? - Mahalanobis distance.
??? manova1???? - One-way multivariate analysis of variance.
??? mdscale???? - Metric and non-metric multidimensional scaling.
??? procrustes? - Procrustes analysis.

?

Decision Tree Techniques.


??? treedisp??? - Display decision tree.
??? treefit???? - Fit data using a classification or regression tree.
??? treeprune?? - Prune decision tree or creating optimal pruning sequence.
??? treetest??? - Estimate error for decision tree.
??? treeval???? - Compute fitted values using decision tree.

?

Hypothesis Tests.


??? ranksum???? - Wilcoxon rank sum test (independent samples).
??? signrank??? - Wilcoxon sign rank test (paired samples).
??? signtest??? - Sign test (paired samples).
??? ztest?????? - Z test.
??? ttest?????? - One sample t test.
??? ttest2????? - Two sample t test.

?

Distribution Testing.


??? jbtest????? - Jarque-Bera test of normality
??? kstest????? - Kolmogorov-Smirnov test for one sample
??? kstest2???? - Kolmogorov-Smirnov test for two samples
??? lillietest? - Lilliefors test of normality

?

Nonparametric Functions.


??? friedman??? - Friedman's test (nonparametric two-way anova).
??? kruskalwallis - Kruskal-Wallis test (nonparametric one-way anova).
??? ksdensity?? - Kernel smoothing density estimation.
??? ranksum???? - Wilcoxon rank sum test (independent samples).
??? signrank??? - Wilcoxon sign rank test (paired samples).
??? signtest??? - Sign test (paired samples).

?

Hidden Markov Models.


??? hmmdecode?? - Calculate HMM posterior state probabilities.
??? hmmestimate - Estimate HMM parameters given state information.
??? hmmgenerate - Generate random sequence for HMM.
??? hmmtrain??? - Calculate maximum likelihood estimates for HMM parameters.
??? hmmviterbi? - Calculate most probable state path for HMM sequence.

?

Statistical Plotting.


??? andrewsplot - Andrews plot for multivariate data.
??? biplot????? - Biplot of variable/factor coefficients and scores.
??? boxplot???? - Boxplots of a data matrix (one per column).
??? cdfplot???? - Plot of empirical cumulative distribution function.
??? ecdfhist??? - Histogram calculated from empirical cdf.
??? fsurfht???? - Interactive contour plot of a function.
??? gline?????? - Point, drag and click line drawing on figures.
??? glyphplot?? - Plot stars or Chernoff faces for multivariate data.
??? gname?????? - Interactive point labeling in x-y plots.
??? gplotmatrix - Matrix of scatter plots grouped by a common variable.
??? gscatter??? - Scatter plot of two variables grouped by a third.
??? hist??????? - Histogram (in MATLAB toolbox).
??? hist3?????? - Three-dimensional histogram of bivariate data.
??? lsline????? - Add least-square fit line to scatter plot.
??? normplot??? - Normal probability plot.
??? parallelcoords - Parallel coordinates plot for multivariate data.
??? probplot??? - Probability plot.
??? qqplot????? - Quantile-Quantile plot.
??? refcurve??? - Reference polynomial curve.
??? refline???? - Reference line.
??? surfht????? - Interactive contour plot of a data grid.
??? wblplot???? - Weibull probability plot.

?

Statistics Demos.


??? aoctool???? - Interactive tool for analysis of covariance.
??? disttool??? - GUI tool for exploring probability distribution functions.
??? polytool??? - Interactive graph for prediction of fitted polynomials.
??? randtool??? - GUI tool for generating random numbers.
??? rsmdemo???? - Reaction simulation (DOE, RSM, nonlinear curve fitting).
??? robustdemo? - Interactive tool to compare robust and least squares fits.

?

File Based I/O.


??? tblread???? - Read in data in tabular format.
??? tblwrite??? - Write out data in tabular format to file.
??? tdfread???? - Read in text and numeric data from tab-delimitted file.
??? caseread??? - Read in case names.
??? casewrite?? - Write out case names to file.

?

Utility Functions.


??? combnk????? - Enumeration of all combinations of n objects k at a time.
??? grp2idx???? - Convert grouping variable to indices and array of names.
??? hougen????? - Prediction function for Hougen model (nonlinear example).
??? statget???? - Get STATS options parameter value.
??? statset???? - Set STATS options parameter value.
??? tiedrank??? - Compute ranks of sample, adjusting for ties.
??? zscore????? - Normalize matrix columns to mean 0, variance 1.

?

事件發(fā)生的機率雖然可以估計,但是必須由過去發(fā)生的事實去統(tǒng)計,並以此預測未來將發(fā)生的事件。這種運作的方式在模擬的作業(yè)上有實際的困難,故僅能依靠亂數(shù)產(chǎn)生器模擬事件之發(fā)生。亂數(shù)的特性是發(fā)生時間或次數(shù)均屬隨機性,每次發(fā)生的機率均依循其特定的模式,常用者為均勻分佈的亂數(shù)型態(tài),每一事件出現(xiàn)之機率是相等的。利用亂數(shù)產(chǎn)生之數(shù)值可以模擬機械之故障率、顧客來訪的時機、旅客搭機的安排,甚或穀物乾燥中心作業(yè)的動線規(guī)劃等。

在MATLAB中常用的亂數(shù)產(chǎn)生指令為rand,這是個依均勻機率分配出現(xiàn)之原則產(chǎn)生一個或一組在[0, 1]間之亂數(shù),每次呼叫其值均不一樣。這些數(shù)值雖有其範圍,但產(chǎn)生之後可以作適當調(diào)整,以符合實際所需。所產(chǎn)生之數(shù)值可為單一數(shù)字,亦可為特定之矩陣,其大小由其後之參數(shù)界定之。若僅輸入一個參數(shù)(n),其所產(chǎn)生之結(jié)果為一nxn之方矩陣;若為(m, n)則得 m x n之矩陣。

例如:


>>rand?
ans =
??? 0.9501?

>>rand(2,4)???
ans =
??? 0.2311??? 0.4860??? 0.7621??? 0.0185
??? 0.6068??? 0.8913??? 0.4565??? 0.8214?


除均勻機率分配產(chǎn)生者外,亦可使用randn之指令,其產(chǎn)生之亂數(shù)則係依常態(tài)分態(tài)的型式。例如:

>> randn(5,6)

ans =

??? 0.8956??? 0.5689?? -0.2340??? 0.6232??? 0.2379??? 0.3899
??? 0.7310?? -0.2556??? 0.1184??? 0.7990?? -1.0078??? 0.0880
??? 0.5779?? -0.3775??? 0.3148??? 0.9409?? -0.7420?? -0.6355
??? 0.0403?? -0.2959??? 1.4435?? -0.9921??? 1.0823?? -0.5596
??? 0.6771?? -1.4751?? -0.3510??? 0.2120?? -0.1315??? 0.4437

利用常態(tài)分佈之指令randn可以產(chǎn)生正負值,其值也可能大於一。其他分配所需之亂數(shù)產(chǎn)生器可以參考如下之各項指令。注意每項指令之後可能需要相關(guān)之參數(shù),讀者可以利用help之輔助指令獲得所需之資訊。

亂數(shù)產(chǎn)生器(Random Number Generators)


依據(jù)實際之用途,可以有不同之亂數(shù)指令可以應用,

??? 均勻分佈亂數(shù)rand
??? 常態(tài)分佈亂數(shù)randn
 ? 常態(tài)亂數(shù)(Normal (Gaussian))-normrnd.
??? 貝他亂數(shù)(Beta)-betarnd.
??? 二項式亂數(shù)(Binomial)- binornd.
??? X's亂數(shù)(Chi square)-chi2rnd.
??? 極值亂數(shù)(Extreme value)-evrnd.
??? 指數(shù)亂數(shù)(Exponential)-exprnd.
??? F亂數(shù)(F)-frnd.
??? 伽瑪亂數(shù)(Gamma)-gamrnd.
??? 伽瑪亂數(shù)( Gamma (unit scale))-randg.
??? 幾何亂數(shù)(Geometric)-geornd.
??? 超幾何亂數(shù)(Hypergeometric)-hygernd.
??? 反威夏亂數(shù)(Inverse Wishart)-iwishrnd.
??? 常態(tài)對數(shù)亂數(shù)(Lognormal)-lognrnd.
??? 多變異亂數(shù)(Multivariate normal)-mvnrnd.
??? 多變T亂數(shù)( Multivariate t)-mvtrnd.
??? 負二項式亂數(shù)(Negative binomial)-nbinrnd .
??? 非中央F亂數(shù)(Noncentral F)-ncfrnd .
??? 非中央亂數(shù)(Noncentral)-nctrnd.
??? 非中央X's亂數(shù)(Noncentral Chi-square)-ncx2rnd.
??? 波義松亂數(shù)(Poisson)-poissrnd.
??? 定群樣本亂數(shù)(finite population)-randsample .
??? 魏利亂數(shù)( Rayleigh )-raylrnd .
??? T亂數(shù)(T)- trnd .
??? 離散均勻亂數(shù)( Discrete uniform)-unidrnd .
??? 均勻亂數(shù)(Uniform)-unifrnd .
??? 魏伯亂數(shù)(Weibull)-wblrnd? .
??? 威夏亂數(shù)矩陣( Wishart)-wishrnd .

上述諸多指令,可以實際統(tǒng)計上之需要實際選取應用。例如:以亂數(shù)的理念,產(chǎn)生一群數(shù)字作排例時,可用randperm(n)指令,以產(chǎn)生不同之排列組合:

>> randperm(7)
ans =
???? 1???? 3???? 7???? 6???? 2???? 5???? 4

>> randperm(7)
ans =
???? 7???? 6???? 3???? 5???? 4???? 2???? 1

>> randperm(7)
ans =
???? 4???? 1???? 7???? 2???? 6???? 3???? 5

其結(jié)果每次顯示不同組合,均為整數(shù)。此指令因此可以模擬擲骰子,若改用randperm(6),則成為每趟擲六次,顯現(xiàn)排列數(shù)據(jù)。這個指令與rand指令不同,若使用rand則必須先化成整數(shù),然後利用程式進行模擬擲骰之動作,其過程比較複雜。

MATLAB亦提供一示範指令,randtool。此指令可以示範各種型態(tài)之分配機率,利用下拉式清單可以選擇所需的分配型式,並直接設定所需之範圍與變異量。其結(jié)果以對應之質(zhì)方圖顯示。

>> randtool


由於不同之分佈,有不同的指令,不容易記清楚。不過,MATLAB提供一個通用的指令,其型式如下:

random('name', a1,a2,m,n)

這個指令則可在name處置放各種指令名稱,其後可依需要再加上對應之參數(shù),其中a1、a2為該函數(shù)所需之參數(shù),而最後之[m,n]則為最後之矩陣大小,變成一個相當有彈性的指令,例如:

random('normal',0,1,3,3)?

ans =
?? -0.1867??? 2.1832??? 1.0668
??? 0.7258?? -0.1364??? 0.0593
?? -0.5883??? 0.1139?? -0.0956?

其中之'name'名稱可以為 'beta', 'binomial', 'chi-square' 或 'chi2', 'extreme value' 或 'ev', 'exponential', 'f', 'gamma', 'geometric', 'hypergeometric' 或 'hyge', 'lognormal', 'negative binomial' 或 'nbin', 'noncentral f' 或 'ncf', 'noncentral t' 或 'nct', 'noncentral chi-square' 或 'ncx2', 'normal', 'poisson', 'rayleigh', 't', 'discrete uniform' 或 'unid', 'uniform', 'weibull' 或 'wbl'等等。部份字母符合亦可,且大小寫不拘。如此可以免去背記上述指令之苦。

?

常態(tài)分配之亂數(shù)中,不像均勻分佈的型態(tài),基本上它會集中在某區(qū)域,故最常發(fā)生的事件會集中在平均值附近。在均勻分佈型態(tài)應用上,常需設定上下界,讓其出現(xiàn)限定在特定範圍內(nèi);常態(tài)分配則無明顯的上下界限,且由於集中存在平均值附近,故與均值間之距離會有正負差。其表示方法如下:


y = μ+σx


其中μ為平均值,而σ為其標準差,亦即常態(tài)分配之亂數(shù)值。在MATLAB中以randn之指令函數(shù)產(chǎn)生標準差值σ。此指令係以標準差為1,而平均值為0作成亂數(shù)結(jié)果,故若有某一群組分佈平均為10,標準差為5,且樣本為200點,則產(chǎn)生後必須乘以所需之標準差5再加上平均值10,其做法如下表示:

y=10+5*randn(1, 200)
>>hist(10+5*randn(1,200))???


這是一個交談式的指令,可以直接變換參數(shù)並得其結(jié)果:


>> disttool

這是一個展示各種分配及機率分佈之指令。沒有任何參數(shù),但可以在圖中更改所要的參數(shù)。開始時可選分佈種類,如常態(tài)、貝他、二項式、指數(shù)、二極值、F、伽瑪、幾何、常態(tài)對數(shù)等等。選定其一後立即會顯示。其後選擇所需函數(shù)名稱。僅有機率分配(pdf)、累積分配函數(shù)(cdf)等二項。上兩項選定後,即可輸入X值,求得對應之機率,或選定機率求得對應之X值。亦可利用滑鼠直接在圖上點出位置,顯示其對應之X值與機率。此外,在其左下方有參數(shù)設定值,可以利用滑尺為之。

單向變異分析(ANOVA1)旨在尋找不同類別下之資料是否具有相同均值,亦即決定不同類別下,各量測值是否具有差異特性。在線性模式下,單向變異分析是最簡單的狀況。任何量測值可以矩陣表示如下:


yij = αjεij

其中,yij中每行j代表不同類別,並具有類別內(nèi)之平均值 。在工具箱中存有hogg.mat資料,下載後可以作為執(zhí)行anova1函數(shù)之用,其過程如下:

>> load hogg
>> hogg

hogg =
24??? 14??? 11???? 7??? 19
15???? 7???? 9???? 7??? 24
21??? 12???? 7???? 4??? 19
27??? 17??? 13???? 7??? 15
33??? 14??? 12??? 12??? 10
23??? 16??? 18??? 18??? 20

ANOVA指令之型式如下:

P = ANOVA1(X,GROUP,DISPLAYOPT)

此指令會將X組合之樣本矩陣中之各行視為比較之組別,以此決定組別間之平均值是否會相等。GROUP為分組向量,其長度應與X之行數(shù)相同,其內(nèi)容即為各組之名稱。指令左邊為結(jié)果之機率,其值愈小表示相等之機率愈小,亦即顯著性之差異愈大。GROUP向量可以容許使用字串或細胞陣列,但必須以行表示。若不使用組別名稱則可以空格。

DISPLAYOPT為顯示開關(guān),可用 'on' (預設值)或 'off'決定是否顯示結(jié)果圖。若要文字輸出ANOVA表,則可在左邊加一參數(shù)如ANOVATAB:

[P,ANOVATAB, stats] = ANOVA1(...)

參數(shù)ANOVATAB為細胞陣列,stats則為一些統(tǒng)計參數(shù)。

例:

>> [p,tbl,stats]=anova1(hogg)
p =
0.00011971
tbl =
'Source'???? 'SS'??????? 'df'??? 'MS'??????? 'F'???????? 'Prob>F'
'Columns'??? [?? 803]??? [ 4]??? [200.75]??? [9.0076]??? [0.00011971]
'Error'????? [557.17]??? [25]??? [22.287]????????? []????????????? []
'Total'????? [1360.2]??? [29]????????? []????????? []????????????? []
stats =
gnames: [5x1 char]
????? n: [6 6 6 6 6]
source: 'anova1'
? means: [23.833 13.333 11.667 9.1667 17.833]
???? df: 25
????? s: 4.7209


在本例中,hogg之資料實際上各行代表不同處理之細菌繁殖情形。利用ANOVA變異分析表之結(jié)果亦可做 F-test,以證明不同處理是否仍然具有相同的結(jié)果。此例所得之p值約為 0.0001,這是一個相當小的值。換句話說,此種結(jié)果強烈顯示不同的處理其結(jié)果之差異是顯著的。若就或然率考慮,每10,000次實驗當中只有一次有結(jié)果相同的機會。對於研究者而言,這是一個很大的鼓舞,因為至少由統(tǒng)計分析結(jié)果可證明:不同處理比較下,應有很大的差異性的。當然,上述p值要正確,其基本假設是各項變異應獨立,且屬常態(tài)分配,其變異值也是固定。這些差異性由圖也可以看出端倪。



上述之指令,若要加上各行之名稱,則可另以如下之指令為之:

>> [p,tbl,stats]=anova1(hogg,{'A1','A2','A3','A4','A5'}')

其執(zhí)行結(jié)果與前述無異,但在最後之組別圖中會出現(xiàn)分組之組名。

?

?

11.7.1Multiple Comparisons

有些時候光比較眾處理之間有無差異仍然不足,有時必須知道那一對有顯著的差異。當然此時也可以藉助一連串的t-test逐對比較,但這種方式仍有缺點。在做t-test時,通常會先計算t的過程中與一特定值比較。此值之選定通常是認為各處理之均值超過此值很小(例如5%)。當均值有顯著差異時,統(tǒng)計的數(shù)據(jù)超過該值的機率反而變大,容易造成誤判。利用MULTCOMPARE指令可以解決此問題,並可執(zhí)行平均值之多重比較。其語法如下:

?? COMPARISON = MULTCOMPARE(stats)

參數(shù)stats為統(tǒng)計資料,可以直接得自 anova1, anova2, anovan, aoctool等函數(shù)之輸出值,其輸出參數(shù)COMPARISON則是一個矩陣,每一列中具有五行,第 1-2行是被比較之兩樣本編號,第 3-5 則分別為其差異之最低值、估計值及最高值。

>>COMPARISON = MULTCOMPARE(stats)
COMPARISON =
?????????? 1??????????? 2?????? 2.4953???????? 10.5?????? 18.505
?????????? 1??????????? 3?????? 4.1619?????? 12.167?????? 20.171
?????????? 1??????????? 4?????? 6.6619?????? 14.667?????? 22.671
?????????? 1??????????? 5????? -2.0047??????????? 6?????? 14.005
?????????? 2??????????? 3????? -6.3381?????? 1.6667?????? 9.6714
?????????? 2??????????? 4????? -3.8381?????? 4.1667?????? 12.171
?????????? 2??????????? 5????? -12.505???????? -4.5?????? 3.5047
?????????? 3??????????? 4????? -5.5047????????? 2.5?????? 10.505
?????????? 3??????????? 5????? -14.171????? -6.1667?????? 1.8381
?????????? 4??????????? 5????? -16.671????? -8.6667???? -0.66193

單變異與雙變異分析上不同的地方是後者具有兩組類別,其所定義的特性不相同。例如某汽車公司有兩家工廠,每家工廠均生產(chǎn)相同的三車種。此時生產(chǎn)之汽車之里程數(shù)作合理的比較時,其差異可能有兩項:其一是甲乙兩工廠間之差異,其二是車型間之差異。 此時工廠別與車型別兩者均會影響汽車之里程數(shù)。其間差異可能來自工廠間之製造,也可能是因車種設計或規(guī)格上的問題,後者之問題可能與工廠無關(guān)。 此外,某工廠也可能對生產(chǎn)某一車種有獨到的製造能力(可能因為有較佳的生產(chǎn)線吧!),卻製造其他車種則與另一工廠不相上下,這種效應稱為加減性,或兩項類別間之交互作用所產(chǎn)生之影響。

雙變異分析(ANOVA)是線性模式之特殊狀況。就汽車之里程數(shù)可以表示如下:


  yijk=μ+α.ji.ijijk

其中,yijk為里程數(shù)之觀察資料(列指標為i,行指標為j,重複標為k)。μ整個矩陣之平均值。α.j為由於車型j與μ差異之均值,βi.則是工廠別i與μ差異之均值,γij為交互作用項,其在行向之和或列向之和為零。εijk 為整個隨機之變異量。

由上述汽車之製造例可知,其觀察之里程數(shù)變成一種矩陣的型式,具有j行與i列的組別項,此時由行與列構(gòu)成之組別或稱為處理(Treatment),對應於行列方向之交叉點即為細胞(Cell),每個細胞位置必須重複置放樣本觀察數(shù),或稱為重複數(shù)( repetition)。若以矩陣表示,此重複數(shù)必須置於k方向。

ANOVA2為分析雙變異數(shù)之指令。其格式如下:

  ANOVA2(X, REPS, DISPLAYOPT)

輸入?yún)?shù)X為觀察之資料矩陣,其行與列均需二項以上,以作為比較之基準。在各行的資料代表一類別;各列者則為另一類別。若每一對應細胞有多個觀察資料,則由REPS指定。若REPS=3,代表每個細胞位置有三筆資料,這些資料依列之類別依序置於列中。因此若REPS=3,則在列中每三筆資料屬其中一組,下三筆屬第二組,如此類推,其總列數(shù)因而應為三(REPS)的倍數(shù)。另外參數(shù)DISPLAYOPT之定義與前述指令之用法相同。

輸出參數(shù)則表示如下:

[P, TABLE, stats] = ANOVA2(...)

其中 P為p值之向量,而 TABLE則為細胞陣列,包括ANOVA表的內(nèi)容。stats則為分析後之統(tǒng)計資料,可以作為其他指令如MULTCOMPARE函數(shù)之輸入。

若有非平衡變異資料則可採用ANOVAN指令進行分析,詳細情形可以參考手冊或輔助資料。
茲以工具箱中存有mileage.mat之車型生產(chǎn)資料為例,可以作為本節(jié)之分析示範:

>>load mileage
>>mileage
mileage =
??????? 33.3???????? 34.5???????? 37.4
??????? 33.4???????? 34.8???????? 36.8
??????? 32.9???????? 33.8???????? 37.6
?????? 32.6???????? 33.4???????? 36.6
?????? 32.5???????? 33.7???????? 37
 ??????? 33???????? 33.9???????? 36.7

此為三種車型(行)及二家工廠製造之汽車每加侖里程資料,其reps=3,故第1-3列代表第一家工廠製造之車輛,第4-6列代表第二家工廠。將相關(guān)參數(shù)代入anova2指令,可得結(jié)果為:

>>[p,tbl,stats] = anova2(mileage,3)
p =
2.4278e-010??? 0.0039095????? 0.84106

tbl =
'Source'???????? 'SS'??????? 'df'??? 'MS'???????? 'F'????????? 'Prob>F'
'Columns'???? [53.351]??? [ 2]??? [ 26.676]??? [ 234.22]??? [2.4278e-010]
'Rows'??????? [ 1.445]??? [ 1]??? [? 1.445]??? [ 12.688]??? [? 0.0039095]
'Interaction' [? 0.04]??? [ 2]??? [?? 0.02]??? [0.17561]??? [??? 0.84106]
'Error'?????? [1.3667]??? [12]??? [0.11389]?????????? []?????????????? []
'Total'?????? [56.203]??? [17]?????????? []?????????? []?????????????? []
stats =
? source: 'anova2'
sigmasq: 0.11389
colmeans: [32.95 34.017 37.017]
??? coln: 6
rowmeans: [34.944 34.378]
??? rown: 9
?? inter: 1
??? pval: 0.84106
????? df: 12




除文字及統(tǒng)計數(shù)據(jù)外,尚有一個標準的 ANOVA表,其中第一行為行,其平方和(SS)、自由度(df)及均方值 (MS=SS/df), F 檢定及或然率 p-值等。 F檢定可以檢驗車型間、工廠間及車型X工廠之交互作用(經(jīng)過調(diào)整後之增加效應)產(chǎn)生之里程數(shù)是否相同。由於 p-值在車型間之效應幾乎近於零(至小數(shù)四位),表示車型對里程數(shù)之變異甚大,有顯著之影響。經(jīng)由 F檢定結(jié)果,因車型而發(fā)生同平均值之情況,其機率為萬分之一。若使用MULTCOMPARE指令比較:

>>COMPARISON = MULTCOMPARE(stats)

Note:? Your model includes an interaction term.? A test of main
effects can be difficult to interpret when the model includes
interactions.
COMPARISON =
??????? 1??????????? 2????? -1.5865????? -1.0667???? -0.54686
??????? 1??????????? 3????? -4.5865????? -4.0667????? -3.5469
??????? 2??????????? 3????? -3.5198?????????? -3????? -2.4802


利用multcompare函數(shù)比較結(jié)果,三車型之差異有明顯之不同。而 p-值對工廠間之機率值為 0.0039,這也是一個非常顯著之差異。顯然,某一工廠製造的汽車,其里程數(shù)是比另一廠為高。由其 p-值得知,只有千分之四的機率兩工廠製造的汽車之里程數(shù)才會相同。就工廠X車型間交互作用之影響而言,則不顯著。其 p-值僅 0.8411,亦即結(jié)果中,可能百分八四機率會出現(xiàn)無交互作用之影響。當然,此處顯示之 p-值若要正確,基本上整個樣本之分佈必須獨立、常態(tài)分配並有固定的變異常數(shù)。

?

?

11.7.3多變異分析ANOVAN

在變異分析的過程中,有些資料並非刻意形成,而是因調(diào)查或分析產(chǎn)生的特定資訊,無法事先規(guī)劃作出等量分組或先經(jīng)實驗分組。例如,有一些汽車銷售資訊,在銷售單中可能有不同車種、型號、排氣量,甚至出產(chǎn)國別等。若就其分類有些數(shù)量並不一定對等,或處於非平衡狀態(tài)。多變異分析可以同時處理平衡或非平衡資料,其功能與anova1、anova2略有不同。此模式之架構(gòu)如下:


其中有相乘之兩項或三項者,即為其間之交互作用(interaction)產(chǎn)生之效應。在所需之X資料向量中,應有其對應之類別資訊。依其類別資訊可以區(qū)分為若干層次,此時可將要設為層次之行歸納在一個細胞陣列 group中。其內(nèi)可有N層之名單,以此指認X資料之歸屬。在group細胞中之名單可以是一向量、文字陣列或文字細胞陣列,以大括符集合起來,惟其個數(shù)需與X中之項目相符。其相關(guān)語法如下:

p = anovan(x,group)
p = anovan(x,group,'Param1',val1,'Param2',val2,...)
[p,table] = anovan(...)
[p,table,stats] = anovan(...)
[p,table,stats,terms] = anovan(...)

其參數(shù)定義大體上與anova1、anova2等指令相同,每個 GROUP中之變數(shù)則必須具有與X等量之元素。在分析後之變異數(shù)分析表中,將列出GROUP中之變數(shù)名稱。另外增加之參數(shù)'Paramx',valx,則可說明如下:







'Paramx'valx
'sstype'1, 2, 或 3(預設值)平方和之型式
'varnames'組別名稱預設為 'X1', 'X2', 'X3', ..., 'XN' ,但可利用此參數(shù)定義新名稱
'display'顯示ANOVA表,'on' | 'off'
'random'以一向量指定群組參數(shù)是否為隨機。
'alpha'信任水準(預設值為0.05或95%)
'model'決定採用何種函數(shù),如 'linear'其p值僅計算N主效應(預設值); 'interaction'除主效應外尚計算兩者間之交互效應; 'full'計算所有主效應及交互效應。若為整數(shù)k(k<=N),則表示計算至k層級。如k=3,表示主效應+兩者間及三者間之交互效應。 若為零壹值矩陣之定義,如[1 0 0] A主效應、[0 1 0] B主效應、[0 0 1 ] C主效應、[1 1 0] AB交互效應、[1 0 1] AC交互效應、[0 1 1] BC交互效應、[1 1 1] ABC交互效應等。


在輸出項當中,TERMS矩陣可作為下一次執(zhí)行ANOVA時,其 MODEL輸入之引數(shù)。若MODEL參數(shù)不採用隨機型態(tài),其變異分析表T 中之欄位將為TERMS 、SS、DF、奇偶數(shù)、均方、F檢測及P值等。若具隨機功能時,則會顯示TERMS之型態(tài)(即固定或隨機)、期望均方、F值之除數(shù)均方、F值之除數(shù)自由度、除數(shù)定義、變異項估計值、低限值及高限值等。

輸出參數(shù)STATS 的結(jié)構(gòu)包括使用MULTCOMPARE函數(shù)之參數(shù)及如下項:

coeffs????? 估計係數(shù)
coeffnames? 各係數(shù)之名稱
vars??????? 組別參數(shù)值矩陣
resid????? 適配模式下之餘數(shù)。

具隨機效應時,則另有下列欄位:

ems???????? 期望均方值。
denom?????? 分母定義。
rtnames???? 隨機項名稱。
varest????? 變異項估計值(每一隨機項一個)。
varci?????? 各變異項之信任區(qū)間

實際上anova2與anovan也可以通用,只是其安排的語法略有不同而已。anova2是將重複的觀察值置於列向,用reps來計算其每組之重複次數(shù),而anovan則必須另外建相關(guān)的行作為指示其階數(shù)之用。例如下面的例子,假設屬anova2之資料,且是二重複,即reps=2:

X = [25 17 21; 28 18 65; 45 5 60; 42 10 95]

X =
25??? 17??? 21
28??? 18??? 65
45???? 5??? 60
42??? 10??? 95
p=anova2(X,2)
p =
0.0175??? 0.1930??? 0.2321




上述資料格式若要改為適用anovan,則必須建立對應之行向量。首先建立一個對應三欄的位置m1;其次再建立對應重複的群組m2:

m1=repmat(1:3,4,1)
m2=[ones(2,3);2*ones(2,3)]

m1 =
1???? 2???? 3
1???? 2???? 3
1???? 2???? 3
1???? 2???? 3
m2 =
1???? 1???? 1
1???? 1???? 1
2???? 2???? 2
2???? 2???? 2

其次,再將X,m1,m2以行向?qū)饋?#xff0c;如:

X=X(:);m1=m1(:);m2=m2(:);
[X m1 m2]
ans =
25???? 1???? 1
28???? 1???? 1
45???? 1???? 2
42???? 1???? 2
17???? 2???? 1
18???? 2???? 1
5???? 2???? 2
10???? 2???? 2
21???? 3???? 1
65???? 3???? 1
60???? 3???? 2
95???? 3???? 2

此時就可以利用anovan執(zhí)行分析了。

anovan(X, {m1 m2},2)
ans =
0.0175
0.1930
0.2321

其結(jié)果與上述anova1相同。看起來使用naovan好像比較複雜,那是因為轉(zhuǎn)換為anova1之型式複雜的綠故,若準備資料時就朝naovan之需要處理,則會簡單許多。指令中之參數(shù)2即為'model'=2的意思,表示分析結(jié)果包括交互作用之影響。若不考慮交互作用,則可以免去,其結(jié)果將完全與anova1相同。

下面是一個具三方變異分析之例子,設y為觀察之資料,其餘三組g1、g2、g3為對應之組合項,各組內(nèi)僅有兩種內(nèi)容,g1為[1 2]、g2['hi' 'lo']、g3['june' 'may']。由三組內(nèi)容可知,其項目可為數(shù)值,可為文字,但文字必須使用文字字串。

y = [52.7 57.5 45.9 44.5 53.0 57.0 45.9 44.0]';
g1 = [1 2 1 2 1 2 1 2];
g2 = {'hi';'hi';'lo';'lo';'hi';'hi';'lo';'lo'};
g3 = {'may'; 'may'; 'may'; 'may'; 'june'; 'june'; 'june'; 'june'};

利用上述資料進行三方變異分析如下:

>>p = anovan(y, {g1 g2 g3})
p =
0.4174
0.0028
0.9140





輸出項P值包括三個主項所產(chǎn)生之效應,由其變異分析表可知分為x1、x2、x3對應g1、g2、g3等三組。由於沒有第三項參數(shù),故僅列示主效率部份。若需其交互效應則必須加入2之參數(shù)在最後一項,或令'model' = 2。三個P值則對應個別零假設(Null Hypothesis)Ho1、Ho2、Ho3之機率值。基本上P值趨近零,對於零假設的成分愈小。例如,Ho2之對應P值0.0028已足夠小到可以認定該組合項下之某一樣本平均值無法與另一組有顯著之差別。在進行此項研斷之前,通常必須先選定P值之顯著水準門檻值,常用的的有0.005或 0.001,依問題的特性而定。

上述的結(jié)果並未考慮交互效應,若要增加此項功能,必須在參數(shù)方面設定,即令'mode'=2,或直接使用2即可。例如:

>>[p,table]=anovan(y, {g1 g2 g3},'model',2)

p =
0.0347
0.0048
0.2578
0.0158
0.1444
0.5000
table =
'Source'? 'Sum Sq.' 'd.f.' 'Singular?' 'Mean Sq.'??? 'F'?????? 'Prob>F'
'X1'???? [? 3.7812] [?? 1] [?? 0]??? [? 3.7812]?? 336.1111]??? [0.0347]
'X2'???? [199.0013] [?? 1] [?? 0]??? [199.0013]? [1.7689e+004] [0.0048]
'X3'???? [? 0.0612] [?? 1] [?? 0]??? [? 0.0612]? [???? 5.4444] [0.2578]
'X1*X2'? [ 18.3012] [?? 1] [?? 0]??? [ 18.3012]? [1.6268e+003] [0.0158]
'X1*X3'? [? 0.2112] [?? 1] [?? 0]??? [? 0.2112]? [??? 18.7778] [0.1444]
'X2*X3'? [? 0.0113] [?? 1] [?? 0]??? [? 0.0113]? [????????? 1] [0.5000]
'Error'? [? 0.0113] [?? 1] [?? 0]??? [? 0.0113]????????? []??? []
'Total'? [221.3788] [?? 7] [?? 0]??????????? []????????? []??? []



結(jié)果會有三項值,其內(nèi)容大致相同。P值之首三項為主效應,其餘三項為交互效應。在名稱上均預設為x1、x2、x3,此名稱亦可自行設定:

>>p= anovan(y,{g1 g2 g3},'varnames',{'time' 'limit' 'month'},'model',2);




>> aoctool
Note: 6 observations with missing values have been removed.

AOCTOOL 適合單向變異分析(ANOCVA)及繪製變量圖。其執(zhí)行格式如下:

AOCTOOL(X,Y,G,ALPHA)

上式之輸入?yún)?shù):X,Y為對應之資料,其中X為已知值,Y為回應值。G則是作為分組之變數(shù)。而ALPHA則是信任水準,其預設值為0.05,亦即其信任水準為95%。

執(zhí)行這個程式後,將出現(xiàn)三個圖表,其中一個為交談曲線圖,可以調(diào)整參數(shù)使其產(chǎn)生預測值變化,其二為為變異分析表(ANOVA),另一為參數(shù)值預測。

除上述輸入?yún)?shù)外, 各參數(shù)尚可指定名稱,其圖表亦可開關(guān)。其指令型式如下:

AOCTOOL(X,Y,G,ALPHA,XNAME,YNAME,GNAME,DISPLAYOPT,MODEL)

其中XNAME,YNAME,GNAME分別為X軸、Y軸及組別名稱。 DISPLAYOPT 則控制圖形是否顯示('on'為預設值,或否'off'), MODEL 則控制最初適配時之選項,包括:

1??? 'same mean' (共同均值)
2??? 'separate means' (分開均值)
3??? 'same line' (共同預測線)
4??? 'parallel lines' (平行預測線)
5??? 'separate lines' (分離預測線)

此外,此一指令亦容許有輸出,其型式如下:

[H,ATAB,CTAB,STATS] = AOCTOOL(...)

其中 H代表圖中各線之握把, ATAB 與 CTAB則為ANOVA表及係數(shù)表之內(nèi)容,而 STATS則為多重均值比較之結(jié)果,亦即為 MULTCOMPARE 函數(shù)執(zhí)行的結(jié)果。

下面的例子中為利用load carsmall下載最近十餘年(1970-1982)的汽車相關(guān)資料。其中包括車重(Weight)、加侖里程數(shù)(MPG)及年份(Model_Year)。其資料內(nèi)容如下,請?zhí)貏e注意在加侖里程之資料中有六點之資料從缺。

>> Weight'
3504???? 3693???? 3436???? 3433???? 3449???? 4341???? 4354
4312???? 4425???? 3850???? 3090???? 4142???? 4034???? 4166
3850???? 3563???? 3609???? 3353???? 3761???? 3086???? 2372
2833???? 2774???? 2587???? 2130???? 1835???? 2672???? 2430
2375???? 2234???? 2648???? 4615???? 4376???? 4382???? 4732
2464???? 2220???? 2572???? 2255???? 2202???? 4215???? 4190
3962???? 4215???? 3233???? 3353???? 3012???? 3085???? 2035
2164???? 1937???? 1795???? 3651???? 3574???? 3645???? 3193
1825???? 1990???? 2155???? 2565???? 3150???? 3940???? 3270
2930???? 3820???? 4380???? 4055???? 3870???? 3755???? 2605
2640???? 2395???? 2575???? 2525???? 2735???? 2865???? 3035
1980???? 2025???? 1970???? 2125???? 2125???? 2160???? 2205
2245???? 1965???? 1965???? 1995???? 2945???? 3015???? 2585
2835???? 2665???? 2370???? 2950???? 2790???? 2130???? 2295
2625???? 2720


>> MPG'
18.0000?? 15.0000?? 18.0000?? 16.0000?? 17.0000?? 15.0000?? 14.0000?? 14.0000
14.0000?? 15.0000?????? NaN?????? NaN?????? NaN?????? NaN?????? NaN?? 15.0000
14.0000?????? NaN?? 15.0000?? 14.0000?? 24.0000?? 22.0000?? 18.0000?? 21.0000
27.0000?? 26.0000?? 25.0000?? 24.0000?? 25.0000?? 26.0000?? 21.0000?? 10.0000
10.0000?? 11.0000??? 9.0000?? 28.0000?? 25.0000?? 25.0000?? 26.0000?? 27.0000
17.5000?? 16.0000?? 15.5000?? 14.5000?? 22.0000?? 22.0000?? 24.0000?? 22.5000
29.0000?? 24.5000?? 29.0000?? 33.0000?? 20.0000?? 18.0000?? 18.5000?? 17.5000
29.5000?? 32.0000?? 28.0000?? 26.5000?? 20.0000?? 13.0000?? 19.0000?? 19.0000
16.5000?? 16.5000?? 13.0000?? 13.0000?? 13.0000?? 28.0000?? 27.0000?? 34.0000
31.0000?? 29.0000?? 27.0000?? 24.0000?? 23.0000?? 36.0000?? 37.0000?? 31.0000
38.0000?? 36.0000?? 36.0000?? 36.0000?? 34.0000?? 38.0000?? 32.0000?? 38.0000
25.0000?? 38.0000?? 26.0000?? 22.0000?? 32.0000?? 36.0000?? 27.0000?? 27.0000
44.0000?? 32.0000?? 28.0000?? 31.0000


>> Model_Year'
70??? 70??? 70??? 70??? 70??? 70??? 70??? 70??? 70??? 70??? 70??? 70??? 70??? 70
70??? 70??? 70??? 70??? 70??? 70??? 70??? 70??? 70??? 70??? 70??? 70??? 70??? 70
70??? 70??? 70??? 70??? 70??? 70??? 70??? 76??? 76??? 76??? 76??? 76??? 76??? 76
76??? 76??? 76??? 76??? 76??? 76??? 76??? 76??? 76??? 76??? 76??? 76??? 76??? 76
76??? 76??? 76??? 76??? 76??? 76

再分享一下我老師大神的人工智能教程吧。零基礎!通俗易懂!風趣幽默!還帶黃段子!希望你也加入到我們?nèi)斯ぶ悄艿年犖橹衼?#xff01;https://blog.csdn.net/jiangjunshow

總結(jié)

以上是生活随笔為你收集整理的MATLAB统计与回归的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯,歡迎將生活随笔推薦給好友。

国产亚洲精品久久久久久网站 | 日韩簧片在线观看 | 国产精品免费看久久久8精臀av | 久久综合在线 | 91女人18片女毛片60分钟 | 99热精品国产一区二区在线观看 | 91资源在线免费观看 | 亚洲黄色大片 | 久草亚洲视频 | 国产精品久久久久久久久久ktv | 色999视频 | 免费看污污视频的网站 | 9在线观看免费高清完整 | 天天做天天看 | 成年人电影免费看 | 91视频在线免费下载 | 免费观看性生活大片3 | 9免费视频 | 在线视频观看成人 | 91污污 | 欧美性色黄大片在线观看 | 成人91在线观看 | 国产午夜一区 | 国产精品久久久久免费观看 | 在线观看完整版免费 | 香蕉在线影院 | 不卡国产在线 | av丝袜制服 | 中文字幕在线日本 | 岛国精品一区二区 | 亚洲欧洲av在线 | 在线网址你懂得 | 国产91在线播放 | 午夜久久久久久久 | 国产精品久久久久毛片大屁完整版 | 国产伦精品一区二区三区四区视频 | 久久免费看a级毛毛片 | 天天插天天狠天天透 | 亚洲国产成人精品久久 | 中文字幕在线播放视频 | 国产在线2020 | 91人人网 | 88av网站 | 欧美一进一出抽搐大尺度视频 | 欧美一二区在线 | 蜜臀av.com| 国产999精品 | 亚洲精品久久久久中文字幕m男 | 国产大尺度视频 | 色视频网站在线观看一=区 a视频免费在线观看 | 久久免费精品 | 久久免费美女视频 | av蜜桃在线 | 高清免费av在线 | 一二三区在线 | 亚洲日本一区二区在线 | 91中文字幕网 | 亚洲精品乱码久久久久久久久久 | av看片在线| 久久亚洲专区 | 美女天天操 | 久久99精品国产 | 久久精品79国产精品 | 婷婷电影在线观看 | 91新人在线观看 | 亚洲精品午夜久久久久久久 | 成人动漫精品一区二区 | av丝袜在线 | 久久天天躁夜夜躁狠狠躁2022 | 免费成人黄色片 | 九九热有精品 | 色99视频 | 免费在线观看av网址 | 亚洲理论在线观看电影 | 美女视频黄是免费的 | 久久久精品 | 亚洲免费精彩视频 | 狠狠的干狠狠的操 | 色吊丝在线永久观看最新版本 | 精品国产一区二区三区免费 | 最近中文字幕高清字幕免费mv | 在线免费观看视频a | 九九热精品视频在线观看 | 精品毛片久久久久久 | 欧美一区,二区 | 伊色综合久久之综合久久 | 国产手机视频 | 一区二区三区动漫 | 亚洲三区在线 | 欧美日韩一级久久久久久免费看 | 在线免费观看黄网站 | 国产成人高清在线 | 精品久久片 | 97国产精品亚洲精品 | av综合在线观看 | 国产精品入口66mio女同 | 中字幕视频在线永久在线观看免费 | 国产色视频网站 | 欧美性生活一级片 | 丁香视频 | 中文字幕永久免费 | 日韩中文字幕免费 | 中文字幕一区二区三区视频 | 久草在线欧美 | 日韩午夜网站 | 五月情婷婷 | 狠狠88综合久久久久综合网 | 欧美 亚洲 另类 激情 另类 | 成人一级片免费看 | 偷拍久久久 | 久久久久久久久久免费 | 91在线看片 | 四虎国产精品免费观看视频优播 | 亚洲激情中文 | 亚洲精品小视频 | 91麻豆精品国产91久久久无限制版 | 国内久久看 | 在线观看一区二区精品 | 九九色在线 | 久久成年人网站 | 日韩欧美视频免费在线观看 | 日韩中文字 | 91在线公开视频 | 天天干,天天插 | 久久久精品网站 | 日韩视频1| 蜜臀av性久久久久蜜臀aⅴ四虎 | 一区三区视频 | 国产精品久久久久久久久久三级 | 欧美一区二区日韩一区二区 | 日本精品免费看 | 日韩精品高清不卡 | 久久五月婷婷丁香 | 一本色道久久精品 | 三级黄色网络 | 久久av电影| 亚洲禁18久人片 | 久久手机在线视频 | 久久男人中文字幕资源站 | 亚洲综合精品在线 | 欧美视频在线观看免费网址 | 91在线看黄 | 99视频在线精品国自产拍免费观看 | 国产美女精品视频 | 亚洲免费av片 | 91免费观看视频网站 | 99久久婷婷国产综合亚洲 | 免费黄a| 91精品1区| 超碰在线人人爱 | 日韩字幕在线观看 | 9797在线看片亚洲精品 | 99久国产 | 国产中文字幕视频在线观看 | 亚欧日韩av | 又长又大又黑又粗欧美 | 亚洲h色精品 | 亚洲理论片在线观看 | 国产精品11| 国产精品久久久久永久免费看 | 欧美日韩大片在线观看 | 国产精品二区在线 | 久日精品 | 国产精品岛国久久久久久久久红粉 | 国产一级一片免费播放放a 一区二区三区国产欧美 | 开心婷婷色 | 国产精品久久久久久一区二区 | 国内精品亚洲 | 91在线一区| 国产精品第10页 | 午夜 久久 tv| 国产在线无 | 日本三级不卡视频 | av片一区| 久久精品黄 | 日韩视频专区 | 免费观看性生活大片 | 不卡av免费在线观看 | 亚洲国产精品小视频 | 毛片网在线观看 | 日韩电影在线观看一区 | 亚洲在线视频播放 | 国内毛片毛片 | 一区二区高清在线 | 国内精品久久久久久久久 | 黄色毛片视频免费观看中文 | 狠狠干网址 | 久久久久免费看 | 99精品久久久久久久久久综合 | av中文字幕免费在线观看 | 精品电影一区二区 | 亚洲黄色一级电影 | 看黄色91| 91香蕉视频 mp4 | av免费在线播放 | 国产精品乱码一区二区视频 | japanese黑人亚洲人4k | 亚洲影视九九影院在线观看 | 国产精品大片免费观看 | 久久老司机精品视频 | 九九热精品视频在线观看 | 夜夜躁天天躁很躁波 | 日本aaa在线观看 | 久久美女视频 | 在线看国产视频 | 人人爽人人爽人人爽人人爽 | 91夫妻视频 | 日本爱爱免费视频 | 国产精品第7页 | 国产小视频免费在线观看 | 狠狠色狠狠色综合日日小说 | 四虎在线免费 | 国产一区二区三区免费观看视频 | 久草视频国产 | 狠狠狠狠狠狠天天爱 | 五月天婷婷丁香花 | 日韩中文在线电影 | 国产亚洲视频中文字幕视频 | 伊人午夜 | 日韩av电影中文字幕在线观看 | 日韩,精品电影 | 五月婷婷香蕉 | 精品久久久久久久久久久久 | 天堂在线一区二区三区 | 国产一区网址 | 精品久久久精品 | av电影 一区二区 | 99免在线观看免费视频高清 | 欧美精品国产综合久久 | 制服丝袜一区二区 | 国产亚洲精品久久网站 | 久久兔费看a级 | 青草视频在线免费 | 中文字幕精品一区二区三区电影 | 91色一区二区三区 | 天堂av高清 | 久久免费视频播放 | 手机看片久久 | 天天激情综合网 | 免费观看一级一片 | 欧美一区二区三区在线播放 | 亚洲精品视频网址 | 欧美天天射 | 西西4444www大胆无视频 | 亚州性色 | 日本精品视频免费 | 久久99影院 | 国产一级电影免费观看 | 最近中文字幕大全中文字幕免费 | 日韩在线观看一区二区三区 | 免费a视频在线观看 | 亚欧日韩av | 91在线观看黄 | 99久久精品视频免费 | 国产精品久久久久久久免费 | 欧美最新另类人妖 | 成人精品在线 | 欧美久久久久久久久久久久久 | 天天草视频 | 亚洲午夜精品福利 | 国产在线精品一区二区三区 | 亚洲精品综合一区二区 | 久久国产剧场电影 | 精品日韩在线一区 | 精品亚洲va在线va天堂资源站 | 欧美日韩成人一区 | 亚洲片在线资源 | 久久人人插 | 日韩在线播放av | 成人一级片免费看 | 成年人在线免费看视频 | 国产精品久久网站 | 精品国产伦一区二区三区观看说明 | 干狠狠| 亚洲天堂网在线视频 | 亚洲免费一级 | 国产精品99久久久久人中文网介绍 | 免费黄色av. | 日韩在线电影一区 | 日韩欧美在线观看一区 | 亚洲精品中文在线观看 | 国产精品午夜免费福利视频 | 特级黄录像视频 | 91免费黄视频 | 亚洲欧美日本A∨在线观看 青青河边草观看完整版高清 | 综合在线观看色 | 免费精品在线 | 在线免费视频 你懂得 | 国产免费观看久久黄 | 在线观看成人毛片 | 国产不卡高清 | 色五月色开心色婷婷色丁香 | 日本久久不卡视频 | 国产精品视频免费看 | 国产精品av免费观看 | 三级av在线 | 欧美日韩成人一区 | 色网站在线免费观看 | 综合黄色网| 三级av在线免费观看 | 黄色a一级片 | 麻豆国产精品永久免费视频 | 在线观看亚洲精品 | 91看片在线免费观看 | 91看片在线播放 | 黄色网在线免费观看 | 天天操夜夜操夜夜操 | www.亚洲视频| 国产免费观看久久黄 | 美女网站视频色 | 色99网 | 亚洲欧美视频 | 国产黄在线免费观看 | 久操中文字幕在线观看 | 精品在线播放视频 | 69亚洲乱 | 91精品国产91p65| 极品国产91在线网站 | 91超级碰碰 | 亚洲成人精品国产 | 国内精品久久久久久久久久清纯 | 免费无遮挡动漫网站 | 国产一级一级国产 | 久久大片| 黄色成人影视 | 亚洲精品小视频 | 中文字幕一区二区三区久久 | 久久免费视频6 | av中文电影 | 亚洲国产小视频在线观看 | 精品爱爱 | 欧美激情视频一区二区三区 | 国产精品九九热 | 成人欧美一区二区三区黑人麻豆 | 亚洲激精日韩激精欧美精品 | 久久久精品视频网站 | 国产一二三精品 | 麻豆va一区二区三区久久浪 | 五月色丁香 | 欧美精品中文字幕亚洲专区 | 四虎在线观看精品视频 | 麻豆久久精品 | 人人爽人人舔 | 911香蕉视频| 99久久精品国产一区二区成人 | 国产精品一区免费观看 | 99久久精品免费看国产一区二区三区 | 在线免费视频你懂的 | 国产在线一区二区三区播放 | 少妇bbb搡bbbb搡bbbb | 九色激情网 | 最近中文字幕第一页 | 日韩美一区二区三区 | 综合网伊人 | 国产裸体视频网站 | 国产一级免费在线观看 | 国产成人福利在线 | 欧美三级高清 | 国产精品毛片一区二区三区 | 欧美片一区二区三区 | 国产小视频网站 | 日日成人网 | 伊人久久国产精品 | 正在播放国产91 | 日韩特级黄色片 | 色视频网站在线观看一=区 a视频免费在线观看 | 亚洲波多野结衣 | 亚洲影视九九影院在线观看 | aaa毛片视频 | 精品a视频| 狠狠色丁香婷婷综合橹88 | 91九色性视频 | 国产高清一 | 日韩最新av | av日韩精品 | 国产日韩欧美在线一区 | 久久国产经典 | 国内精品视频一区二区三区八戒 | 又粗又长又大又爽又黄少妇毛片 | 国产区第一页 | 日韩欧美在线视频一区二区 | 久久免费国产电影 | 人人讲 | 国产精品自产拍在线观看中文 | av在线观| 91色九色 | 亚洲黄色在线免费观看 | 久久国产手机看片 | 视频在线99re | 亚洲一二三久久 | 成人性生交视频 | 亚洲精品黄色 | 视频在线观看国产 | 人人爽久久久噜噜噜电影 | 免费网站黄色 | 久久精品79国产精品 | 亚洲免费永久精品国产 | 黄色片网站av | 国产精品久久久久久久久久新婚 | 最近中文字幕高清字幕在线视频 | 一级精品视频在线观看宜春院 | 97在线影院 | 久久精品久久精品久久精品 | 色九九在线| 久久久久久久久爱 | 五月婷婷播播 | 国产一级a毛片视频爆浆 | 国产成人久久久77777 | 香蕉影视 | 亚洲人人射 | 色吊丝在线永久观看最新版本 | 日韩夜夜爽| 欧美乱码精品一区二区 | 日日操天天操狠狠操 | 免费精品| 九九热精品视频在线播放 | 欧美精品久久久久久久久久 | 91日韩精品 | 狠狠操狠狠干2017 | 久久综合九色九九 | 国产在线观看国语版免费 | 狂野欧美激情性xxxx | 久久国内精品99久久6app | 亚洲精品18日本一区app | 2019中文最近的2019中文在线 | 日韩二区三区 | 97在线视频免费观看 | 在线观看免费黄色 | 波多野结衣在线观看一区二区三区 | 深爱激情综合 | 国产精品欧美久久久久天天影视 | 欧美午夜理伦三级在线观看 | 91桃色国产在线播放 | 91在线免费观看网站 | 视频精品一区二区三区 | 在线v| 开心丁香婷婷深爱五月 | 久久极品 | 97综合视频| 国产精品日韩欧美一区二区 | 久久国产高清视频 | 高清日韩一区二区 | 日韩精品视频在线免费观看 | 国产一级久久 | 97av视频在线观看 | 91麻豆高清视频 | 玖玖在线播放 | 色天天久久 | 免费看一级 | 91热爆视频| 成人av在线影院 | 国产拍揄自揄精品视频麻豆 | 麻豆免费观看视频 | 欧美精品久久久久久久 | 色播亚洲婷婷 | 色综久久 | 极品久久久久久久 | 中文字幕视频免费观看 | 国产视频日韩视频欧美视频 | 亚a在线 | 激情丁香久久 | 成人免费视频观看 | 男女男视频| 精品产品国产在线不卡 | 国产一级特黄毛片在线毛片 | 在线看黄色的网站 | 国产色女人 | 久草精品视频在线观看 | 国产一区二区在线免费播放 | 日本乱视频 | 日韩在线播放视频 | 国产一区二区不卡视频 | 丁香久久久 | 成人午夜在线电影 | 狠狠久久婷婷 | 日韩精品久久一区二区 | 久久久久久久综合色一本 | 欧美日韩视频在线观看免费 | av电影一区二区三区 | 国产麻豆果冻传媒在线观看 | 九九导航 | 国产精品视频全国免费观看 | 婷婷六月丁 | 久久99国产精品免费 | 在线超碰av| 91视频免费看片 | 一区二区三区高清在线观看 | 色综合久久88色综合天天人守婷 | 丁香婷婷久久久综合精品国产 | 中文在线免费观看 | 99精品国产99久久久久久福利 | 日韩二区三区在线 | 网站免费黄色 | 国产精品免费一区二区三区 | 欧美日韩国产一区 | 人人爽人人 | 精品久久国产精品 | 日日爽| 成人免费在线网 | 欧美日韩网址 | 国产成人精品一区二区三区在线观看 | 国产精品电影在线 | 国产精品福利视频 | 亚洲国产精品一区二区久久,亚洲午夜 | 天天插天天操天天干 | 国产精品com | 91日韩精品视频 | 免费男女羞羞的视频网站中文字幕 | 又黄又刺激视频 | 精品在线一区二区 | 国产精品高潮呻吟久久av无 | 中文字幕在线观看完整版 | 久久高清视频免费 | 国产精品久久久久久影院 | 亚洲综合网站在线观看 | 国产一级一片免费播放放a 一区二区三区国产欧美 | 麻豆视频在线观看免费 | 中文字幕在线观看一区二区三区 | 日韩免费在线网站 | 国产亚洲精品女人久久久久久 | 久久久国产一区二区 | 久久国产免费看 | 日本性动态图 | 国产精品国内免费一区二区三区 | 国产麻豆精品久久 | 五月婷婷视频在线 | 色在线免费观看 | 午夜视频在线观看一区二区三区 | 欧美日韩免费在线观看视频 | 亚洲精品av在线 | 久久精品毛片 | 久草在线视频首页 | 偷拍精偷拍精品欧洲亚洲网站 | 精品一区精品二区高清 | 国产精品久久久久久一区二区三区 | 精品国产视频在线观看 | 国产明星视频三级a三级点| 久久综合桃花 | 成人黄色免费在线观看 | 国产精品视频久久久 | 伊人电影在线观看 | 岛国一区在线 | 欧美国产精品一区二区 | 高清精品在线 | 天天射一射 | 一级久久精品 | 久久五月天色综合 | 久久毛片高清国产 | 亚洲黄色三级 | 中文字幕第| 黄色精品一区 | 日韩在线电影 | 久久精品免费播放 | 日本h视频在线观看 | 夜夜操天天操 | 久久理论视频 | 国产精品mv | 久久av在线播放 | 99久久夜色精品国产亚洲 | 欧美高清成人 | 五月天堂网 | 黄色国产成人 | 婷婷久操 | 亚洲第一香蕉视频 | av字幕在线 | 在线观看久草 | 国产伦精品一区二区三区照片91 | 福利视频第一页 | 国产午夜精品一区二区三区 | 天天狠狠 | 在线观看成人国产 | 久久久久久久久久久网站 | 亚洲国产经典视频 | 国产精品久久中文字幕 | www.久久免费视频 | 日韩二区在线 | 99精品黄色片免费大全 | 国产精品理论视频 | av日韩精品 | 精品国产综合区久久久久久 | 九九九电影免费看 | 在线观看av大片 | 国产中文字幕第一页 | 亚洲一区二区精品在线 | 中文久久精品 | 国产精品国产三级在线专区 | 亚洲无线视频 | 精品欧美在线视频 | 欧美日韩免费视频 | 日韩在线| 国产精品久久久久久久久久直播 | 97av影院 | 国产高清在线视频 | 99中文字幕视频 | 久久久综合 | 日日摸日日爽 | 在线观看深夜福利 | 国产91在线免费视频 | 日本精品一二区 | 国产码电影 | 亚洲激情视频 | 日韩一区在线免费观看 | 精品久久久久久亚洲综合网站 | 国产精选视频 | 8x成人免费视频 | 欧美精品久久久久性色 | 黄色软件视频大全免费下载 | 国产精品久久精品国产 | 丁香六月伊人 | 超级碰碰碰免费视频 | 欧美激情综合五月色丁香 | 黄色片毛片 | 色婷婷综合成人av | 色综合久久88色综合天天6 | 日韩av视屏| 国产精品9区 | 国产破处精品 | 久久精品亚洲国产 | 国产91九色蝌蚪 | 国产精品美女www爽爽爽视频 | 中文字幕一区二区三区四区在线视频 | 成人九九视频 | 免费观看高清 | 日韩在线观看中文字幕 | 婷婷av电影 | 欧美a性 | 激情在线网址 | av成人免费 | 国产三级视频 | 久久国产精品99久久久久久老狼 | 日韩在线观看你懂的 | 97夜夜澡人人双人人人喊 | av一本久道久久波多野结衣 | 国产一级视屏 | 日本韩国精品一区二区在线观看 | 人人爽久久涩噜噜噜网站 | 99免费在线观看视频 | 亚洲午夜久久久久久久久久久 | 在线电影 一区 | 一区二区不卡高清 | 日日夜夜精品 | 日韩在线视频免费播放 | 91色偷偷| 欧美激情va永久在线播放 | 久久黄色免费观看 | 黄色一集片 | 久久夜色电影 | 国产精品美女免费视频 | 精品国产精品久久 | 十八岁免进欧美 | 在线观看国产一区二区 | 狠狠干天天操 | 久久99视频免费观看 | 美女视频黄频大全免费 | 日本激情视频中文字幕 | 久久黄色免费观看 | av一级网站 | 91麻豆精品国产91久久久久久久久 | 久久精品婷婷 | 激情久久伊人 | 精品电影一区 | 亚洲国产精品一区二区久久hs | 在线观看91精品国产网站 | 久久精品国产成人 | 久久久网址 | 久久久国产影院 | 五月开心激情 | 久久久国产精品一区二区中文 | 免费观看黄色av | 99免费在线观看视频 | 激情狠狠干 | 黄色app网站在线观看 | 国产在线一线 | 久久9精品 | 国产成人久久av | 99久久精品国产免费看不卡 | 久草视频在线免费看 | 国产一区二区午夜 | 精品视频免费观看 | 日韩欧美精品在线观看视频 | 久久99国产精品二区护士 | 日韩精品免费一区二区三区 | 中文字幕视频网站 | 日韩一区二区三区在线观看 | 91亚洲精品国偷拍自产在线观看 | 日韩欧美在线高清 | 人人天天夜夜 | 2020天天干天天操 | 偷拍福利视频一区二区三区 | www.干| 日韩中文三级 | 日韩欧美一区二区在线 | 日韩欧美电影网 | 人人超碰免费 | 亚洲精品视频一 | 日日夜夜精品视频天天综合网 | 国产精品 日韩 | 91麻豆网站 | 久久综合一本 | 欧美了一区在线观看 | 天堂av网址 | 久久精品美女视频 | 欧美一区在线观看视频 | 黄色不卡av | 精品高清美女精品国产区 | 波多野结衣精品视频 | 久草色在线观看 | 免费观看性生交大片3 | 成人一区二区在线观看 | 日韩电影中文 | 日韩欧美精品一区 | 视频三区 | 91桃色国产在线播放 | 日批网站在线观看 | 狠狠综合网 | 精品国产一区二区三区四区在线观看 | 很黄很黄的网站免费的 | 久操操| 亚洲一区 影院 | 成人午夜av电影 | 亚洲精品毛片一级91精品 | 国产精品久久久久四虎 | 天天干天天拍天天操天天拍 | 最近中文字幕国语免费高清6 | 久久高视频 | 亚洲国产精品成人女人久久 | 韩日精品中文字幕 | av888.com| 狠狠色伊人亚洲综合网站野外 | 精品一区久久 | 日韩av成人在线 | 国产成人一区三区 | 久久九九久久 | 久要激情网 | 免费网站看v片在线a | 国产黑丝一区二区 | 亚洲精品国产精品国 | 激情av网址 | 久久久伦理| 国产福利在线免费观看 | 偷拍精偷拍精品欧洲亚洲网站 | 国产精品美女网站 | 27xxoo无遮挡动态视频 | 人人爽人人爽人人爽 | 久久久性 | 亚洲国产精品99久久久久久久久 | 插久久 | 丰满少妇高潮在线观看 | 国产理论片在线观看 | 一级黄色大片 | 黄色福利网| 亚洲午夜精品福利 | av一级免费 | 日韩在线观看一区二区三区 | 中文免费在线观看 | 涩涩色亚洲一区 | 999电影免费在线观看 | 97超碰人人爱 | 久久99精品久久久久蜜臀 | 韩日电影在线 | 在线视频国产区 | 91视频高清 | 婷婷色社区 | 国产精久久久久久妇女av | 久久国产三级 | 亚洲黄色小说网址 | 国产在线a视频 | 久久高清毛片 | 夜夜躁狠狠躁日日躁视频黑人 | 处女av在线 | 久久精品久久精品久久精品 | wwxxxx日本| 黄色大片国产 | 欧美性生活免费看 | 欧美日韩不卡一区二区三区 | 成人毛片在线观看视频 | 国产一区二区在线观看免费 | 日韩专区一区二区 | 人人看人人爱 | 777xxx欧美 | 国产精品一区二区在线 | 色婷婷综合久久久中文字幕 | 国产在线p | av九九九| 日韩一区二区三区免费视频 | 成人小视频在线 | 精品久久久久国产 | 亚洲天天干| 国产成人免费av电影 | 精品一区 在线 | 久久高清毛片 | 九九综合久久 | 天天爽天天射 | 久草在线这里只有精品 | 在线观看久久久久久 | 国产精品国产三级国产不产一地 | 最新av免费 | 欧美性色黄 | 国产一区二区精品91 | 操一草 | 中文字幕丰满人伦在线 | 91成人在线免费观看 | 天天天天爽 | 久久精品国产一区二区三 | 黄免费在线观看 | 久操中文字幕在线观看 | 中文字幕av全部资源www中文字幕在线观看 | 久草在线视频中文 | 韩国三级av在线 | 最新av在线免费观看 | 91人人射 | 免费男女羞羞的视频网站中文字幕 | 久久毛片视频 | www.com黄| 韩日电影在线 | 91在线最新 | 99热国产在线 | 色噜噜日韩精品一区二区三区视频 | 91亚洲视频在线观看 | 国产又粗又猛又黄 | 国产一级二级av | 久久国产a| 日本一区二区高清不卡 | 在线免费视频 你懂得 | 国偷自产视频一区二区久 | 国产午夜精品一区二区三区在线观看 | 天天做天天爱天天爽综合网 | 欧美在线视频第一页 | 欧美性色综合网 | 69热国产视频 | 日本福利视频在线 | 久久99国产综合精品免费 | 亚洲97在线 | 成人av免费在线看 | 久久www免费人成看片高清 | 日韩欧美综合视频 | 国产1区2 | 日韩超碰 | 午夜精品久久久久久久99 | 日韩在线免费小视频 | 欧美一级特黄aaaaaa大片在线观看 | 91高清免费看 | 日韩网页| 亚洲国产精品人久久电影 | 日韩欧美电影网 | 久久黄色免费观看 | 超级碰碰碰免费视频 | 亚洲理论在线观看 | 国产精品久久久久久久妇 | 最新日韩在线观看 | 久久婷婷五月综合色丁香 | 精品久久电影 | 四虎海外影库www4hu | 免费亚洲一区二区 | 国产视频69 | 婷婷网址 | 亚洲精品久久久蜜桃直播 | 视频一区二区视频 | 精品一区二区在线看 | 亚洲精品在线免费播放 | 黄色软件网站在线观看 | 特级西西444www高清大视频 | 久久综合狠狠综合 | 中文字幕日韩免费视频 | 国产不卡在线观看 | 亚洲婷久久 | 久久成人亚洲欧美电影 | 久草视频资源 | 精品二区视频 | 久草在线播放视频 | 国产高清中文字幕 | 一区二区三区 中文字幕 | 免费av片在线 | 国产亚洲精品久久久久久久久久久久 | 91麻豆精品国产 | 999久久久久久久久 69av视频在线观看 | 亚洲综合成人在线 | 国产二区免费视频 | 午夜影视一区 | 人人干人人爽 | 免费av高清 | 91亚洲精品久久久 | 福利视频网站 | 亚洲精品av中文字幕在线在线 | 美女视频黄在线 | 2022久久国产露脸精品国产 | 国产日韩在线看 | 91看片网址| 99久久精品免费看国产免费软件 | 欧美精品小视频 | 国产美女视频黄a视频免费 久久综合九色欧美综合狠狠 | 国产视频在线免费 | 天天干,夜夜爽 | 精品福利视频在线观看 | 免费色网站 | 日本mv大片欧洲mv大片 | 西西人体www444| 久草免费手机视频 | 日本成址在线观看 | 五月天六月丁香 | a级片网站 | 人人澡人摸人人添学生av | 亚洲电影免费 | 中文字幕黄色 | 欧美精品视 | 亚洲电影一区二区 | 综合网婷婷 | 婷婷中文在线 | 免费在线观看成人av | 91精品欧美一区二区三区 | 福利电影久久 | 丁香婷婷综合五月 | 999久久久久久久久6666 | 91丨九色丨蝌蚪丨对白 | 成人免费在线播放视频 | 九九免费精品 | 国产视频九色蝌蚪 | 九九视频在线观看视频6 | 欧美专区亚洲专区 | 久久99精品久久久久久三级 | 国产精品精 | 成人午夜电影久久影院 | 蜜臀av性久久久久蜜臀aⅴ四虎 | 欧美婷婷综合 | 69国产精品视频免费观看 | 婷婷午夜天 | 国内外成人免费在线视频 | 日本h视频在线观看 | 欧美一区二区三区免费看 | 久久免费大片 | 欧美二区视频 | 久久亚洲成人网 | 天天操天天操天天操天天操 | 免费观看黄色12片一级视频 | 国产视频九色蝌蚪 | 激情小说久久 | 最近更新好看的中文字幕 | 欧美老女人xx| 丁香婷婷激情 | 国产一区二区中文字幕 | 99精品视频在线播放观看 | 狠狠色丁香婷婷综合基地 | 香蕉视频18 | 精品嫩模福利一区二区蜜臀 | 日韩精品免费一线在线观看 | 91av中文| 精品久久久一区二区 | 久久艹在线观看 | 色97在线 | 国产精品嫩草55av | 在线观看完整版 | 亚洲伊人网在线观看 | 天天操人人要 | 亚洲精品乱码白浆高清久久久久久 | av在线播放不卡 | 亚洲天天综合网 | 久久精品综合网 | 日韩av在线免费播放 | 免费一级片在线观看 | 亚洲精品视频免费观看 | 日本中文乱码卡一卡二新区 | av一本久道久久波多野结衣 | 国产精品一区二区三区在线播放 | 伊人www22综合色 | 在线视频中文字幕一区 | 久久在线观看视频 | 久产久精国产品 | 午夜视频一区二区三区 | 国产大片免费久久 | 亚洲在线视频网站 | 欧美一级性生活片 | 五月婷婷在线观看视频 | 成人午夜在线观看 | 亚洲精品视频免费在线观看 | 亚洲综合导航 | 手机在线黄色网址 | 久久久亚洲成人 | 免费看黄色小说的网站 | 国产在线不卡视频 | 久久久精品国产一区二区 | 成人a级黄色片 | 97超在线视频| 久久黄色影视 | 一区精品在线 | 免费在线观看午夜视频 | 99r在线精品 | 国产在线观看你懂的 | 91爱爱视频 | 色综合久久久网 | 国产欧美精品一区二区三区四区 | 狠狠色丁香九九婷婷综合五月 | 亚洲免费在线 | а中文在线天堂 | 黄在线免费看 | 日韩丝袜视频 | 久久精品精品电影网 | 欧美巨大荫蒂茸毛毛人妖 | 人人狠狠综合久久亚洲 |