精华内容
下载资源
问答
我要提问
  • 多元正态分布的参数估计 可用于多元正态的参数估计 及 均值向量和协差阵检验 的R语言函数总结 千次阅读 多人点赞
    2021-12-24 23:28:19

    一、多元正态的参数估计

    1.1 样本均值

            在R语言中，均值通常用函数mean()得到，但是mean()只能计算一维变量的样本均值，在面对多元随机变量的样本时，假设我们以数据框的形式保存样本，我们有以下方法可以得到样本均值：

    1. 对多元样本的每一个分量用mean()函数，可以用apply()或sapply()函数
    2. 以数据框类型保存的样本，可以用summary()函数返回各个变量的各项描述性数据，其中包括均值

    例1.1：

            计算有割草机家庭的样本均值向量

    有割草机家庭无割草机家庭
    x1x2x1x2
    20.09.225.09.8
    28.58.417.610.4
    21.610.821.68.6
    20.510.414.410.2
    29.011.828.08.8
    36.79.616.48.8
    36.08.819.88.0
    27.611.222.09.2
    23.010.015.88.2
    31.010.411.09.4
    17.011.017.07.0
    27.010.021.07.4

     注：在输入数据时，通常会用一个新的变量(假设命名为Y)来表示每个观测所属的组，称为分组变量，这个变量在R中通常要转换成因子

    > data4.1=read.csv('Table4-1.csv')
> head(data4.1)
###                y是分组变量，y=1表示有割草机家庭
     x1   x2 y
1  20.0  9.2 1
2  28.5  8.4 1
3  21.6 10.8 1
4  20.5 10.4 1
5  29.0 11.8 1
6  36.7  9.6 1
###
> apply(data4.1[data4.1$y==1,],2,mean)[1:2]    #用apply()函数运算
###
      x1       x2 
26.49167 10.13333 
###
> summary(data4.1[data4.1$y==1,])           #用summary()获取各分量的样本均值
###
       x1              x2              y    
 Min.   :17.00   Min.   : 8.40   Min.   :1  
 1st Qu.:21.32   1st Qu.: 9.50   1st Qu.:1  
 Median :27.30   Median :10.20   Median :1  
 Mean   :26.49   Mean   :10.13   Mean   :1   #该行为均值
 3rd Qu.:29.50   3rd Qu.:10.85   3rd Qu.:1  
 Max.   :36.70   Max.   :11.80   Max.   :1 
###
    • apply()用法：apply(A,margin,fun,...) 

            apply()函数用来对矩阵或数据框的每行或每列进行指定函数的运算。其中A为矩阵或数据框；margin指定对行或对列进行运算，当margin=1时对行进行运算，当margin=2时对列进行运算；fun是指定的函数 

    • summary()用法：summary(object,...)

            summary()多用于获取项目的摘要，包含部分信息。当object为数据框时，会返回各个变量的五数（最小值，下四分位数，中位数，上四分为数，最大值）和均值

    1.2 样本协差阵

            在R中，样本协差阵的获取非常简便， 对数据框使用cov()函数即可

    例1.2：

            继上题，计算有割草机组的样本协差阵

    > cov(data4.1[data4.1$y==1,][,1:2])
###
          x1        x2
x1 39.182652 -1.969697
x2 -1.969697  1.020606
###
    • cov()用法：cov(x,y=NULL,...)

            当指定cov()的参数x和y，且两者都为一维向量时，会返回两个向量的样本协方差；而未指定参数y，且x为矩阵或数据框时，会返回以x每一列作为变量样本的协差阵

    1.3  样本相关阵

            获取样本相关阵的函数是cor()，其用法与cov()相同，两个一维向量返回相关系数；数据框返回相关阵

    二、各类检验

    2.1 正态性检验

            正态性检验即检验样本是否来自正态总体的检验，原假设都为来自正态总体。正态性的检验方法有许多种，此处介绍小样本量(3~50)时所用的夏皮洛-威尔克检验。R中的夏皮洛-威尔克检验的函数为shapiro.test()

            shapiro.test()一次只能对一维变量进行正态性检验，当面对多元随机变量的样本时，有以下方法

    1. 我们可以对其每一个分量都进行一次正态性检验，当所有分量都检验得出服从正态分布后，可以认为该多元随机变量服从多元正态分布
    2. 运用mvnormtest包内的mshapiro.test()函数进行多元正态性检验

            实现时可能会用到的函数有：

    • sapply()，对每个分量进行指定的检验
    • tapply()，对以分组变量指定的不同组别分别进行指定的检验

    例2.1：

            继上题，对不同类型家庭的随机向量数据进行正态性检验

    > sapply(data4.1[,-3],shapiro.test)    #对各分量进行正态性检验，但是未分组
###
          x1                            x2                           
statistic 0.9654387                     0.9880936                    
p.value   0.5568611                     0.9897171                    
method    "Shapiro-Wilk normality test" "Shapiro-Wilk normality test"
data.name "X[[i]]"                      "X[[i]]" 
###
> tapply(data4.1[,1],data4.1$y,shapiro.test)
###                            对分组后的x1进行正态性检验
$`0`

        Shapiro-Wilk normality test

data:  X[[i]]
W = 0.98551, p-value = 0.9971


$`1`

        Shapiro-Wilk normality test

data:  X[[i]]
W = 0.95332, p-value = 0.6859

###
> tapply(data4.1[,2],data4.1$y,shapiro.test)
###                            对分组后的x2进行正态性检验
$`0`

        Shapiro-Wilk normality test

data:  X[[i]]
W = 0.97557, p-value = 0.9596


$`1`

        Shapiro-Wilk normality test

data:  X[[i]]
W = 0.98262, p-value = 0.992
###
##对有割草机家庭的随机向量数据进行正态性检验
> mshapiro.test(t(as.matrix(data4.1[data4.1$y==1,-3])))

        Shapiro-Wilk normality test

data:  Z
W = 0.96877, p-value = 0.8975

##对无割草机家庭的随机向量数据进行正态性检验
> mshapiro.test(t(as.matrix(data4.1[data4.1$y==0,-3])))

        Shapiro-Wilk normality test

data:  Z
W = 0.98001, p-value = 0.9837

             检验结果为均服从正态分布

    • sapply()：sapply(X,Fun,...)

            sapply()用于对X的每个分量进行Fun函数运算，X应该是矩阵或数据框

    • tapply()：tapply(X,Index,Fun=NULL,...)

            tapply()用于对以分组变量Index指示的每个组中对应的X的数据进行Fun函数运算

    • mshapiro.test()：mshapiro.test(U)

            mshapiro.test()用于进行多元的夏皮洛-威尔克正态性检验，需要注意U只能是数据矩阵，当遇到用数据框存储的数据时要用as.matrix()转化为矩阵，且这个函数默认变量的数据按行排放，通常我们需要对矩阵再进行一次转置

            另外，可以画出Q-Q图查看样本的正态性，常用的函数有qqnorm()和qqline()

    • qqnorm(x)，其中x为一维变量的样本，当画出的散点图越接近一条斜线，其正态性越强
    • qqline(x)，其中x为一维变量的样本，当画出散点的Q-Q图后，添加点所靠近的斜线，该斜线的斜率为标准差,截距为均值

    2.2  均值向量的检验

            一维的均值检验有很多，若样本服从正态分布我们可以用t.test()单个总体或双总体的t检验；若不服从正态分布，我们可以用wilcox.test()进行秩和检验，用法与t.test()类似；当遇到多个总体时，若各个变量的方差相差不大，我们可以用将各个变量的数据放到一列，然后用一个分组变量表示数据属于哪个变量，运用aov()进行方差分析，从而进行多总体的均值检验

            当遇到多元随机变量的均值检验时，我们有以下方法：

    1. 对每个分量进行均值检验，通过正态性检验的用t检验，未通过正态性检验的用秩和检验
    2. 对通过多元正态检验的数据，运用ICSNP包中的HotellingsT2()函数进行均值向量的检验
    3. 多总体时，若协差阵齐性检验通过，可以用manova()进行多元方差分析

     例2.2.1：

            继上题，检验有割草机家庭和无割草机家庭的向量均值是否相同

    #上题已得出各类家庭的数据均通过正态检验

> t.test(x1~y,data=data4.1)            #对x1分量进行t检验

        Welch Two Sample t-test

data:  x1 by y
t = -3.2508, df = 20.458, p-value = 0.003919
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -12.073229  -2.643437
sample estimates:
mean in group 0 mean in group 1 
       19.13333        26.49167 

> t.test(x2~y,data=data4.1)            #对x2分量进行t检验

        Welch Two Sample t-test

data:  x2 by y
t = -3.1203, df = 21.956, p-value = 0.004991
alternative hypothesis: true difference in means between group 0 and group 1 is not equal to 0
95 percent confidence interval:
 -2.1918725 -0.4414609
sample estimates:
mean in group 0 mean in group 1 
       8.816667       10.133333 

> library(ICSNP)
> attach(data4.1)
> HotellingsT2(cbind(x1,x2)~y)            #霍特林T方检验，用于多元正态向量

        Hotelling's two sample T2-test

data:  cbind(x1, x2) by y
T.2 = 12.257, df1 = 2, df2 = 21, p-value = 0.000297
alternative hypothesis: true location difference is not equal to c(0,0)

            对各分量进行t检验的结果是：有割草机家庭与无割草机家庭的两个分量都不相等 

            对向量整体进行霍特林T方检验的结果是：有割草机家庭与无割草机家庭的向量均值不相等

             该题也可以用多元方差分析，只是水平数为2，其实也是通过对每个分量进行检验

    例2.2.2：

            现有如下表所示各省统计数据，试检验它们的均值向量是否等于  (1081,1822,115,179)

    序号

    省份

    工资性收入

    家庭性收入

    财产性收入

    转移性收入

    1

    北京

    4524.25

    1778.33

    588.04

    455.64

    2

    天津

    2720.85

    2626.46

    152.88

    79.64

    3

    河北

    1293.50

    1988.58

    93.74

    105.81

    4

    山西

    1177.94

    1563.52

    62.70

    86.49

    5

    内蒙古

    504.46

    2223.26

    73.05

    188.10

    6

    辽宁

    1212.20

    2163.49

    113.24

    201.28

    注：以上为部分表格，表格全部内容在此不展示

    > data3=read.csv('Table_0.csv',encoding='UTF-8')    #读取数据
> mu_bar=c(1081,1882,115,179)
> rownames(data3)=data3[,1]            #将省名赋给数据框行名
> data3=data3[,-1]                     #去除省名一列
###
假设通过了正态性检验
###
> HotellingsT2(data3,mu=mu_bar)

        Hotelling's one sample T2-test

data:  data3
T.2 = 1.8443, df1 = 4, df2 = 27, p-value = 0.1494
alternative hypothesis: true location is not equal to c(1081,1882,115,179)

            检验结果为各省统计数据的均值向量等于(1081,1822,115,179)

    例2.2.3：

    在数据New drug.xls中，各变量的意义为drug（药），取值1表示对病人给以新药，取值2表示对病人给以安慰剂，resp1-resp3是治疗后病人三个时点的呼吸状况，pulse1-pulse3是病人三个时点的脉搏。试分析这两方法的各次重复测定均值向量是否有显著差异？

    drugresp1resp2resp3pulse1pulse2pulse3
    13.43.33.32.22.12.1
    13.43.43.32.22.12.2
    13.33.43.42.32.42.3
    23.33.33.32.82.92.7
    23.23.33.42.62.72.7
    23.23.23.22.72.92.7

     注：以上为部分表格，表格全部内容在此不展示

            题目要求检验不用用药的组之间，向量(resp1,resp2,resp3,pulse1,pulse2,pulse3)的均值是否相等。因为drug只有2个水平，可以对每个分量进行t检验，但是分量比较多会比较麻烦；也可以用多元方差分析，查看结果也是对每个分量的检验，不过需要先进行协差阵检验；用霍特林T2检验会比较简单。

    > data_drug=read.csv('new drug.csv',encoding='UTF-8')
> names(data_drug)[1]='drug'        #UTF-8格式的csv文件读取后，第一列的名字会有变动，此处改回
> attach(data_drug)
###
对每个变量进行正态性检验后
得知随机向量不服从多元正态分布
因此不能用t检验和霍特林T方检验，不过可以对每个分量进行秩和检验
假设数据通过了协差阵检验
接下来进行多元方差分析
###
> modle_drug=manova(cbind(resp1,resp2,resp3,pulse1,pulse2,pulse3)~drug,data=data_drug)
> summary.aov(modle_drug)
 Response resp1 :
            Df   Sum Sq  Mean Sq F value   Pr(>F)   
drug         1 0.040833 0.040833  14.412 0.003507 **
Residuals   10 0.028333 0.002833                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response resp2 :
            Df   Sum Sq  Mean Sq F value   Pr(>F)   
drug         1 0.040833 0.040833  14.412 0.003507 **
Residuals   10 0.028333 0.002833                    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response resp3 :
            Df   Sum Sq   Mean Sq F value Pr(>F)
drug         1 0.020833 0.0208333  3.0488 0.1114
Residuals   10 0.068333 0.0068333               

 Response pulse1 :
            Df  Sum Sq Mean Sq F value    Pr(>F)    
drug         1 0.65333 0.65333      70 7.936e-06 ***
Residuals   10 0.09333 0.00933                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response pulse2 :
            Df  Sum Sq Mean Sq F value    Pr(>F)    
drug         1 1.08000 1.08000  79.024 4.623e-06 ***
Residuals   10 0.13667 0.01367                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

 Response pulse3 :
            Df  Sum Sq Mean Sq F value    Pr(>F)    
drug         1 0.75000 0.75000  64.286 1.155e-05 ***
Residuals   10 0.11667 0.01167                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

            只有resp3的检验结果为相同，其他都为不同，所以认为两方法的各次重复测定均值向量有显著差异

    • t.test()用法:

    t.test(x,y=NULL,alternative=c("two.sided","less","greater"),mu=0,...)或

    t.test(formula,data,...)

    指定x未指定y时，进行单总体的t检验，可以指定mu检验其是否与mu相同

    alternative指定双边检验或左尾检验或右尾检验

    也可以用formula类的参数，即x~y的类型，y是分组变量，会对不同组的x进行双总体t检验

    • wilcox.test()用法：

    wilcox.test()用法与t.test()相同，此处不赘述

    • HotellingsT2()用法:

    HotellingsT2(X,Y=NULL,mu=NULL,test="f",...)或

    HotellingsT2(formula,...)

    X与Y为矩阵或数据框，未指定Y时进行单总体的检验，可以指定mu检验其是否与mu相同

    test参数指定近似统计量，默认为f，即F近似，可以指定"chi"，即卡方近似

    可以用formula类参数，与先前用法相同，但是HotellingsT2()没有data参数

    • manova()用法:

    manova(formula,data,...)

    manova()的formula参数用法aov()类似，manova()返回的是多元方差分析的模型，将其赋给某个变量，然后用aov.summary()函数可以看每个变量的检验

    2.3  协差阵检验

            在进行多元方差分析前需要进行协差阵齐性检验，协差阵检验可以用heplots包内的boxM()函数。

    例2.3：

            继有无割草机家庭数据，检验两组的协差阵是否有差异

    > boxM(data4.1[,-3],group=data4.1[,3])

        Box's M-test for Homogeneity of Covariance Matrices

data:  data4.1[, -3]
Chi-Sq (approx.) = 0.99346, df = 3, p-value = 0.8028

            检验结果为两组协差阵相同

    • boxM()用法：

    boxM(formula,data,...)或

    boxM(Y,group,...)

    formula类参数的用法与之前的函数相同

    Y是数据矩阵或数据框，group是指定的分组变量

    boxM()函数进行的是协差阵齐性检验，在分组变量的水平数大于2时也可以使用

    三、小结

    总结本文提到的函数和应用场景

    参数估计正态性检验
    函数应用场景函数应用场景
    mean()计算一维变量的样本均值shapiro.test()小样本正态性检验
    apply()对矩阵或数据框的行或列进行运算mshapiro.test()多元小样本正态性检验
    sapply()对矩阵或数据框的每个变量进行运算sapply()对每个变量进行指定运算或检验
    summary()对数据框使用时返回每个变量的统计描述tapply()对以分组变量指定的不同组别分别进行指定的运算或检验
    cov()获取协方差或协差阵qqnorm()画Q-Q图
    cor()获取相关系数或相关阵qqline()在Q-Q图中添加正态标准线
    均值向量检验协差阵检验
    函数应用场景函数应用场景
    t.test()正态样本的单双总体均值检验boxM()协差阵齐性检验
    wilcox.test()非正态样本的单双总体均值检验
    HotellingsT2()多元正态样本的单双总体均值检验
    aov()方差齐性情况下的方差分析
    manova()协差阵齐性下的多元方差分析
    aov.summary()获取方差分析模型的检验结果
    r语言
    更多相关内容

  • 【多元统计分析】04.多元正态分布的参数估计

     千次阅读 2020-10-25 18:46:36
    多元正态分布具有两个参数——均值向量与自协方差函数，与数理统计一样，可以用抽样的方式定义一些统计量对它们进行参数估计。在这里，我们使用极大似然估计的方法，用样本均值和样本离差阵对它们进行估计。

    文章目录

    四、多元正态分布的参数估计

    1.多元正态分布的估计量

    对于多元正态分布 N p ( μ , Σ ) N_p(\mu,\Sigma) Np​(μ,Σ)，其参数只有两个——均值向量 μ \mu μ与自协方差矩阵 Σ \Sigma Σ，要对其进行估计，就要从总体中抽取简单随机样本。记抽取样本的容量为 n n n，每一个样本分别是 X ( α ) = ( x α 1 , ⋯   , x α p ) X_{(\alpha)}=(x_{\alpha1},\cdots,x_{\alpha p}) X(α)​=(xα1​,⋯,xαp​)，将样本纵向排列，得到样本数据阵
    X = [ x 11 ⋯ x 1 p ⋮ ⋮ x n 1 ⋯ x n p ] . X=\begin{bmatrix} x_{11} & \cdots & x_{1p} \\ \vdots & & \vdots \\ x_{n1} & \cdots & x_{np} \end{bmatrix}. X=⎣⎢⎡​x11​⋮xn1​​⋯⋯​x1p​⋮xnp​​⎦⎥⎤​.
    从样本数据阵出发，可以获得以下统计量：

    1. 样本均值 X ˉ \bar X Xˉ，这是对每个维度求均值，得到的一个 p p p维向量
      X ˉ = 1 n ∑ α = 1 n X ( α ) = ( x ˉ 1 , ⋯   , x ˉ p ) ′ = 1 n X ′ 1 n . \bar X=\frac 1n\sum_{\alpha=1}^n X_{(\alpha)}=(\bar x_1,\cdots ,\bar x_p)'=\frac 1nX'\boldsymbol 1_n. Xˉ=n1​α=1∑n​X(α)​=(xˉ1​,⋯,xˉp​)′=n1​X′1n​.
      这里 x ˉ i \bar x_i xˉi​是对第 i i i个分量的平均，即
      x ˉ i = 1 n ∑ α = 1 n x α i . \bar x_i=\frac 1n\sum_{\alpha=1}^n x_{\alpha i}. xˉi​=n1​α=1∑n​xαi​.

    2. 样本离差阵 A A A，可以类比一维随机变量中的 ∑ i = 1 n ( x i − x ˉ ) 2 \sum_{i=1}^n (x_i-\bar x)^2 ∑i=1n​(xi​−xˉ)2，即
      A = ∑ α = 1 n ( X ( α ) − X ˉ ) ( X ( α ) − X ˉ ) ′ A=\sum_{\alpha=1}^n(X_{(\alpha)}-\bar X)(X_{(\alpha)}-\bar X)' A=α=1∑n​(X(α)​−Xˉ)(X(α)​−Xˉ)′
      这样， A A A是一个 p × p p\times p p×p对角阵，它的第 ( i , j ) (i,j) (i,j)元，其实就是
      a i j = ∑ α = 1 n ( x α i − x ˉ i ) ( x α j − x ˉ j ) . a_{ij}=\sum_{\alpha=1}^n (x_{\alpha i}-\bar x_i)(x_{\alpha j}-\bar x_j). aij​=α=1∑n​(xαi​−xˉi​)(xαj​−xˉj​).
      由此，还可以得到
      A = X ′ X − n X ˉ X ˉ ′ = X ′ [ I n − 1 n 1 n 1 n ′ ] X . A=X'X-n\bar X\bar X'=X'\left[I_n-\frac 1n\boldsymbol 1_n\boldsymbol 1_n' \right] X. A=X′X−nXˉXˉ′=X′[In​−n1​1n​1n′​]X.
      这个式子用来计算离差阵更为方便。

    3. 样本协方差阵 S S S，可以类比一维随机变量中的样本方差，即
      S = 1 n − 1 A , S=\frac 1{n-1}A, S=n−11​A,
      其 ( i , i ) (i,i) (i,i)元是变量 X i X_i Xi​的样本方差，即
      s i i = 1 n − 1 ∑ α = 1 n ( x α i − x ˉ i ) 2 . s_{ii}=\frac 1{n-1}\sum_{\alpha=1}^n (x_{\alpha i}-\bar x_i)^2. sii​=n−11​α=1∑n​(xαi​−xˉi​)2.
      类似一维中样本方差的定义，也有
      S ∗ = 1 n ∑ α = 1 n ( x α i − x ˉ i ) 2 . S^*=\frac 1n\sum_{\alpha=1}^n(x_{\alpha i}-\bar x_i)^2. S∗=n1​α=1∑n​(xαi​−xˉi​)2.

    4. 样本相关阵 R R R，自然是由样本相关系数 r i j r_{ij} rij​构成的 p × p p\times p p×p矩阵，即
      R = s i j s i i s j j = a i j a i i a j j . R=\frac{s_{ij}}{\sqrt{s_{ii}s_{jj}}}=\frac{a_{ij}}{\sqrt{a_{ii}a_{jj}}}. R=sii​sjj​ ​sij​​=aii​ajj​ ​aij​​.

    有了这些统计量，我们就可以对总体的参数 μ , Σ \mu,\Sigma μ,Σ进行估计，使用的方法是最大似然估计。

    2.最大似然估计

    最大似然估计指的是，以使获得样本的出现几率最大的那组参数估计量，作为参数的点估计量。与一元情形类似，可以建立似然函数的概念。使用拉直运算，对 V e c ( X ′ ) {\rm Vec}(X') Vec(X′)的密度函数建立似然函数，称为样本 X ( i ) X_{(i)} X(i)​的似然函数（对数似然函数）。
    L ( μ , Σ ) = ∏ α = 1 n 1 ( 2 π ) p / 2 ∣ Σ ∣ 1 / 2 exp ⁡ [ − 1 2 ( x ( α ) − μ ) ′ Σ − 1 ( x ( α ) − μ ) ] = 1 ( 2 π ) n p / 2 ∣ Σ ∣ n / 2 exp ⁡ [ − 1 2 ∑ α = 1 n ( x ( α ) − μ ) ′ Σ − 1 ( x ( α ) − μ ) ] l ( μ , Σ ) = − n p 2 ln ⁡ ( 2 π ) + n 2 ln ⁡ ∣ Σ − 1 ∣ − 1 2 ∑ α = 1 n ( x ( α ) − μ ) ′ Σ − 1 ( x ( α ) − μ ) \begin{aligned} L(\mu,\Sigma)=&\prod_{\alpha=1}^n \frac{1}{(2\pi)^{p/2}|\Sigma|^{1/2}}\exp\left[-\frac12(x_{(\alpha)}-\mu)'\Sigma^{-1}(x_{(\alpha)}-\mu) \right] \\ =&\frac{1}{(2\pi)^{np/2}|\Sigma|^{n/2}}\exp\left[-\frac12\sum_{\alpha=1}^n(x_{(\alpha)}-\mu)'\Sigma^{-1}(x_{(\alpha)}-\mu) \right]\\ l(\mu,\Sigma)=&-\frac{np}2\ln(2\pi)+\frac n2\ln |\Sigma^{-1}|-\frac12\sum_{\alpha=1}^n(x_{(\alpha)}-\mu)'\Sigma^{-1}(x_{(\alpha)}-\mu) \end{aligned} L(μ,Σ)==l(μ,Σ)=​α=1∏n​(2π)p/2∣Σ∣1/21​exp[−21​(x(α)​−μ)′Σ−1(x(α)​−μ)](2π)np/2∣Σ∣n/21​exp[−21​α=1∑n​(x(α)​−μ)′Σ−1(x(α)​−μ)]−2np​ln(2π)+2n​ln∣Σ−1∣−21​α=1∑n​(x(α)​−μ)′Σ−1(x(α)​−μ)​

    要求其极大似然估计，需要对矩阵 Σ \Sigma Σ，向量 μ \mu μ求导（参见矩阵微商），得
    d l ( μ , Σ ) d μ = 1 2 ∑ α = 1 n ( Σ − 1 + ( Σ − 1 ) ′ ) ( x ( α ) − μ ) = Σ − 1 ( ∑ α = 1 n ( x ( α ) − μ ) ) = n Σ − 1 ( X ˉ − μ ) . d l ( μ , Σ ) d Σ − 1 = − n 2 Σ − 1 2 ∑ α = 1 n ( x ( α ) − μ ) ( x ( α ) − μ ) ′ = − 1 2 ( n Σ − A ) . \frac{{\rm d}l(\mu,\Sigma)}{{\rm d}\mu}=\frac12\sum_{\alpha=1}^n(\Sigma^{-1}+(\Sigma^{-1})')(x_{(\alpha)}-\mu)=\Sigma^{-1}(\sum_{\alpha=1}^n(x_{(\alpha)}-\mu))=n\Sigma^{-1}(\bar X-\mu).\\ \frac{{\rm d}l(\mu,\Sigma)}{{\rm d}\Sigma^{-1}}=-\frac n2\Sigma-\frac12\sum_{\alpha=1}^n(x_{(\alpha)}-\mu)(x_{(\alpha)}-\mu)'=-\frac12(n\Sigma-A). dμdl(μ,Σ)​=21​α=1∑n​(Σ−1+(Σ−1)′)(x(α)​−μ)=Σ−1(α=1∑n​(x(α)​−μ))=nΣ−1(Xˉ−μ).dΣ−1dl(μ,Σ)​=−2n​Σ−21​α=1∑n​(x(α)​−μ)(x(α)​−μ)′=−21​(nΣ−A).
    所以
    μ ^ = X ˉ , Σ ^ = A n . \hat \mu=\bar X,\quad \hat\Sigma = \frac An. μ^​=Xˉ,Σ^=nA​.

    用到的矩阵微商结论：对于对称阵 A A A与列向量 x x x，有
    d ln ⁡ ∣ A ∣ d A = A − 1 , d x ′ A x d A = x x ′ , d x ′ A x d x = ( A + A ′ ) x . \frac{{\rm d}\ln |A|}{{\rm d}A}=A^{-1},\\ \frac{{\rm d}x'Ax}{{\rm d}A}=xx',\\ \frac{{\rm d}x'Ax}{{\rm d}x}=(A+A')x. dAdln∣A∣​=A−1,dAdx′Ax​=xx′,dxdx′Ax​=(A+A′)x.

    如果在已知 μ = μ 0 \mu=\mu_0 μ=μ0​的情况下，依照以上过程，就可以得到
    Σ ^ = 1 n ∑ α = 1 n ( x ( α ) − μ 0 ) ( x ( α ) − μ 0 ) ′ . \hat \Sigma=\frac{1}{n}\sum_{\alpha=1}^n(x_{(\alpha)}-\mu_0)(x_{(\alpha)}-\mu_0)'. Σ^=n1​α=1∑n​(x(α)​−μ0​)(x(α)​−μ0​)′.
    所以，我们要找到 ( μ , Σ ) (\mu,\Sigma) (μ,Σ)的估计，就需要计算 ( X ˉ , A ) (\bar X,A) (Xˉ,A)，接下来对它们进行性质讨论。

    3.最大似然估计的性质

    ( X ˉ , A ) (\bar X,A) (Xˉ,A)的分布具有类似一元统计中 X ˉ \bar X Xˉ和 S 2 S^2 S2的性质。

    定理：设 X ˉ \bar X Xˉ和 A A A分别是 p p p元正态总体 N p ( μ , Σ ) N_p(\mu,\Sigma) Np​(μ,Σ)的样本均值向量和样本离差阵，则有

    1. X ˉ ∼ N p ( μ , 1 n Σ ) \bar X\sim N_p(\mu,\frac1n\Sigma) Xˉ∼Np​(μ,n1​Σ)；
    2. A = d ∑ t = 1 n Z t Z t ′ A\stackrel {\rm d}=\sum\limits_{t=1}^n Z_tZ_t' A=dt=1∑n​Zt​Zt′​，其中 Z 1 , ⋯   , Z n − 1 Z_1,\cdots,Z_{n-1} Z1​,⋯,Zn−1​独立同 N p ( 0 , Σ ) N_p(0,\Sigma) Np​(0,Σ)分布；
    3. X ˉ \bar X Xˉ和 A A A相互独立；
    4. P { A > 0 } = 1 ⇔ n > p {\rm P}\{A>0\}=1\Leftrightarrow n>p P{A>0}=1⇔n>p。

    前三个性质的证明方式也与一元情况类似，设 X X X为从多元正态总体中抽取的 n × p n\times p n×p样本数据阵， Γ \Gamma Γ是 n n n阶正交阵，形式如同
    Γ = [ r 11 ⋯ r 1 n ⋮ ⋮ r ( n − 1 ) 1 ⋯ r ( n − 1 ) n 1 / n ⋯ 1 / n ] = ( r i j ) n × n . \Gamma=\begin{bmatrix} r_{11} & \cdots & r_{1n} \\ \vdots & & \vdots \\ r_{(n-1)1} & \cdots & r_{(n-1)n} \\ 1/\sqrt n & \cdots & 1/\sqrt n \end{bmatrix}=(r_{ij})_{n\times n}. Γ=⎣⎢⎢⎢⎡​r11​⋮r(n−1)1​1/n ​​⋯⋯⋯​r1n​⋮r(n−1)n​1/n ​​⎦⎥⎥⎥⎤​=(rij​)n×n​.
    令
    Z = [ Z 1 ′ ⋮ Z n ′ ] = Γ [ X ( 1 ) ′ ⋮ X ( n ) ′ ] = Γ X . Z=\begin{bmatrix} Z_1' \\ \vdots \\ Z_n' \end{bmatrix} = \Gamma\begin{bmatrix} X_{(1)}' \\ \vdots \\ X_{(n)}' \end{bmatrix}=\Gamma X. Z=⎣⎢⎡​Z1′​⋮Zn′​​⎦⎥⎤​=Γ⎣⎢⎡​X(1)′​⋮X(n)′​​⎦⎥⎤​=ΓX.
    即 Z i ′ = ( r i 1 , ⋯ r i n ) X Z_i'=(r_{i1},\cdots r_{in})X Zi′​=(ri1​,⋯rin​)X，
    Z i = ( X ( 1 ) , ⋯   , X ( n ) ) [ r i 1 ⋮ r i n ] , i = 1 , ⋯   , n . Z_i=(X_{(1)},\cdots,X_{(n)})\begin{bmatrix} r_{i1} \\ \vdots \\ r_{in} \end{bmatrix},\quad i=1,\cdots,n. Zi​=(X(1)​,⋯,X(n)​)⎣⎢⎡​ri1​⋮rin​​⎦⎥⎤​,i=1,⋯,n.
    因为 Z i Z_i Zi​是 X ( 1 ) , ⋯   , X ( n ) X_{(1)},\cdots,X_{(n)} X(1)​,⋯,X(n)​的线性组合，所以 Z i Z_i Zi​也是 p p p维正态向量，且
    E Z i = ∑ α = 1 n r i α E ( X ( α ) ) = { n ∑ α = 1 n r i α r n α μ = 0 , t ≠ n ; ∑ α = 1 n 1 n μ = n μ , t = n . C o v ( Z α , Z β ) = ∑ i = 1 n r α i r β i Σ = { O , α ≠ β ; Σ , α = β . {\rm E}Z_i=\sum_{\alpha=1}^n r_{i\alpha}{\rm E}(X_{(\alpha)})=\left\{ \begin{array}l \sqrt{n}\sum\limits_{\alpha=1}^n r_{i\alpha}r_{n\alpha}\mu=0,&t\ne n;\\ \sum\limits_{\alpha=1}^n \frac 1{\sqrt n}\mu=\sqrt n \mu,&t=n. \end{array} \right.\\ {\rm Cov}(Z_\alpha,Z_{\beta})=\sum_{i=1}^nr_{\alpha i}r_{\beta i}\Sigma=\left\{ \begin{array}l O,&\alpha\ne \beta;\\ \Sigma,&\alpha=\beta. \end{array} \right. EZi​=α=1∑n​riα​E(X(α)​)=⎩⎪⎨⎪⎧​n ​α=1∑n​riα​rnα​μ=0,α=1∑n​n ​1​μ=n ​μ,​t​=n;t=n.​Cov(Zα​,Zβ​)=i=1∑n​rαi​rβi​Σ={O,Σ,​α​=β;α=β.​
    而显然 Z n = n X ˉ Z_n=\sqrt n\bar X Zn​=n ​Xˉ，且 Z n ∼ N p ( n μ , Σ ) Z_n\sim N_p(\sqrt n\mu,\Sigma) Zn​∼Np​(n ​μ,Σ)，所以 X ˉ ∼ N p ( μ , Σ / n ) \bar X\sim N_p(\mu,\Sigma/n) Xˉ∼Np​(μ,Σ/n)。而
    ∑ α = 1 n Z α Z α ′ = ( Z 1 , ⋯   , Z n ) [ Z 1 ⋮ Z n ] = Z ′ Z = X ′ X , ∑ α = 1 n − 1 Z α Z α ′ = X ′ X − Z n Z n ′ = X ′ X − n X ˉ X ˉ ′ = A . \sum_{\alpha=1}^nZ_{\alpha}Z_{\alpha}'=(Z_1,\cdots,Z_n)\begin{bmatrix} Z_1\\ \vdots \\ Z_n \end{bmatrix}=Z'Z=X'X,\\ \sum_{\alpha=1}^{n-1}Z_{\alpha}Z_{\alpha}'=X'X-Z_nZ_n'=X'X-n\bar X\bar X'=A. α=1∑n​Zα​Zα′​=(Z1​,⋯,Zn​)⎣⎢⎡​Z1​⋮Zn​​⎦⎥⎤​=Z′Z=X′X,α=1∑n−1​Zα​Zα′​=X′X−Zn​Zn′​=X′X−nXˉXˉ′=A.
    可以注意到， A A A是 Z 1 , ⋯   , Z n − 1 Z_1,\cdots,Z_{n-1} Z1​,⋯,Zn−1​的函数， X ˉ \bar X Xˉ是 Z n Z_n Zn​的函数，又因为 Z 1 , ⋯   , Z n Z_1,\cdots,Z_n Z1​,⋯,Zn​互相独立，所以 X ˉ \bar X Xˉ和 A A A相互独立。至于第四个性质，只需要记住，样本够多就能保证 A A A的非负定性即可。

    除此以外， X ˉ , A \bar X,A Xˉ,A作为 μ , Σ \mu,\Sigma μ,Σ的最大似然估计原型，还具有以下的性质：

    1. 无偏性： X ˉ \bar X Xˉ是 μ \mu μ的无偏估计， A / n A/n A/n不是 Σ \Sigma Σ的无偏估计，但 S = A / ( n − 1 ) S=A/(n-1) S=A/(n−1)是 Σ \Sigma Σ的无偏估计。
    2. 有效性： X ˉ \bar X Xˉ和 S S S是 μ , Σ \mu,\Sigma μ,Σ的一致最小方差无偏估计，即 X ˉ , S \bar X,S Xˉ,S是 μ , Σ \mu,\Sigma μ,Σ的有效估计量。
    3. 相合性：当 n → ∞ n\to \infty n→∞时， X ˉ , Σ ^ = A / n \bar X,\hat \Sigma=A/n Xˉ,Σ^=A/n是 μ , Σ \mu,\Sigma μ,Σ的强相合估计，即随着抽样数的增加，它们总会收敛于参数。
    4. 充分性： X ˉ , Σ ^ \bar X,\hat \Sigma Xˉ,Σ^是 μ , Σ \mu,\Sigma μ,Σ的充分统计量。

    最大似然估计满足对参数函数依然适用的性质，即对于 μ , Σ \mu,\Sigma μ,Σ的最大似然估计 μ ^ , Σ ^ \hat \mu,\hat \Sigma μ^​,Σ^，参数的函数 φ ( μ , Σ ) \varphi(\mu,\Sigma) φ(μ,Σ)的最大似然估计还是 φ ( μ ^ , Σ ^ ) \varphi(\hat \mu,\hat \Sigma) φ(μ^​,Σ^)。

    回顾总结

    1. 参数估计中，最重要的两个统计量是样本均值 X ˉ \bar X Xˉ与样本离差阵 A A A，它们与样本数据阵 X X X的关系分别是
      X ˉ = 1 n X ′ 1 n , A = X ′ X − n X ˉ X ˉ ′ = X ′ [ I n − 1 n 1 n 1 n ′ ] X . \bar X=\frac 1nX'\boldsymbol 1_n,\\ A=X'X-n\bar X\bar X'=X'\left[I_n-\frac1n\boldsymbol 1_n\boldsymbol 1_n' \right]X. Xˉ=n1​X′1n​,A=X′X−nXˉXˉ′=X′[In​−n1​1n​1n′​]X.
      还有相关的统计量如 S = A / ( n − 1 ) , S ∗ = A / n S=A/(n-1),S^*=A/n S=A/(n−1),S∗=A/n和样本相关阵 R , r i j = a i j / a i i a j j R,r_{ij}=a_{ij}/\sqrt{a_{ii}a_{jj}} R,rij​=aij​/aii​ajj​ ​。

    2. N p ( μ , Σ ) N_p(\mu,\Sigma) Np​(μ,Σ)的参数 μ , Σ \mu,\Sigma μ,Σ的最大似然估计分别是 μ ^ = X ˉ , Σ ^ = A / n \hat \mu=\bar X,\hat \Sigma=A/n μ^​=Xˉ,Σ^=A/n，一般可以由 X ˉ , A \bar X,A Xˉ,A估计出 μ , Σ \mu,\Sigma μ,Σ估计出。

    3. 关于 X ˉ , A \bar X,A Xˉ,A的性质，有 X ˉ , A \bar X,A Xˉ,A相互独立，且
      X ˉ ∼ N p ( μ , Σ / n ) , A = d ∑ α = 1 n − 1 Z α Z α ′ , Z α ∼ i . i . d . N p ( 0 , Σ ) . \bar X\sim N_p(\mu,\Sigma/n),\\ A\stackrel {\rm d}=\sum_{\alpha=1}^{n-1} Z_\alpha Z_{\alpha}',\quad Z_\alpha\stackrel {\rm i.i.d.}\sim N_p(0,\Sigma). Xˉ∼Np​(μ,Σ/n),A=dα=1∑n−1​Zα​Zα′​,Zα​∼i.i.d.Np​(0,Σ).

    4. 在无偏性方面， X ˉ \bar X Xˉ是 μ \mu μ的无偏估计， S = A / ( n − 1 ) S=A/(n-1) S=A/(n−1)是 Σ \Sigma Σ的无偏估计。

    5. 在有效性方面， X ˉ , S \bar X,S Xˉ,S是 μ , Σ \mu,\Sigma μ,Σ的最小方差无偏估计，即有效估计。

    6. 在相合性方面， X ˉ , A / n \bar X,A/n Xˉ,A/n是 μ , Σ \mu,\Sigma μ,Σ的强相合估计。

    7. 对于参数函数 φ ( μ , Σ ) \varphi(\mu,\Sigma) φ(μ,Σ)，它的最大似然估计是 φ ( X ˉ , A / n ) \varphi(\bar X,A/n) φ(Xˉ,A/n)。

    收起
    展开全文

  • 多元正态分布的参数估计

     2021-10-04 07:09:53
     在实际应用中，多元正态分布中均值向量，和协差阵。通常是未知的，需由样本来估计，而参数的估计方法很多，这里用最常见的最大似然估计法给出其估计量，并借助一元统计中学过的估计量性质指出这里给出的估计量也...

    引言

     在实际应用中，多元正态分布中均值向量，和协差阵。通常是未知的，需由样本来估计，而参数的估计方法很多，这里用最常见的最大似然估计法给出其估计量，并借助一元统计中学过的估计量性质指出这里给出的估计量也满足通常要求的性质。

    多元样本的概念及表示法

     多元分析研究的总体是多元总体，从多元总体中随机抽取 n n n个个体 X ( 1 ) , X ( 2 ) , ⋯   , X ( n ) X_{(1)},X_{(2)},\cdots,X_{(n)} X(1)​,X(2)​,⋯,X(n)​，若 X ( 1 ) , X ( 2 ) , ⋯   , X ( n ) X_{(1)},X_{(2)},\cdots,X_{(n)} X(1)​,X(2)​,⋯,X(n)​相互独立且与总体同分布，则称 X ( 1 ) , X ( 2 ) , ⋯   , X ( n ) X_{(1)},X_{(2)},\cdots,X_{(n)} X(1)​,X(2)​,⋯,X(n)​为该总体的一个多元随机样本，简称为简单样本。每个 X ( n ) = ( X a 1 , X a 2 , ⋯   , X a p ) ⊤ ( a = 1 , 2 , ⋯   , n ) X_{(n)}=(X_{a1},X_{a2},\cdots,X_{ap})^{\top}(a=1,2,\cdots,n) X(n)​=(Xa1​,Xa2​,⋯,Xap​)⊤(a=1,2,⋯,n)称为一个样品，其中， X a j X_{aj} Xaj​为第 a a a个样品对第 j j j个指标的观测值， 显然每个样品都是 p p p维向量，将 n n n个样品对 p p p项指标进行观测，将全部观测结果用一个 n × p n \times p n×p阶矩阵 X X X表示： X = [ X 11 X 12 ⋯ X 1 p X 21 X 22 ⋯ X 2 p ⋮ ⋮ ⋮ ⋮ X n 1 X n 2 ⋯ X n p ] = [ X ( 1 ) ⊤ X ( 2 ) ⊤ ⋮ X ( n ) ⊤ ] X=\left[\begin{array}{cccc}X_{11}&X_{12}&\cdots&X_{1p}\\X_{21}&X_{22}&\cdots &X_{2p}\\ \vdots&\vdots& \vdots & \vdots \\ X_{n1}&X_{n2} & \cdots& X_{np}\end{array}\right]=\left[\begin{array}{c}X_{(1)}^{\top}\\X_{(2)}^{\top}\\ \vdots \\ X_{(n)}^{\top}\end{array}\right] X=⎣⎢⎢⎢⎡​X11​X21​⋮Xn1​​X12​X22​⋮Xn2​​⋯⋯⋮⋯​X1p​X2p​⋮Xnp​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​X(1)⊤​X(2)⊤​⋮X(n)⊤​​⎦⎥⎥⎥⎤​

    多元样本的数字特征

    定义1： 设 X ( 1 ) , ⋯   , X ( n ) X_{(1)},\cdots,X_{(n)} X(1)​,⋯,X(n)​为来自 p p p元总体的样本，其中 X ( a ) = ( X a 1 , ⋯   , X a p ) ⊤ , a = 1 , 2 , ⋯   , n , X_{(a)}=(X_{a1},\cdots,X_{ap})^{\top},a=1,2,\cdots,n, X(a)​=(Xa1​,⋯,Xap​)⊤,a=1,2,⋯,n,
    （1）样本均值向量定义为： X ˉ = Δ 1 n ∑ a = 1 X ( a ) = ( X ˉ 1 , X ˉ 2 , ⋯   , X ˉ p ) ⊤ \bar{X}\stackrel{\Delta}{=}\frac{1}{n}\sum\limits_{a=1}X_{(a)}=(\bar{X}_1,\bar{X}_2,\cdots,\bar{X}_p)^{\top} Xˉ=Δn1​a=1∑​X(a)​=(Xˉ1​,Xˉ2​,⋯,Xˉp​)⊤ ∵ 1 n ∑ a = 1 n X ( a ) = 1 n [ [ X 11 X 12 ⋮ X 1 p ] + [ X 21 X 22 ⋮ X 2 p ] + ⋯ + [ X n 1 X n 2 ⋮ X n p ] ] = 1 n [ X 11 + X 21 + ⋯ + X n 1 X 12 + X 22 + ⋯ + X n 2 ⋮ X 1 p + X 2 p + ⋯ + X n p ] = [ X ˉ 1 X ˉ 2 ⋮ X ˉ p ] \begin{aligned}\because \frac{1}{n}\sum\limits_{a=1}^nX_{(a)}&=\frac{1}{n}\left[\left[\begin{array}{c}X_{11}\\X_{12}\\\vdots\\X_{1p}\end{array}\right]+\left[\begin{array}{c}X_{21}\\X_{22}\\\vdots\\X_{2p}\end{array}\right]+\cdots+\left[\begin{array}{c}X_{n1}\\X_{n2}\\\vdots\\X_{np}\end{array}\right]\right]\\&=\frac{1}{n}\left[\begin{array}{c}X_{11}+X_{21}+\cdots+X_{n1}\\X_{12}+X_{22}+\cdots+X_{n2}\\\vdots\\X_{1p}+X_{2p}+\cdots+X_{np} \end{array}\right]\\&=\left[\begin{array}{c}\bar{X}_1\\\bar{X}_2\\ \vdots\\\bar{X}_p\end{array}\right]\end{aligned} ∵n1​a=1∑n​X(a)​​=n1​⎣⎢⎢⎢⎡​⎣⎢⎢⎢⎡​X11​X12​⋮X1p​​⎦⎥⎥⎥⎤​+⎣⎢⎢⎢⎡​X21​X22​⋮X2p​​⎦⎥⎥⎥⎤​+⋯+⎣⎢⎢⎢⎡​Xn1​Xn2​⋮Xnp​​⎦⎥⎥⎥⎤​⎦⎥⎥⎥⎤​=n1​⎣⎢⎢⎢⎡​X11​+X21​+⋯+Xn1​X12​+X22​+⋯+Xn2​⋮X1p​+X2p​+⋯+Xnp​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​Xˉ1​Xˉ2​⋮Xˉp​​⎦⎥⎥⎥⎤​​
    （2）样本离差阵定义为： S p × p = Δ ∑ a = 1 n ( X ( a ) − X ˉ ) ( X ( a ) − X ˉ ) ⊤ = ( S i j ) p × p S_{p \times p}\stackrel{\Delta}{=}\sum\limits_{a=1}^n(X_{(a)}-\bar{X})(X_{(a)}-\bar{X})^{\top}=(S_{ij})_{p \times p} Sp×p​=Δa=1∑n​(X(a)​−Xˉ)(X(a)​−Xˉ)⊤=(Sij​)p×p​ ∵ ∑ a = 1 ( X ( a ) − X ˉ ) ( X ( a ) − X ˉ ) ⊤ = ∑ a = 1 n [ [ X a 1 − X ˉ 1 X a 2 − X ˉ 2 ⋮ X a p − X ˉ p ] ( X a 1 − X ˉ 1 , X a 2 − X ˉ 2 , ⋯   , X a p − X ˉ ) ] = ∑ [ ( X a 1 − X ˉ 1 ) 2 ( X a 1 − X ˉ 1 ) ( X a 2 − X ˉ 2 ) ⋯ ( X a 1 − X ˉ 1 ) ( X a p − X ˉ p ) ( X a 2 − X ˉ 2 ) ( X a 1 − X ˉ 1 ) ( X a 2 − X ˉ 2 ) 2 ⋯ ( X a 2 − X ˉ 2 ) ( X a p − X ˉ p ) ⋮ ⋮ ⋮ ( X a p − X ˉ p ) ( X a 1 − X ˉ 1 ) ( X a p − X ˉ p ) ( X a 2 − X ˉ 2 ) ⋯ ( X a p − X ˉ p ) 2 ] = [ S 11 S 12 ⋯ S 1 p S 21 S 22 ⋯ S 2 p ⋮ ⋮ ⋮ S p 1 S p 2 ⋯ S p p ] = ( S i j ) p × p \begin{aligned}\because &\sum\limits_{a=1}(X_{(a)}-\bar{X})(X_{(a)}-\bar{X})^{\top}\\&=\sum\limits_{a=1}^n\left[\left[\begin{array}{c}X_{a1}-\bar{X}_1\\X_{a2}-\bar{X}_2\\\vdots\\ X_{ap}-\bar{X}_p\end{array}\right](X_{a1}-\bar{X}_1,X_{a2}-\bar{X}_2,\cdots,X_{ap}-\bar{X})\right]\\&=\sum\left[\begin{array}{cccc}(X_{a1}-\bar{X}_1)^2 & (X_{a1}-\bar{X}_1)(X_{a2}-\bar{X}_2)& \cdots & (X_{a1}-\bar{X}_1)(X_{ap}-\bar{X}_p)\\(X_{a2}-\bar{X}_2)(X_{a1}-\bar{X}_1)&(X_{a2}-\bar{X}_2)^2 &\cdots& (X_{a2}-\bar{X}_2)(X_{ap}-\bar{X}_p) \\ \vdots&\vdots&&\vdots\\ (X_{ap}-\bar{X}_p)(X_{a1}-\bar{X}_1) &(X_{ap}-\bar{X}_p)(X_{a2}-\bar{X}_2)&\cdots &(X_{ap}-\bar{X}_p)^2 \end{array}\right]\\&=\left[\begin{array}{cccc}S_{11}&S_{12}&\cdots&S_{1p}\\S_{21}&S_{22}&\cdots & S_{2p}\\\vdots &\vdots && \vdots\\S_{p1}&S_{p2}&\cdots&S_{pp}\end{array}\right]=(S_{ij})_{p \times p}\end{aligned} ∵​a=1∑​(X(a)​−Xˉ)(X(a)​−Xˉ)⊤=a=1∑n​⎣⎢⎢⎢⎡​⎣⎢⎢⎢⎡​Xa1​−Xˉ1​Xa2​−Xˉ2​⋮Xap​−Xˉp​​⎦⎥⎥⎥⎤​(Xa1​−Xˉ1​,Xa2​−Xˉ2​,⋯,Xap​−Xˉ)⎦⎥⎥⎥⎤​=∑⎣⎢⎢⎢⎡​(Xa1​−Xˉ1​)2(Xa2​−Xˉ2​)(Xa1​−Xˉ1​)⋮(Xap​−Xˉp​)(Xa1​−Xˉ1​)​(Xa1​−Xˉ1​)(Xa2​−Xˉ2​)(Xa2​−Xˉ2​)2⋮(Xap​−Xˉp​)(Xa2​−Xˉ2​)​⋯⋯⋯​(Xa1​−Xˉ1​)(Xap​−Xˉp​)(Xa2​−Xˉ2​)(Xap​−Xˉp​)⋮(Xap​−Xˉp​)2​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​S11​S21​⋮Sp1​​S12​S22​⋮Sp2​​⋯⋯⋯​S1p​S2p​⋮Spp​​⎦⎥⎥⎥⎤​=(Sij​)p×p​​
    （3）样本协差阵定义为： V p × p = Δ 1 n S = 1 n ∑ a = 1 n ( X ( a ) − X ˉ ) ( X ( a ) − X ˉ ) ⊤ = ( v i j ) p × p V_{p \times p}\stackrel{\Delta}{=}\frac{1}{n}S=\frac{1}{n}\sum\limits_{a=1}^n(X_{(a)}-\bar{X})(X_{(a)}-\bar{X})^{\top}=(v_{ij})_{p \times p} Vp×p​=Δn1​S=n1​a=1∑n​(X(a)​−Xˉ)(X(a)​−Xˉ)⊤=(vij​)p×p​ ∵ 1 n S = 1 n ∑ a = 1 n ( X ( a ) − X ˉ ) ( X ( a ) − X ˉ ) ⊤ = [ 1 n ∑ a = 1 n ( X a i − X ˉ i ) ( X a j − X j ) ] p × p = [ v i j ] p × p \because \begin{aligned}\frac{1}{n}S&=\frac{1}{n}\sum\limits_{a=1}^{n}(X_{(a)}-\bar{X})(X_{(a)}-\bar{X})^{\top}\\&=\left[\frac{1}{n}\sum\limits_{a=1}^n(X_{ai}-\bar{X}_i)(X_{aj}-X_{j})\right]_{p \times p}\\&=[v_{ij}]_{p \times p}\end{aligned} ∵n1​S​=n1​a=1∑n​(X(a)​−Xˉ)(X(a)​−Xˉ)⊤=[n1​a=1∑n​(Xai​−Xˉi​)(Xaj​−Xj​)]p×p​=[vij​]p×p​​
    （4）样本相关阵定义为： R p × p = Δ ( r i j ) p × p R_{p \times p}\stackrel{\Delta}{=}(r_{ij})_{p \times p} Rp×p​=Δ(rij​)p×p​，其中， r i j = v i j v i i v j j = s i j s i i s j j r_{ij}=\frac{v_{ij}}{\sqrt{v_{ii}}\sqrt{v_{jj}}}=\frac{s_{ij}}{\sqrt{s_{ii}}\sqrt{s_{jj}}} rij​=vii​ ​vjj​ ​vij​​=sii​ ​sjj​ ​sij​​样本均值向量和离差阵也可用 X X X直接表示如下： X ˉ p × 1 = 1 n X ⊤ 1 n \bar{X}_{p \times 1}=\frac{1}{n}X^{\top}1_n Xˉp×1​=n1​X⊤1n​其中， 1 n = ( 1 , 1 , ⋯   , 1 ) ⊤ 1_n=(1,1,\cdots,1)^{\top} 1n​=(1,1,⋯,1)⊤ X ˉ = 1 n X ⊤ 1 n = 1 n [ X 11 X 21 ⋯ X n 1 X 12 X 22 ⋯ X n 2 ⋮ ⋮ ⋮ X 1 p X 2 p ⋯ X n p ] [ 1 1 ⋮ 1 ] = 1 n [ X 11 + X 21 + ⋯ + X n 1 X 12 + X 22 + ⋯ + X n 2 ⋮ X 1 p + X 2 p + ⋯ + X n p ] = [ X ˉ 1 X ˉ 2 ⋮ X ˉ p ] \begin{aligned}\bar{X}&=\frac{1}{n}X^{\top}1_n\\&=\frac{1}{n}\left[\begin{array}{cccc}X_{11}&X_{21}&\cdots&X_{n1}\\X_{12}&X_{22}&\cdots&X_{n2}\\\vdots&\vdots&&\vdots\\X_{1p}&X_{2p}&\cdots&X_{np}\end{array}\right]\left[\begin{array}{c}1\\1\\\vdots\\1\end{array}\right]\\&=\frac{1}{n}\left[\begin{array}{c}X_{11}+X_{21}+\cdots+X_{n1}\\X_{12}+X_{22}+\cdots+X_{n2}\\\vdots\\ X_{1p}+X_{2p}+\cdots+X_{np}\end{array}\right]=\left[\begin{array}{c}\bar{X}_1\\\bar{X}_2\\\vdots\\\bar{X}_p\end{array}\right]\end{aligned} Xˉ​=n1​X⊤1n​=n1​⎣⎢⎢⎢⎡​X11​X12​⋮X1p​​X21​X22​⋮X2p​​⋯⋯⋯​Xn1​Xn2​⋮Xnp​​⎦⎥⎥⎥⎤​⎣⎢⎢⎢⎡​11⋮1​⎦⎥⎥⎥⎤​=n1​⎣⎢⎢⎢⎡​X11​+X21​+⋯+Xn1​X12​+X22​+⋯+Xn2​⋮X1p​+X2p​+⋯+Xnp​​⎦⎥⎥⎥⎤​=⎣⎢⎢⎢⎡​Xˉ1​Xˉ2​⋮Xˉp​​⎦⎥⎥⎥⎤​​ S = X ⊤ ( I − 1 n 1 n 1 n ⊤ ) X S=X^{\top}(I-\frac{1}{n}1_n1_n^{\top})X S=X⊤(I−n1​1n​1n⊤​)X其中， I n = [ 1 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 1 ] I_n=\left[\begin{array}{ccc}1&\cdots&0\\ \vdots & \ddots & \vdots \\ 0 & \cdots &1\end{array}\right] In​=⎣⎢⎡​1⋮0​⋯⋱⋯​0⋮1​⎦⎥⎤​ ∵ S = ∑ a = 1 n ( X ( a ) − X ˉ ) ( X ( a ) − X ˉ ) ⊤ = X ⊤ X − n X ˉ X ˉ ⊤ = X ⊤ X − 1 n X ⊤ 1 n 1 n ⊤ X = X ⊤ ( I n − 1 n 1 n 1 n ⊤ ) X \begin{aligned}\because S&=\sum\limits_{a=1}^n(X_{(a)}-\bar{X})(X_{(a)}-\bar{X})^{\top}\\&=X^{\top}X-n \bar{X}\bar{X}^{\top}\\&=X^{\top}X-\frac{1}{n}X^{\top}1_n1_n^{\top}X\\&=X^{\top}(I_n-\frac{1}{n}1_n1_n^{\top})X\end{aligned} ∵S​=a=1∑n​(X(a)​−Xˉ)(X(a)​−Xˉ)⊤=X⊤X−nXˉXˉ⊤=X⊤X−n1​X⊤1n​1n⊤​X=X⊤(In​−n1​1n​1n⊤​)X​

    μ \mu μ和 Σ \Sigma Σ的最大似然估计及基本性质

     通过样本来估计总体的参数叫做参数估计，参数估计的原则和方法是很多的，这里用常见的且具有很多优良性质的最大似然法给出 μ \mu μ和 Σ \Sigma Σ的估计量。
     设 X ( 1 ) , X ( 2 ) , ⋯   , X ( n ) X_{(1)},X_{(2)},\cdots,X_{(n)} X(1)​,X(2)​,⋯,X(n)​是来自正态总体 N p ( μ , Σ ) N_p(\mu,\Sigma) Np​(μ,Σ)容量为 n n n的样本，每个样品： X ( a ) = ( X a 1 , X a 2 , ⋯   , X a p ) ⊤ X_(a)=(X_{a1},X_{a2},\cdots,X_{ap})^{\top} X(​a)=(Xa1​,Xa2​,⋯,Xap​)⊤， a = 1 , ⋯   , n a=1,\cdots,n a=1,⋯,n，样本阵 X X X为： X = [ X 11 X 12 ⋯ X 1 p X 21 X 22 ⋯ X 2 p ⋮ ⋮ ⋱ ⋮ X n 1 X n 2 ⋯ X n p ] X=\left[\begin{array}{cccc}X_{11}&X_{12}&\cdots&X_{1p}\\X_{21}&X_{22}&\cdots &X_{2p}\\\vdots&\vdots&\ddots&\vdots \\ X_{n1}&X_{n2}&\cdots &X_{np}\end{array}\right] X=⎣⎢⎢⎢⎡​X11​X21​⋮Xn1​​X12​X22​⋮Xn2​​⋯⋯⋱⋯​X1p​X2p​⋮Xnp​​⎦⎥⎥⎥⎤​则用最大似然估计求出 μ \mu μ和 Σ \Sigma Σ的估计量分别为： μ ^ = X ˉ , Σ ^ = 1 n S \hat{\mu}=\bar{X},\quad \hat{\Sigma}=\frac{1}{n}S μ^​=Xˉ,Σ^=n1​S μ \mu μ和 Σ \Sigma Σ的估计量有如下基本性质：
    （1） E ( X ˉ ) = μ E(\bar{X})=\mu E(Xˉ)=μ，即 X ˉ \bar{X} Xˉ是 μ \mu μ的无偏估计； E ( 1 n S ) = n − 1 n Σ E(\frac{1}{n}S)=\frac{n-1}{n}\Sigma E(n1​S)=nn−1​Σ，即 1 n S \frac{1}{n}S n1​S不是 Σ \Sigma Σ的无偏估计；而 E ( 1 n − 1 S ) = Σ E(\frac{1}{n-1}S)=\Sigma E(n−11​S)=Σ，即 1 n − 1 S \frac{1}{n-1}S n−11​S是 Σ \Sigma Σ的无偏估计；
    （2） X ˉ , 1 n − 1 S \bar{X},\frac{1}{n-1}S Xˉ,n−11​S分别是 μ , Σ \mu,\Sigma μ,Σ的有效估计。
    （3） X ˉ , 1 n S \bar{X},\frac{1}{n}S Xˉ,n1​S分别是 μ , Σ \mu,\Sigma μ,Σ的一致估计（相合估计）。

    定理1： 设 X ˉ \bar{X} Xˉ和 S S S分别正态总体 N p ( μ , Σ ) N_p(\mu,\Sigma) Np​(μ,Σ)的样本均值向量和离差阵，则
    （1） X ˉ ∼ N p ( μ , 1 n Σ ) \bar{X}\sim N_p(\mu,\frac{1}{n}\Sigma) Xˉ∼Np​(μ,n1​Σ)；
    （2）离差阵 S S S可以写为： S = ∑ a = 1 n − 1 Z a Z a ⊤ , S=\sum\limits_{a=1}^{n-1}Z_aZ_a^{\top}, S=a=1∑n−1​Za​Za⊤​,其中 Z 1 , ⋯   , Z n − 1 Z_1,\cdots,Z_{n-1} Z1​,⋯,Zn−1​独立同分布与 N p ( 0 , Σ ) N_p(0,\Sigma) Np​(0,Σ)；
    （3） X ˉ \bar{X} Xˉ与 S S S相互独立；
    （4） S S S为正定阵的充要条件是 n > p n>p n>p。

    收起
    展开全文

  • 3-25 浅谈多元正态分布的参数估计

     2021-03-25 14:25:27
    再求 ，使得似然函数最大]： 最大似然估计的性质： 由于充分利用了中的信息，极大似然估计 和 有着非常好的性质(估计的精确度很高），同时具有： 无偏性（需要修正），有效性，相合性，渐进正态性，且是 及 的...

    数字特征天团四人组：

    均值向量\small \overline{X} , 离差阵 A ，协方差阵 S ,相关阵 R 

    其中离差阵除以自由度就是协方差阵，相关阵的元素则由相关系数组成。

     

    \small \mu 及 \tiny \sum 的极大似然估计

    似然函数：将数据阵 X 进行拉直运算后形成的np维向量的联合密度函数被称为样本\small x_{(i)}的似然函数(关于\small \mu 及 \tiny \sum ）

    极大似然估计：认为样本来自使样本出现概率最大的总体！

    了解思想后证明不难[因为存在两个变量，所以需要先固定\small \widehat{\mu}，再求\tiny \widehat{\sum} ，使得似然函数最大]：

    最大似然估计的性质：

    由于充分利用了\small X_{(i)}中的信息，极大似然估计 \small \widehat{\mu} 和 \tiny \widehat{\sum} 有着非常好的性质(估计的精确度很高），同时具有：

    无偏性（\tiny \widehat{\sum}需要修正），有效性，相合性，渐进正态性，且是\small \mu 及 \tiny \sum的充分统计量。（数理统计的总结上会着重聊到，回头补链接= =）

     

     

     

     

     

    收起
    展开全文

  • 多元高斯分布

     2018-07-27 21:40:20
    多元高斯分布函数的python程序，供大家学习整理和使用
    收起

  • 多元正态分布的参数估计PPT教案.pptx

     2021-10-04 06:04:42
    多元正态分布的参数估计PPT教案.pptx
    收起

  • 多元正态分布的参数估计PPT学习教案.pptx

     2021-10-05 06:24:28
    多元正态分布的参数估计PPT学习教案.pptx
    收起

  • 多元正态分布的参数估计与检验学习教案.pptx

     2021-11-20 08:33:18
    多元正态分布的参数估计与检验学习教案.pptx
    收起

  • 多元正态分布的参数估计与检验PPT课件.pptx

     2021-10-11 21:00:25
    多元正态分布的参数估计与检验PPT课件.pptx
    收起

  • 多元正态分布的参数估计与检验PPT学习教案.pptx

     2021-10-02 21:31:00
    多元正态分布的参数估计与检验PPT学习教案.pptx
    收起

  • 多元正态分布最大似然估计

     千次阅读 2020-09-15 12:48:15
    多元正态分布的概率密度函数 N维随机向量 ...多元正态分布的最大似然估计 我们对均值求偏导 针对上面的矩阵求导，给出如下证明： 下面两个是常用的两个公式： 第一个公式证明，和上面的类似。 第二个小编
    收起

  • 多元统计分析第01讲--多元正态分布及参数估计（随机向量，多元正态分布定义，条件分布和独立性）

     2021-01-17 19:03:48
    第二章、多元正态分布及参数估计这一讲主要是给出概率论中若干概念向高维的推广2.1随机向量一、随机向量的联合分布、边缘分布和条件分布1、多元数据 维随机向量： ，其中每个 都是随机变量随机矩阵： ，其中每个 都...
    收起

  • 机器学习系列（二）多元正态分布

     千次阅读 2021-03-07 11:11:14
    (Univariate normal distribution)， ，则其概率密度函数为：整个分布可以仅用均值及方差来刻画如果变量之间不相关，则它们相互独立经典统计检验通常基于正态分布假设正态分布可以模拟大量自然现象多元正态分布多元...
    收起

  • 多元正态分布的极大似然估计

     万次阅读 多人点赞 2018-06-20 13:53:12
    多元正态分布的极大似然估计 1. 一元正态分布的密度函数 一元正态分布的密度函数表示为： f(x)=1(2π)−−−−√σe−(x−μ)22σ2f(x)=1(2π)σe−(x−μ)22σ2f(x) = \frac{1}{\sqrt {(2 \pi)} \sigma} e^{...
    收起

  • 【多元统计分析】02.多元正态分布

     千次阅读 2020-10-24 11:30:38
    本文讨论了多元正态分布的定义，重点讨论多元正态分布的独立性、回归与最佳预测等问题。
    收起

  • MATLAB实现多元正态Copula分布

     2022-04-18 15:15:33
    此篇文章主要想尝试一下用MATLAB现成的函数实现多元正态Copula。感受一下多维Copula是个什么感觉。 1 与正态Copula有关的MATLAB函数 1.1 copulafit函数 copulafit函数用来根据样本观测数据估计Copula函数中的未知...
    收起
    matlab

  • 正态分布时的贝叶斯估计

     千次阅读 2021-09-20 17:23:55
    正态分布时的贝叶斯估计 前导知识：【贝叶斯估计】 下面以最简单的一维正态分布模型为例来说明贝叶斯估计的应用。假设模型的均值μ\muμ是待估计的参数，方差为σ2\sigma^2σ2为已知，分布密度写为： ρ(x∣μ)=12π...
    收起
    加权平均 密度函数

  • 03 ，seaborn 颜色 ： 正态分布图，多元正态分布，核密度估计图

     2020-07-08 17:16:45
    1 ，多元正态分布 ： 数据 代码 ： x, y = np.random.multivariate_normal([0, 0], [[1, -.5], [-.5, 1]], size=300).T 2 ，核密度估计图 ：sns.kdeplot 含义 ： 根据已有数据，画图连线，推测数据走向 3 ，...
    收起
    python

  • 【python数据分析】正态分布、正态性检验与相关性分析

     千次阅读 2021-02-22 22:08:25
    正态分布、正态性检验与相关性分析1 正态分布2 正态性检验2.1 直方图初判2.2 QQ图2.3 K-S检验2.3.1推导过程2.3.2 直接一行代码调用3 相关性分析3.1 图示初判3.2 Pearson相关系数3.2.1 计算推导3.2.2 代码一步到位3.3...
    收起
    python 算法 相关性分析

  • 混合正态分布参数极大似然估计的EM算法.pdf

     2011-07-15 10:40:24
    混合正态分布参数极大似然估计的EM算法.pdf
    收起

  • 两个多元正态分布的KL散度、巴氏距离和W距离

     千次阅读 2021-07-24 00:07:40
    ©PaperWeekly 原创 ·作者|苏剑林单位|追一科技研究方向|NLP、神经网络正态分布是最常见的连续型概率分布之一。它是给定均值和协方差后的最大熵分布（参考《“熵”不...
    收起
    机器学习 人工智能

  • 正态分布下的最大似然估计

     千次阅读 2021-09-19 15:11:18
    前导知识：【最大似然参数估计的求解】 本文仅以单变量正态分布情况下估计其均值和方差为例来说明最大似然估计的用法。 单变量正态分布的形式为： ρ(x∣θ)=12πσe−12(x−μσ)2(1) \rho(x|\theta)=\frac{1}{\...
    收起

  • 多元正太分布下贝叶斯估计法的参数估计

     2019-10-18 11:03:17
    这里给出多元正太分布下的贝叶斯估计的推导过程。（本人比较懒，直接上图趴。。。） 至此推导出了均值与协方差矩阵测估计量；可以与一元正太分布的结果做对比， 发现有很高的相似性！数学之美就在于此嘿嘿嘿~ ...
    收起
    模式识别

  • mvncond2(Mean, Sigma, ind, values)：多元正态分布条件均值和协方差矩阵的向量估计。-matlab开发

     2021-05-29 20:55:28
    此函数提供多元正态分布条件期望和协方差矩阵的矢量化估计。 均值是一个矩阵，其中行表示期望向量。 Sigma 是协方差矩阵。 Ind 是第一个无条件参数的索引。 值是条件值的矩阵，其中的行对应于平均行。
    收起
    matlab

  • 【ML学习笔记】17：多元正态分布下极大似然估计最小错误率贝叶斯决策

     千次阅读 2018-01-03 10:36:18
    简述多元正态分布下的最小错误率贝叶斯如果特征的值向量服从d元正态分布，即其概率密度函数为： 即其分布可以由均值向量和对称的协方差矩阵 唯一确定。如果认为样本的特征向量在类内服从多元正态分布： 即...
    收起
    机器学习 正态分布 最小错误率贝叶斯

  • 概率笔记12——多维正态分布的最大似然估计

     万次阅读 多人点赞 2019-08-19 19:33:18
    我们在前面的章节中见识过二维正态分布，(X,Y)服从参数为μ1, μ2, σ1, σ2, ρ的二维正态分布，记作(X, Y)~N(μ1, μ2, σ1, σ2, ρ)，它的密度函数： 　其中μ1是第1维度的均值，σ12是第1维度的方差，ρ是将...
    收起

  • （多元）偏正态分布、正态分布、对数正态分布的随机数的产生（R语言）

     千次阅读 2021-07-24 17:13:15
    目录0引言1、函数名2、示例2.1正态分布随机数2.2偏正态分布2.3对数正态分布写在最后的话 0引言 最近在看偏正态分布相关的东西，偏正态分布的定义形式还是挺多样的，在偏态分布及其数字特征（R语言可视化）中我介绍的...
    收起
    R语言

  • 《机器学习笔记（三）：多元线性回归与正态分布最大似然估计》

     千次阅读 2019-01-29 17:37:40
    回归问题普遍讨论的是多元线性回归，考虑多个特征可以得到更精确的模型，这其中涉及中心极限定理，正态分布，概率密度函数和最大似然估计。 (一)背景——多元线性回归 1. 本质上就是算法（公式）变换为了多元一次...
    收起
    机器学习

  • 拓端tecdat|R语言混合正态分布极大似然估计和EM算法

     千次阅读 2020-04-27 22:30:14
    为了在统计过程中发现更多有趣的结果，我们将解决极大似然估计没有简单分析表达式的情况。举例来说，如果我们混合了各种分布，
    收起
    R语言 EM算法

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 5,477
精华内容 2,190
热门标签