精华内容
下载资源
问答
  • 贝叶斯判别算法matlab的实现,详细地介绍了贝叶斯算法
  • 贝叶斯判别分析,用于降维和分类,包含多个程序,亲测可以有实例
  • 贝叶斯判别matlab

    2020-05-22 09:00:24
    matlab实现的贝叶斯距离判别法,程序设计使用的
  • 贝叶斯判别

    千次阅读 2019-10-18 17:13:13
    贝叶斯判别 根据概率判别规则,有: 若P(ω1 | x) > P(ω2 | x),则x∈ω1 若P(ω1 | x) < P(ω2 | x),则x∈ω2 由贝叶斯定理,后验概率P(ωi | x)可由类别ωi的先验概率P(ωi)和x的条件概率密度p(x | ω...

    贝叶斯判别

    根据概率判别规则,有:

    P(ω1 | x) > P(ω2 | x)xω1 

    P(ω1 | x) < P(ω2 | x)xω2

    由贝斯定理,后验概率P(ωi | x)可由类别ωi的先验概率P(ωi)和x的条件概率密度p(x | ωi)来计算,即:

    这里p(x | ωi)也称为似然函数。将该式代入上述判别式,有:

    p(x | ω1)P(ω1) > p(x | ω2)P(ω2)

    xω1

    p(x | ω1)P(ω1) < p(x | ω2)P(ω2)

    xω2

    xω1

    ,则xω2

    其中,l12称为似然比,P(ω2)/P(ω1)=θ21称为似然比的判决阈值,此判别称为贝叶斯判别。

     

     

    展开全文
  • 贝叶斯判别程序

    2013-05-17 13:42:20
    数学判别 discriminiant.bayes (TrnX1, TrnX2, rate=1, TstX = NULL, var.equal = FALSE){ if (is.null(TstX) == TRUE) TstX(TrnX1,TrnX2) if (is.vector(TstX) == TRUE) TstX(as.matrix(TstX)) else if (is....
  • 【多元统计分析】14.贝叶斯判别

    千次阅读 2020-11-11 14:16:48
    贝叶斯判别法是基于平均错判损失最小的前提下,定义的判别法。

    十四、贝叶斯判别法

    1.贝叶斯判别的定义

    贝叶斯判别的定义,是找到一个错判平均损失最小的判别准则,这句话虽然简单,但还有一些概念需要解析,接下来我们假设有 k k k个总体 G 1 , ⋯   , G k G_1,\cdots,G_k G1,,Gk

    首先,错判损失指的是将属于某类的实体错判为其他类,在实际生活中会导致的损失。比如考虑今天会不会下雨的判别,这决定了你出门是否带雨伞,如果今天实际上出太阳,但你判断今天会下雨,这将导致你需要多承受一把雨伞的重量,带来了一定的损失;但如果今天实际上下雨,但你判断今天会出太阳,这将导致你承担被雨淋的痛苦或者等伞的无聊,也带来损失。这两种损失给你造成的体验是否一样?显然下雨错判为晴天的损失更大一些。而在实际的问题中,不同情况的错判损失也很可能不同,因此有必要加以区分。

    使用判别法 D D D将第 i i i类的样本错判为第 j j j类,错判损失记作 L ( j ∣ i ; D ) = L ( j ∣ i ) L(j|i;D)=L(j|i) L(ji;D)=L(ji),一般错判损失可比较而不可量化,但在应用贝叶斯判别法的情况下必须量化。量化方式可以是经验赋值,对所有错判损失给一个大致的判断;而如果不同类别的错判损失大致相同,则定义 L ( j ∣ i ) = 1 − δ i j L(j|i)=1-\delta_{ij} L(ji)=1δij

    其次,既然是错判平均损失,就存在一种平均准则。使用算术平均是否合适呢?事实上是不合适的。首先,每种错判的发生可能不一样,假设实体来自 i i i类,在观测前使用某判别准则将其判断到 j j j类的概率是固定的,即 P ( j ∣ i ) P(j|i) P(ji),这样,如果实体来自 i i i类,则此时的错判损失是
    r i ( D ) = ∑ j = 1 k P ( j ∣ i ) L ( j ∣ i ) r_i(D)=\sum_{j=1}^k P(j|i)L(j|i) ri(D)=j=1kP(ji)L(ji)
    但也不能够直接将属于每一类的错判损失求算术平均,因为实体来自每一类的概率本身就不同,这称为先验概率

    先验概率代表了出现类别的概率分布,这是在没有任何样本信息时能做出的关于类的直接判断。假设来自第 i i i类的先验概率是 q i q_i qi,那么此时的错判平均损失,实际上是一种关于先验概率的加权平均。现在,我们可以定义判别法 D D D的错判平均损失为
    g ( D ) = ∑ i = 1 k q i ∑ j = 1 k P ( j ∣ i ) L ( j ∣ i ) = d e f ∑ i = 1 k q i r i ( D ) . g(D)=\sum_{i=1}^kq_i\sum_{j=1}^kP(j|i)L(j|i)\stackrel {\rm def}= \sum_{i=1}^kq_ir_i(D). g(D)=i=1kqij=1kP(ji)L(ji)=defi=1kqiri(D).
    这样,贝叶斯判别准则就可以被视为这样的最优化问题:找到一个 D ∗ D^* D,使得 g ( D ∗ ) = min ⁡ D g ( D ) g(D^*)=\min_Dg(D) g(D)=minDg(D)

    2.贝叶斯判别的解

    如何找到使得错判平均最小的判别准则,就是贝叶斯判别的求解问题。现在,我们假设 k k k m m m维总体 G 1 , ⋯   , G k G_1,\cdots,G_k G1,,Gk的先验概率分别为 q 1 , ⋯   , q k q_1,\cdots,q_k q1,,qk,每个 G i G_i Gi的联合密度函数为 f i ( X ) f_i(X) fi(X),错判损失为 L ( j ∣ i ) L(j|i) L(ji)。任何一种判别法 D D D,都将样本空间 R m \R^m Rm划分成 k k k个(连通与否的)区域 { D 1 , ⋯   , D k } \{D_1,\cdots,D_k\} {D1,,Dk},这里 D j D_j Dj表示样本落在被判别到 j j j类的区域。

    据此,我们可以先给出错判概率为 P ( j ∣ i ) P(j|i) P(ji),它表示样品 X X X本身来自密度函数 f i ( X ) f_i(X) fi(X),但落在区域 D j D_j Dj内:
    P ( j ∣ i ; D ) = ∫ D j f i ( X ) d X . P(j|i;D)=\int_{D_j}f_i(X){\rm d}X. P(ji;D)=Djfi(X)dX.
    所以此时的错判平均损失是
    g ( D ) = ∑ i = 1 k q i ∑ j = 1 k L ( j ∣ i ) ∫ D j f i ( X ) d X = ∑ j = 1 k ∫ D j ∑ i = 1 k q i L ( j ∣ i ) f i ( X ) d X = d ∑ j = 1 k ∫ D j h j ( X ) d X . \begin{aligned} g(D)=& \sum_{i=1}^k q_i\sum_{j=1}^k L(j|i)\int_{D_j}f_i(X){\rm d}X \\ =&\sum_{j=1}^k \int_{D_j}\sum_{i=1}^k q_iL(j|i)f_i(X){\rm d}X \\ \stackrel {\rm d}=&\sum_{j=1}^k\int_{D_j} h_j(X){\rm d}X. \end{aligned} g(D)===di=1kqij=1kL(ji)Djfi(X)dXj=1kDji=1kqiL(ji)fi(X)dXj=1kDjhj(X)dX.
    这里定义
    h j ( X ) = d e f ∑ i = 1 k q i L ( j ∣ i ) f i ( X ) , h_j(X)\stackrel {\rm def}= \sum_{i=1}^k q_iL(j|i)f_i(X), hj(X)=defi=1kqiL(ji)fi(X),
    它表示把样品 X X X归到 G j G_j Gj类的平均损失,注意到 h j ( X ) h_j(X) hj(X) D D D无关,对 h j ( X ) h_j(X) hj(X)求和,就得到了错判平均损失。对于贝叶斯判别的解 D ∗ D^* D,要使得 g ( D ∗ ) g(D^*) g(D)是所有 D D D中最小的,所以
    g ( D ∗ ) − g ( D ) = ∑ i = 1 k ∫ D i ∗ h i ( X ) d X − ∑ j = 1 k ∫ D j h j ( X ) d X = ∑ i = 1 k ∑ j = 1 k ∫ D i ∗ ∩ D j [ h i ( X ) − h j ( X ) ] d X ≤ 0. \begin{aligned} & g(D^*)-g(D) \\ =&\sum_{i=1}^k \int_{D_i^*}h_i(X){\rm d}X-\sum_{j=1}^k \int_{D_j}h_j(X){\rm d}X \\ =&\sum_{i=1}^k \sum_{j=1}^k \int_{D_i^*\cap D_j}[h_i(X)-h_j(X)]{\rm d}X \\ \le&0. \end{aligned} ==g(D)g(D)i=1kDihi(X)dXj=1kDjhj(X)dXi=1kj=1kDiDj[hi(X)hj(X)]dX0.
    由此,我们能够得到贝叶斯判别的解是:在所有 h i ( X ) h_i(X) hi(X)中,如果 h j ( X ) h_j(X) hj(X)最小,就将 X X X判别为第 j j j类,在这个判别法的条件下,能够让 g ( D ∗ ) − g ( D ) ≤ 0 g(D^*)-g(D)\le 0 g(D)g(D)0恒成立。

    特别当我们指定错判损失都相等的情况下,如果 h i ( X ) < h j ( X ) h_i(X)<h_j(X) hi(X)<hj(X) h i ( X ) − h j ( X ) < 0 h_i(X)-h_j(X)<0 hi(X)hj(X)<0,那么就有 q j f j ( X ) < q i f i ( X ) q_jf_j(X)<q_if_i(X) qjfj(X)<qifi(X),所以如果在 i , j i,j i,j类中将 X X X判定为 i i i类,就应该让 q i f i ( X ) q_if_i(X) qifi(X)更大,所以错判损失都相等的情况下,贝叶斯判别的解是:在所有 q i f i ( X ) q_if_i(X) qifi(X)中,如果 q j f j ( X ) q_jf_j(X) qjfj(X)最大,就将 X X X判别为第 j j j。在此基础上,如果先验概率都相等,则在所有 f i ( X ) f_i(X) fi(X)中,如果 f j ( X ) f_j(X) fj(X)最大,就将 X X X判别为第 j j j类。

    3.广义马氏距离

    对于正态总体,在错判损失都相等的情况下,
    q i f i ( X ) = q i ( 2 π ) m / 2 ∣ Σ ∣ 1 / 2 exp ⁡ { − 1 2 ( X − μ ( i ) ) ′ Σ − 1 ( X − μ ( i ) ) } , ln ⁡ q i f i ( X ) = C + ln ⁡ q i − 1 2 ln ⁡ ∣ Σ ∣ − 1 2 ( X − μ ( i ) ) ′ Σ − 1 ( X − μ ( i ) ) . q_if_i(X)=\frac{q_i}{(2\pi)^{m/2}|\Sigma|^{1/2}}\exp\left\{-\frac12(X-\mu^{(i)})'\Sigma^{-1}(X-\mu^{(i)}) \right\}, \\ \ln q_if_i(X)=C+\ln q_i-\frac12\ln|\Sigma|-\frac12(X-\mu^{(i)})'\Sigma^{-1}(X-\mu^{(i)}). qifi(X)=(2π)m/2Σ1/2qiexp{21(Xμ(i))Σ1(Xμ(i))},lnqifi(X)=C+lnqi21lnΣ21(Xμ(i))Σ1(Xμ(i)).
    因此,我们定义样本 X X X到总体 G i G_i Gi的广义马氏距离为
    D i 2 ( X ) = d i 2 ( X ) + ln ⁡ ∣ S ∣ − 2 ln ⁡ q i . D_i^2(X)=d^2_i(X)+\ln |S|-2\ln q_i. Di2(X)=di2(X)+lnS2lnqi.
    可以看到,当样本 X X X到总体 G i G_i Gi的广义马氏距离最小的时候,它会被归类到 G i G_i Gi。因此,在每一类都是多元正态总体,且错判损失相等的情况下,用广义马氏距离替代马氏距离,贝叶斯判别的解与直接判别法是一样的。

    回顾总结

    1. 贝叶斯判别法,是在定义了先验概率 q i q_i qi与错判损失 L ( j ∣ i ) L(j|i) L(ji)的情况下,使平均错判损失最小的判别准则。

    2. 在贝叶斯判别的条件下,决定样本 X X X应该判别到某一类的,是某一类的平均错判损失,即
      h j ( X ) = ∑ i = 1 k q i L ( j ∣ i ) f i ( X ) . h_j(X)=\sum_{i=1}^k q_iL(j|i)f_i(X). hj(X)=i=1kqiL(ji)fi(X).
      样本 X X X被判别到平均错判损失最小的一类,即 X ∈ G j ⇔ h j ( X ) ≤ h i ( X ) , ∀ i X\in G_j\Leftrightarrow h_j(X)\le h_i(X),\forall i XGjhj(X)hi(X),i

    3. 如果 L ( j ∣ i ) = 1 − δ i j L(j|i)=1-\delta_{ij} L(ji)=1δij,也就是错判的损失为1,正判的损失为0,那么判别函数可以化简为
      q i f i ( X ) . q_if_i(X). qifi(X).
      样本 X X X被判别到 q i f i ( X ) q_if_i(X) qifi(X)最大的一类,即 X ∈ G j ⇔ q j f j ( X ) ≥ q i f i ( X ) , ∀ i X\in G_j\Leftrightarrow q_jf_j(X)\ge q_if_i(X),\forall i XGjqjfj(X)qifi(X),i

    4. 定义广义马氏距离为
      D i 2 ( X ) = d i 2 ( X ) + ln ⁡ ∣ S ∣ − 2 ln ⁡ ∣ q i ∣ , D_i^2(X)=d_i^2(X)+\ln |S|-2\ln |q_i|, Di2(X)=di2(X)+lnS2lnqi,
      对于错判损失为 L ( j ∣ i ) = 1 − δ i j L(j|i)=1-\delta_{ij} L(ji)=1δij的正态总体判别,基于广义马氏距离的直接判别法就是贝叶斯准则下的最优判别法。

    展开全文
  • 浅谈贝叶斯判别(Bayes)

    千次阅读 2019-02-11 16:51:34
    在现实世界中,由许多客观现象的发生,就每一次观察和测量来说,即使在基本条件保持不变的情况下也具有不确定性。...这就是贝叶斯判别的基本出发点。 简介 给定一个输入x\mathbf{x}x,我们想要确定它属...

    在现实世界中,由许多客观现象的发生,就每一次观察和测量来说,即使在基本条件保持不变的情况下也具有不确定性。只有在大量重复的观察下,其结果才能呈现出某种规律性,即对它们观察到的特征具有统计特性。特征值不再是一个确定的向量,而是一个随机向量。此时,只能利用模式集的统计特性来分类,以使分类器发生错误的概率最小。这就是贝叶斯判别的基本出发点。

    简介

    给定一个输入 x \mathbf{x} x,我们想要确定它属于 w 1 w_1 w1类还是 w 2 w_2 w2类。而根据贝叶斯判别思想,我们需要观察或者说,计算 x x x是来自于 w 1 w_1 w1类的概率大还是来自 w 2 w_2 w2类的概率大。因此,若

    P ( w 1 ∣ x ) &gt; P ( w 2 ∣ x ) P(w_1|\mathbf{x}) &gt; P(w_2 | \mathbf{x}) P(w1x)>P(w2x)

    x ∈ w 1 \mathbf{x} \in w_1 xw1。相反,则 x ∈ w 2 \mathbf{x} \in w_2 xw2

    而由贝叶斯定理可知,后验概率 P ( w i ∣ x ) P(w_i|\mathbf{x}) P(wix)可由类别的先验概率 P ( w i ∣ x ) P(w_i|\mathbf{x}) P(wix) x \mathbf{x} x的先验概率 P ( x ∣ w i ) P(\mathbf{x}|w_i) P(xwi)得到,

    P ( w i ∣ x ) = P ( x ∣ w i ) P ( w i ) P ( x ) = P ( x ∣ w i ) P ( w i ) ∑ i = 1 2 P ( x ∣ w i ) P ( w i ) P(w_i | \mathbf{x}) = \frac{P(\mathbf{x}|w_i)P(w_i)}{P(\mathbf{x})} = \frac{P(\mathbf{x}|w_i)P(w_i)}{\sum_{i = 1}^{2}P(\mathbf{x}|w_i)P(w_i)} P(wix)=P(x)P(xwi)P(wi)=i=12P(xwi)P(wi)P(xwi)P(wi)

    这里的 P ( x ∣ w i ) P(\mathbf{x} | w_i) P(xwi)也被称为似然函数。因此,原先判别式可改写为,若

    P ( x ∣ w 1 ) P ( w 1 ) &gt; P ( x ∣ w 2 ) P ( w 2 ) P(\mathbf{x}|w_1)P(w_1) &gt; P(\mathbf{x}|w_2)P(w_2) P(xw1)P(w1)>P(xw2)P(w2)

    x ∈ w 1 \mathbf{x} \in w_1 xw1。相反,则 x ∈ w 2 \mathbf{x} \in w_2 xw2

    贝叶斯最小风险判别

    前面说到,贝叶斯判别函数的基本思想或者说优化目标是使得分类器发生的错误率最小。然而,实际任务中,不同的错误所带来的风险也是不同的。比如,将一个身患感冒的病人误诊为肺癌和误诊为肺炎所造成的风险是肯定不同的。因此,当考虑到对于某一类的错误判决要比对另一类的判决更为关键时,就需要把最小错误概率的贝叶斯判别做一些修正,提出条件平均风险。

    M M M类问题,我们令 r i j ( x ) r_{ij}(\mathbf{x}) rij(x)表示将本属于 w i w_i wi类的样本判定为 w j w_j wj类的条件风险,则有

    r j ( x ) = 1 P ( x ) ∑ i = 1 M L i j P ( x ∣ w i ) P ( w i ) r_j(\mathbf{x}) = \frac{1}{P(\mathbf{x})}\sum_{i = 1}^{M}L_{ij}P(\mathbf{x}|w_i)P(w_i) rj(x)=P(x)1i=1MLijP(xwi)P(wi)

    这里, P ( x ) P(\mathbf{x}) P(x)是公共项,可以舍去。 L i j L_{ij} Lij为将样本应属于 w i w_i wi类的模式判别成属于 w j w_j wj类的代价。因此,上式可简写为,

    r j ( x ) = ∑ i = 1 M L i j P ( x ∣ w i ) P ( w i ) r_j(\mathbf{x}) = \sum_{i = 1}^{M}L_{ij}P(\mathbf{x}|w_i)P(w_i) rj(x)=i=1MLijP(xwi)P(wi)

    这也是贝叶斯分类器,只是它的判别方法不是按错误概率最小作为标准,而是按平均条件风险作为标准。

    有意思地是,假设代价 L L L判对失分不计分,判错失分记为1分,则有

    L i j = { 0 , if i = j 1 , if i ≠ j L_{ij} = \left\{ \begin{aligned} &amp; 0, \textbf{if} \quad i = j \\ &amp; 1, \textbf{if} \quad i \ne j \\ \end{aligned} \right. Lij={0,ifi=j1,ifi̸=j

    则条件风险可写成
    r j ( x ) = ∑ i = 1 M L i j P ( x ∣ w i ) P ( w i ) = L 1 j P ( x ∣ w 1 ) P ( w 1 ) + ⋯ + L j j P ( x ∣ w j ) P ( w j ) + ⋯ + L M j P ( x ∣ w M ) P ( w M ) = ∑ i = 1 M P ( x ∣ w i ) P ( w i ) − P ( x ∣ w j ) P ( w j ) = P ( x ) − P ( x ∣ w j ) P ( w j ) \begin{aligned} r_j(\mathbf{x}) &amp;= \sum_{i = 1}^{M}L_{ij}P(\mathbf{x} | w_i)P(w_i) \\ &amp;= L_{1j}P(\mathbf{x}|w_1)P(w_1) + \dots + L_{jj}P(\mathbf{x}|w_j)P(w_j) + \dots + L_{Mj}P(\mathbf{x}|w_M)P(w_M) \\ &amp;= \sum_{i = 1}^{M}P(\mathbf{x} | w_i)P(w_i) - P(\mathbf{x} | w_j)P(w_j) \\ &amp;= P(\mathbf{x}) - P(\mathbf{x} | w_j)P(w_j)\\ \end{aligned} rj(x)=i=1MLijP(xwi)P(wi)=L1jP(xw1)P(w1)++LjjP(xwj)P(wj)++LMjP(xwM)P(wM)=i=1MP(xwi)P(wi)P(xwj)P(wj)=P(x)P(xwj)P(wj)

    假如样本 x ∈ w 1 \mathbf{x} \in w_1 xw1,由 r i ( x ) &lt; r j ( x ) r_i(x) &lt; r_j(x) ri(x)<rj(x)可得, P ( x ∣ w i ) P ( w i ) &gt; P ( x ∣ w j ) P ( w j ) P(\mathbf{x} | w_i)P(w_i) &gt; P(\mathbf{x} | w_j)P(w_j) P(xwi)P(wi)>P(xwj)P(wj)
    。不难看出,当误判风险一致时,贝叶斯最小风险判别与一般的贝叶斯判别是一致的。

    正态分布模式的贝叶斯分类器

    通过之前的内容,我们现在对贝叶斯判别已经有了一定的认识。因此,我们不妨进一步地从概率分布的角度分析贝叶斯判别。在实际问题中,常常要研究一个随机变量 ξ \xi ξ取值小于某一数值 x x x的概率,这概率是 x x x的函数,称这种函数为随机变量ξ的分布函数,简称分布函数,记作 F ( x ) F(x) F(x),即 F ( x ) = P ( ξ &lt; x ) ( − ∞ &lt; x &lt; + ∞ ) F(x) = P(\xi &lt; x)(-\infty &lt; x &lt;+ \infty) F(x)=P(ξ<x)(<x<+),由它并可以决定随机变量落入任何范围内的概率。说人话,就是刻画了概率值的分布情况。而常用的有正态分布,均匀分布等。在这里,我们讨论正太分布情况。

    M M M种模式类别的多变量正态类密度函数为例,

    P ( x ∣ w i ) = 1 ( 2 π ) n 2 ∣ C i ∣ 1 2 exp ⁡ ( − 1 2 ( x − m i ) T C i − 1 ( x − m i ) ) , i = 1 , 2 , … , M P(\mathbf{x} | w_i) = \frac{1}{(2\pi)^{\frac{n}{2}}|C_i|^{\frac{1}{2}}} \exp(-\frac{1}{2}(\mathbf{x} - \mathbf{m}_i)^{T}C_i^{-1}(\mathbf{x} - \mathbf{m}_i)), i = 1, 2, \dots, M P(xwi)=(2π)2nCi211exp(21(xmi)TCi1(xmi)),i=1,2,,M

    其中,每一类模式的分布密度都完全被其均值向量 m i \mathbf{m}_i mi和协方差矩阵 C i C_i Ci所规定。

    而我们已知判别函数可写成如下形式,

    d i ( x ) = P ( x ∣ w i ) P ( w i ) , i = 1 , … , M d_i(\mathbf{x}) = P(\mathbf{x} | w_i)P(w_i), i = 1, \dots, M di(x)=P(xwi)P(wi),i=1,,M

    对于正态密度函数,可取自然对数的形式以方便计算(因为自然对数是单调递增的,取对数后不影响相应的分类性能),则有

    d i ( x ) = ln ⁡ P ( x ∣ w i ) + ln ⁡ P ( w i ) , i = 1 , 2 , … , M d_i(\mathbf{x}) = \ln P(\mathbf{x} | w_i) + \ln P(w_i), i = 1, 2, \dots, M di(x)=lnP(xwi)+lnP(wi),i=1,2,,M

    代入正太密度函数,有

    d i ( x ) = ln ⁡ P ( w i ) − n 2 ( ln ⁡ 2 π ) − 1 2 ln ⁡ ∣ C i ∣ − 1 2 ( x − m i ) T C i − 1 ( x − m i ) , i = 1 , 2 , … , M d_i(\mathbf{x}) = \ln P(w_i) - \frac{n}{2}(\ln 2\pi) - \frac{1}{2}\ln |C_i| - \frac{1}{2}(\mathbf{x} - \mathbf{m}_i)^{T}C_i^{-1}(\mathbf{x} - \mathbf{m}_i), i = 1, 2, \dots, M di(x)=lnP(wi)2n(ln2π)21lnCi21(xmi)TCi1(xmi),i=1,2,,M

    去掉常数项,有

    d i ( x ) = ln ⁡ P ( w i ) − 1 2 ln ⁡ ∣ C i ∣ − 1 2 ( x − m i ) T C i − 1 ( x − m i ) d_i(\mathbf{x}) = \ln P(w_i) - \frac{1}{2}\ln |C_i|- \frac{1}{2}(\mathbf{x} - \mathbf{m}_i)^{T}C_i^{-1}(\mathbf{x} - \mathbf{m}_i) di(x)=lnP(wi)21lnCi21(xmi)TCi1(xmi)

    即为正态分布模式的贝叶斯判别函数。

    总结

    1. 优点

      • 朴素贝叶斯基于贝叶斯定理和属性条件独立性假设,具有坚实的理论基础。
      • 朴素贝叶斯简单易懂,模型参数少。
      • 朴素贝叶斯对于缺失数据,无关属性,噪声点等问题不敏感。
    2. 缺点

      • 现实生活中,样本分布往往不满足属性条件独立性假设,而这恰恰是朴素贝叶斯算法的基本出发点

    参考文献

    黄庆明,《第二章.ppt》

    百度百科

    展开全文
  • 判别分析是用以判别个体所属群体的一种统计方法,它产生于20世纪30年代,近年来,在许多现代自然科学的各个分支和技术部门中,得到广泛的应用. 例如,利用计算机对一个人是否有心脏病进行诊断时,可以取一批没有心脏...

    前言

    判别分析是用以判别个体所属群体的一种统计方法,它产生于20世纪30年代,近年来,在许多现代自然科学的各个分支和技术部门中,得到广泛的应用.

    例如,利用计算机对一个人是否有心脏病进行诊断时,可以取一批没有心脏病的人,测其p个指标的数据,然后再取一批已知患有心胜病的人,同样也测得p个相同指标的数据,利用这些数据建一个判别函数,并求出相应的临界值,这时对于需要进行诊断的人,也同样测其p个指标的数据,将其代入判别函数,求得判别得分,再依判别临界值,即可判断此人是属于有心脏病的那一群体,还是属于没有心脏病的那一群体。又如在考古学中,对化石及文物年代的判断;在地质学中,判断是有矿还是无矿;在量管理中,判断某种产品是合格品,还是不合格品;在植物学中,对于新发现的一种植物,判断其属于那一科.总之判别分析方法在很多学科中有着广泛的应用.

    判别方法有多种,这里主要介绍的是最常用的判别方法,而且是两类群体的判别方法.

    1 距离判别

    在这里插入图片描述
    在这里插入图片描述

    1.1 双群体

    1.1.1 理论推导

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    1.1.2 R语言实现

    discriminiant.distance.R(距离判别函数)

    discriminant.distance <- function(TrnX1,TrnX2,TstX=NULL,var.equal=FALSE){
      # 输入变量处理
      if (is.null(TstX) == TRUE) TstX <- rbind(TrnX1,TrnX2)
      if (is.vector(TstX) == TRUE) 
        TstX <- t(as.matrix(TstX))
      else if (is.matrix(TstX) != TRUE)
        TstX <- as.matrix(TstX)
      if (is.matrix(TrnX1) == FALSE)
        TrnX1 <- as.matrix(TrnX1)
      if (is.matrix(TrnX2) != TRUE)
        TrnX2 <- as.matrix(TrnX2)
      
      # 
      nx <- nrow(TstX)    # 需要用以预测的集合的大小,或者说测试集的大小
      # 生成长度为nx的0向量,用以存储预测的标签
      blong <- matrix(rep(0,nx),nrow=1,byrow=TRUE,dimnames=list("blong",1:nx))
      # 两个群体的均值向量
      mu1 <- colMeans(TrnX1); mu2 <- colMeans(TrnX2) 
      # 两群体同方差
      if (var.equal == TRUE || var.equal == T){
        # 计算混合样本方差
        S <- var(rbind(TrnX1,TrnX2))
        # 到第二群体的马氏距离减去到第一群体的马氏距离——>判别函数W(x)
        W <- mahalanobis(TstX, mu2, S) - mahalanobis(TstX, mu1, S)
      }
      # 两群体异方差
      else{
        S1 <- var(TrnX1); S2 <- var(TrnX2)
        W <- mahalanobis(TstX, mu2, S2) - mahalanobis(TstX, mu1, S1)
      }  
      for (i in 1:nx){
        if (W[i] > 0)
          blong[i] <- 1
        else
          blong[i] <- 2
      }
      blong
    }
    

    在程序中,输入变量Trnx1Trnx2表示X1类、X2类训练样本,其输入格式是数据框,或矩阵(样本按行输入),输入变量TstX是待测样本,其输入格式是数据框,或矩阵(样本按行输入),或向量(一个待测样本).如果不输入Tstx(缺省值),则待测样本为两个训练样本之和,即计算训练样本的回代情况。输入变量var.equal是逻辑变量,var.equal=TRUE表示两个总体的协方差阵相同;否则(缺省值)为不同。函数的输出是由“1”和“2”构成的的一维矩阵,“1”表示待测样本属于X1类,“2”表示待测样本属于X2类。

    在上述程序中,用到Mahalanobis距离函数mahalanobis(),该函数的使用格式为mahalanobis(x, center, cov, inverted=FLASE)
    其中x是由样本数据构成的向量或矩阵(p维),center为样本中心,cov为样本的协方差阵.

    1.1.3 实例分析

    在这里插入图片描述
    在这里插入图片描述
    R语言代码

    classX1 <- data.frame(
      x1=c(6.60, 6.60, 6.10, 6.10, 8.40, 7.2, 8.40, 7.50, 7.50, 8.30, 7.80, 7.80),
      x2=c(39.00,39.00, 47.00, 47.00, 32.00, 6.0, 113.00, 52.00, 52.00,113.00,172.00,172.00),
      x3=c(1.00, 1.00, 1.00, 1.00, 2.00, 1.0, 3.50, 1.00,3.50, 0.00, 1.00, 1.50),
      x4=c(6.00, 6.00, 6.00, 6.00, 7.50, 7.0, 6.00, 6.00,7.50, 7.50, 3.50, 3.00),
      x5=c(6.00, 12.00, 6.00, 12.00, 19.00, 28.0, 18.00, 12.00,6.00, 35.00, 14.00, 15.00),
      x6=c(0.12, 0.12, 0.08, 0.08, 0.35, 0.3, 0.15, 0.16,0.16, 0.12, 0.21, 0.21),
      x7=c(20.00,20.00, 12.00, 12.00, 75.00, 30.0, 75.00, 40.00,40.00,180.00, 45.00, 45.00)
    )
    classX2<-data.frame(
      x1=c(8.40, 8.40, 8.40, 6.3, 7.00, 7.00, 7.00, 8.30,
           8.30, 7.2, 7.2, 7.2, 5.50, 8.40, 8.40, 7.50,
           7.50, 8.30, 8.30, 8.30, 8.30, 7.80, 7.80),
      x2=c(32.0 ,32.00, 32.00, 11.0, 8.00, 8.00, 8.00,161.00,
           161.0, 6.0, 6.0, 6.0, 6.00,113.00,113.00, 52.00,
           52.00, 97.00, 97.00,89.00,56.00,172.00,283.00),
      x3=c(1.00, 2.00, 2.50, 4.5, 4.50, 6.00, 1.50, 1.50,
           0.50, 3.5, 1.0, 1.0, 2.50, 3.50, 3.50, 1.00,
           1.00, 0.00, 2.50, 0.00, 1.50, 1.00, 1.00),
      x4=c(5.00, 9.00, 4.00, 7.5, 4.50, 7.50, 6.00, 4.00,
           2.50, 4.0, 3.0, 6.0, 3.00, 4.50, 4.50, 6.00,
           7.50, 6.00, 6.00, 6.00, 6.00, 3.50, 4.50),
      x5=c(4.00, 10.00, 10.00, 3.0, 9.00, 4.00, 1.00, 4.00,
           1.00, 12.0, 3.0, 5.0, 7.00, 6.00, 8.00, 6.00,
           8.00, 5.00, 5.00,10.00,13.00, 6.00, 6.00),
      x6=c(0.35, 0.35, 0.35, 0.2, 0.25, 0.25, 0.25, 0.08,
           0.08, 0.30, 0.3, 0.3, 0.18, 0.15, 0.15, 0.16,
           0.16, 0.15, 0.15, 0.16, 0.25, 0.21, 0.18),
      x7=c(75.00,75.00, 75.00, 15.0,30.00, 30.00, 30.00, 70.00,
           70.00, 30.0, 30.0, 30.0,18.00, 75.00, 75.00, 40.00,
           40.00,180.00,180.00,180.00,180.00,45.00,45.00)
    )
    
    path <- 'E:\\桌面文档\\学习文件\\大三\\多元统计\\统计建模R语言书籍'
    setwd(path)   # 设置工作路径(work directory)
    source('distanceDiscriminant.R')
    # 方差相同
    discriminant.distance(classX1,classX2,var.equal=TRUE)
    # 方差不同
    discriminant.distance(classX1,classX2)
    

    在认为两总体协方差阵相同的情况下,将训练样本回代进行判别,有三个点判错,分别是第9号样本、第28号样本和第29号样本.

    在认为两总体协方差阵不同的情况下,将训练样本回代进行判别,只有一个点判错,是第9号样本。

    1.2 多群体

    1.2.1 理论推导

    在这里插入图片描述

    1.2.2 R语言实现

    distinguish.distance(多分类问题的距离判别函数)

    distinguish.distance <- function(TrnX,TrnG,TstX=NULL,var.equal=FALSE){
      if (is.factor(TrnG) == FALSE){
        mx <- nrow(TrnX); mg <- nrow(TrnG)
        TrnX <- rbind(TrnX, TrnG)
        TrnG <- factor(rep(1:2,c(mx,mg)))  # 1重复mx遍,2重复mg遍
      }
      if (is.null(TstX) == TRUE) TstX <- TrnX
      if (is.vector(TstX) == TRUE) 
        TstX <-t(as.matrix(TstX))
      else if (is.matrix(TstX) != TRUE)
        TstX <- as.matrix(TstX)
      if (is.matrix(TrnX) != TRUE)
        TrnX <- as.matrix(TrnX)
      
      nx <- nrow(TstX)
      # blong用于存放预测值
      blong <- matrix(rep(0,nx),nrow=1,dimnames=list("blong",1:nx))
      g <- length(levels((TrnG)))    # 计算群体类别个数
      mu <- matrix(0,nrow=g,ncol=ncol(TrnX))
      # 每一个群体都有一个均值
      for(i in 1:g)   
        mu[i,] <- colMeans(TrnX[TrnG == i,])
      print(mu)
      # 计算马氏距离
      D <- matrix(0,nrow=g,ncol=nx)
      if (var.equal == TRUE || var.equal == T){
        for (i in 1:g)    # 样本到每一个类别的马氏距离
          D[i,] <- mahalanobis(TstX,mu[i,],var(TrnX))  # 混合样本方差
      }
      
      else{
        for (i in 1:g)
          D[i,] <- mahalanobis(TstX, mu[i,],var(TrnX[TrnG == i,]))
      }
      print(D)
      for (j in 1:nx){   # 分别判别每一个样本属于哪一个类别
        dmin <- Inf
        for (i in 1:g){    # 遍历每一个类别,找出最小距离
          if (D[i,j] < dmin){
            dmin <- D[i,j];
            blong[j] <- i
          }
        }
      }
      blong
    }
    

    程序分别考虑了总体协方差阵相同和总体协方差阵不同的两种情况。输入变量TrnX表示训练样本,其输入格式是矩阵(样本按行输入),或数据框.TrnG是因子变量,表示输入训练样本的分类情况.输入变量TstX是待测样本,其输入格式是矩阵(样本按行输入),或数据框,或向量(一个待测样本).如果不输入TstX(缺省值),则待测样本为训练样本.输入变量var.equal是逻辑变量,var.equal=TRUE表示计算时认为总体协方差阵是相同的;否则(缺省值)是不同的.函数的输出是由数字构成的的一维矩阵,数字表示相应的类.为了与前一个程序兼容,对于二分类问题,也可以按照discriminiant.distance函数的输入格式输入.

    1.2.3 实例分析

    在这里插入图片描述
    R软件中提供了Iris数据,数据的前四列是数据的四个属性,第五列标明数据属于哪一类.

    代码:

    ### 多群体距离判别
    X <- iris[,1:4]
    G <- gl(3,50)    # 与rep函数有所不同,gl是将数字先重复再接起来
    path <- 'E:\\桌面文档\\学习文件\\大三\\多元统计\\统计建模R语言书籍'
    setwd(path)
    source('distanceDiscriminant2.R')
    distinguish.distance(X,G)
    

    结果:

          1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
    blong 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
          25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
    blong  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
          46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
    blong  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
          67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
    blong  2  2  2  2  3  2  3  2  2  2  2  2  2  2  2  2  2  3  2  2  2
          88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
    blong  2  2  2  2  2  2  2  2  2  2  2  2   2   3   3   3   3   3
          106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
    blong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3
          121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
    blong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3
          136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
    blong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3
    

    从计算结果可以看出,只有第71号样本、第73号样本和第84号样本错判,回代的判别正确率为147/150=98%.

    2 贝叶斯判别

    Bayes判别是假定对研究对象已有一定的认识,这种认识常用先验概率来描述,当取得样本后,就可以用样本来修正已有的先验概率分布,得出后验概率分布,现通过后验概率分布进行各种统计推断。

    2.1 双群体

    2.1.1 理论推导

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    2.1.2 R语言实现

    discriminant.bayes.R(双群体最小ECM贝叶斯判别函数)

    discriminant.bayes <- function(TrnX1,TrnX2,rate=1,TstX=NULL,var.equal=FALSE){
      if (is.null(TstX) == TRUE) TstX <- rbind(TrnX1,TrnX2)
      if (is.vector(TstX) == TRUE) TstX <- t(as.matrix(TstX))
      else if (is.matrix(TstX) != TRUE)
        TstX <- as.matrix(TstX)
      if (is.matrix(TrnX1) != TRUE) TrnX1 <- as.matrix(TrnX1)
      if (is.matrix(TrnX2) != TRUE) TrnX2 <- as.matrix(TrnX2)
      
      nx <- nrow(TstX)
      blong <- matrix(rep(0,nx), nrow=1, byrow=TRUE, dimnames = list("blong",1:nx))
      mu1 <- colMeans(TrnX1); mu2 <- colMeans(TrnX2)
      if (var.equal == TRUE || var.equal == T){
        S <- var(rbind(TrnX1,TrnX2))
        beta <- 2*log(rate)   # W前面的1/2乘到beta上面去了
        W <- mahalanobis(TstX,mu2,S) - mahalanobis(TstX,mu1,S)
      }
      else{
        S1 <- var(TrnX1); S2 <- var(TrnX2)
        beta <- 2*log(rate) + log(det(S1)/det(S2))
        W <- mahalanobis(TstX, mu2, S2) - mahalanobis(TstX, mu1, S1)
      }
      
      for (i in 1:nx){
        if (W[i] > beta)
          blong[i] <- 1
        else
          blong[i] <- 2
      }
      blong
    }
    

    在程序中,输入变量TrnX1Trnx2表示X1类、X2类训练样本,其输入格式是数据框,或矩阵(样本按行输入). r a t e = L ( 1 ∣ 2 ) L ( 2 ∣ 1 ) ⋅ p 2 p 1 rate=\frac{L(1|2)}{L(2|1)}\cdot \frac{p_2}{p_1} rate=L(21)L(12)p1p2,缺省值为1.TstX是待测样本,其输入格式是数据框,或矩阵(样本按行输入),或向量(一个待测样本).如果不输入TstX(缺省值),则待测样本为两个训练样本之和,即计算训练样本的回代情况.输入变量var.equal是逻辑变量,var.equal=TRUE表示认为两总体的协方差阵是相同的;否则(缺省值)是不同的.函数的输出是由“1”和“2”构成的的一维矩阵,“1”表示待测样本属于X1类,“2”表示待测样本属于X2类。`

    2.1.3 实例分析

    在这里插入图片描述
    代码:

    TrnX1<-matrix(
      c(24.8, 24.1, 26.6, 23.5, 25.5, 27.4,
        -2.0, -2.4, -3.0, -1.9, -2.1, -3.1),
      ncol=2)
    TrnX2<-matrix(
      c(22.1, 21.6, 22.0, 22.8, 22.7, 21.5, 22.1, 21.4,
        -0.7, -1.4, -0.8, -1.6, -1.5, -1.0, -1.2, -1.3),
      ncol=2)
    path <- 'E:\\桌面文档\\学习文件\\大三\\多元统计\\统计建模R语言书籍'
    setwd(path)
    source('bayesDiscriminant.R')
    discriminant.bayes(TrnX1,TrnX2,rate=8/6,var.equal=TRUE)
    

    结果:

          1 2 3 4 5 6 7 8 9 10 11 12 13 14
    blong 1 1 1 2 1 1 2 2 2  2  2  2  2  2
    

    第4号样本被错判。

    2.2 多群体

    2.2.1 理论推导

    在这里插入图片描述

    2.2.2 R语言实现

    ditinguish.bayes.R(多群体ECM最小准则贝叶斯分类函数)

    distinguish.bayes <- function(TrnX,TrnG,p=rep(1,length(levels(TrnG))),TstX=NULL,var.equal=FALSE){
      if (is.factor(TrnG) == FALSE){
        mx <- nrow(TrnX); mg <- nrow(TrnG)
        TrnX <- rbind(TrnX, TrnG)
        TrnG <- factor(rep(1:2, c(mx, mg)))
      }
      if (is.null(TstX) == TRUE) TstX <- TrnX
      if(is.vector(TstX) == TRUE) TstX <- t(as.matrix(TstX))
      else if (is.matrix(TstX) != TRUE)
        TstX <- as.matrix(TstX)
      if (is.matrix(TrnX) != TRUE) TrnX <- as.matrix(TrnX)
      
      nx <- nrow(TstX)
      blong <- matrix(0,nrow=1,ncol=nx,dimnames = list("blong",1:nx))
      g <- length(levels(TrnG))
      mu <- matrix(0,nrow=g,ncol=ncol(TrnX))
      for (i in 1:g)
        mu[i,] <- colMeans(TrnX[TrnG==i,])
      D <- matrix(0,nrow=g,ncol=nx)  
      if (var.equal == TRUE || var.equal == T){
        for (i in 1:g){
          d2 <- mahalanobis(TstX, mu[i,],var(TrnX))
          D[i,] <- d2 - 2*log(p[i])
        }
      }
      else{
        for (i in 1:g){
          S <- var(TrnX[TrnG == i,])
          d2 <- mahalanobis(TstX, mu[i,],S)
          D[i,] = d2 - 2*log(p[i]) - log(det(S))
        }
      }
      for (j in 1:nx){
        dmin <- Inf
        for (i in 1:g){
          if (D[i,j] < dmin){
            dmin <- D[i,j]; blong[j] <- i
          }
        }
      }
      blong
    }
    

    程序分别考虑了总体协方差阵相同和协方差阵不同的情况.输入变量Trnx表示训练样本,其输入格式是矩阵(样本按行输入),或数据框.TrnG是因子变量,表示训练样本的分类情况.输入变量p是先验概率,缺省值均为1.输入变量TstX是待测样本,其输入格式是矩阵(样本按行输入),或数据框,或向量(一个待测样本).如果不输入TstX(缺省值),则待测样本为训练样本.输入变量var.equal是逻辑变量,var.equal=TRUE表示认为总体协方差阵是相同的;否则(缺省值)是不同的.函数的输出是由数字构成的的一维矩阵,数字表示相应的类.为了与前面两总体的判别程序兼容,对于二分类问题,也可以按照discriminiant.bayes函数的输入格式输入.

    2.2.3 实例分析

    在这里插入图片描述
    代码:

    ## 多类别
    X <- iris[,1:4]
    G <- gl(3,50)
    path <- 'E:\\桌面文档\\学习文件\\大三\\多元统计\\统计建模R语言书籍'
    setwd(path)
    source('bayesDiscriminant2.R')
    blong <- distinguish.bayes(X,G)
    

    结果:

          1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
    blong 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
          25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
    blong  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
          46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
    blong  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2
          67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87
    blong  2  2  3  2  3  2  3  2  2  2  2  3  2  2  2  2  2  3  2  2  2
          88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105
    blong  2  2  2  2  2  2  2  2  2  2  2  2   2   3   3   3   3   3
          106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
    blong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3
          121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
    blong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3
          136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
    blong   3   3   3   3   3   3   3   3   3   3   3   3   3   3   3
    
    > table(blong,G)
         G
    blong  1  2  3
        1 50  0  0
        2  0 45  0
        3  0  5 50
    

    从计算结果可以看出,只有第69、71、73、78、84号样本错判,回代的判别正确率为145/150=96.67%.

    3 Fisher判别

    Fisher(费歇)判别是按类内方差尽量小,类间方差尽量大的准则来求判别函数的.在这里仅讨论两个总体的判别方法.

    3.1 理论推导

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    3.2 R语言实现

    3.2.1 自己实现

    方法1

    discriminant.fisher.R(双群体Fisher判别函数)

    discriminant.fisher <- function(TrnX1, TrnX2, TstX=NULL){
      if (is.null(TstX) == TRUE) TstX <- rbind(TrnX1,TrnX2)
      if (is.vector(TstX) == TRUE) TstX <- t(as.matrix(TstX))
      else if (is.matrix(TstX) != TRUE)
        TstX <- as.matrix(TstX)
      if (is.matrix(TrnX1) != TRUE) TrnX1 <- as.matrix(TrnX1)
      if (is.matrix(TrnX2) != TRUE) TrnX2 <- as.matrix(TrnX2)
      
      nx <- nrow(TstX)
      blong <- matrix(rep(0,nx),nrow=1,byrow=TRUE,dimnames=list("blong",1:nx))
      
      n1 <- nrow(TrnX1); n2 <- nrow(TrnX2)
      mu1 <- colMeans(TrnX1); mu2 <- colMeans(TrnX2)
      S <- (n1-1)*var(TrnX1) + (n2-1)*var(TrnX2)
      mu <- n1/(n1+n2)*mu1 + n2/(n1+n2)*mu2
      W <- (TstX - rep(1,nx) %o% mu) %*% solve(S,mu2-mu1)
      for (i in 1:nx){
        if (W[i] <= 0)
          blong[i] <- 1
        else
          blong[i] <- 2
      }
      blong
    }
    

    在程序中,输入变量TrnX1Trnx2表示X1类、X2类训练样本,其输入格式是数据框,或矩阵(样本按行输入).TstX是待测样本,其输入格式是数据框,或矩阵(样本按行输入),或向量(一个待测样本).如果不输入TstX(缺省值),则待测样本为两个训练样本之和,即计算训练样本的回代情况.函数的输出是由“1”和“2”构成的的一维矩阵,“1”表示待测样本属于X1类,“2”表示待测样本属于X2类。

    有一个需要注意的地方,rep(1,nx) %o% mu不可改为mu,两个算出来的结果是不一样的:

    (1)情况1> rep(1,nrow(classX1)) %o% mu1
                x1       x2       x3 x4    x5        x6   x7
     [1,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
     [2,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
     [3,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
     [4,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
     [5,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
     [6,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
     [7,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
     [8,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
     [9,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
    [10,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
    [11,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
    [12,] 7.358333 73.66667 1.458333  6 15.25 0.1716667 49.5
    
    此时
    > classX1 - rep(1,nrow(classX1)) %o% mu1
               x1        x2          x3   x4    x5          x6    x7
    1  -0.7583333 -34.66667 -0.45833333  0.0 -9.25 -0.05166667 -29.5
    2  -0.7583333 -34.66667 -0.45833333  0.0 -3.25 -0.05166667 -29.5
    3  -1.2583333 -26.66667 -0.45833333  0.0 -9.25 -0.09166667 -37.5
    4  -1.2583333 -26.66667 -0.45833333  0.0 -3.25 -0.09166667 -37.5
    5   1.0416667 -41.66667  0.54166667  1.5  3.75  0.17833333  25.5
    6  -0.1583333 -67.66667 -0.45833333  1.0 12.75  0.12833333 -19.5
    7   1.0416667  39.33333  2.04166667  0.0  2.75 -0.02166667  25.5
    8   0.1416667 -21.66667 -0.45833333  0.0 -3.25 -0.01166667  -9.5
    9   0.1416667 -21.66667  2.04166667  1.5 -9.25 -0.01166667  -9.5
    10  0.9416667  39.33333 -1.45833333  1.5 19.75 -0.05166667 130.5
    11  0.4416667  98.33333 -0.45833333 -2.5 -1.25  0.03833333  -4.5
    12  0.4416667  98.33333  0.04166667 -3.0 -0.25  0.03833333  -4.5
    
    (2)情况2> classX1 -  mu1
                x1        x2          x3         x4         x5           x6         x7
    1   -0.7583333  38.82833  -5.0000000 -67.666667 -43.500000 -15.13000000  18.541667
    2  -67.0666667 -10.50000 -14.2500000   4.541667   4.641667  -0.05166667  14.000000
    3    4.6416667  39.64167   0.8283333   0.000000 -67.666667 -49.42000000  -3.250000
    4    0.1000000 -26.66667 -48.5000000  -9.250000  10.541667  -7.27833333  11.828333
    5   -6.8500000  30.54167  -5.3583333   7.328333  13.000000 -73.31666667  25.500000
    6    7.0283333   0.00000 -72.6666667 -42.500000  12.750000  -1.15833333  22.641667
    7  -41.1000000  97.75000   2.0416667  -1.358333  17.828333  -5.85000000   1.333333
    8    0.1416667  51.82833  -5.0000000 -67.666667 -37.500000 -15.09000000  38.541667
    9  -66.1666667   2.50000 -11.7500000   6.041667  -1.358333  -0.01166667  34.000000
    10   6.8416667 105.64167  -0.1716667   1.500000 -38.666667 -49.38000000 164.750000
    11   1.8000000  98.33333 -48.5000000 -11.750000  12.541667  -7.14833333  44.828333
    12  -7.4500000 170.54167  -5.8583333   2.828333   9.000000 -73.45666667  -4.500000
    

    换成mu之后结果明显有问题。

    方法2

    discriminant.fisher <- function(TrnX1, TrnX2, TstX=NULL){
      if (is.null(TstX) == TRUE) TstX <- rbind(TrnX1,TrnX2)
      if (is.vector(TstX) == TRUE) TstX <- t(as.matrix(TstX))
      else if (is.matrix(TstX) != TRUE)
        TstX <- as.matrix(TstX)
      if (is.matrix(TrnX1) != TRUE) TrnX1 <- as.matrix(TrnX1)
      if (is.matrix(TrnX2) != TRUE) TrnX2 <- as.matrix(TrnX2)
      
      nx <- nrow(TstX)
      blong <- matrix(rep(0,nx),nrow=1,byrow=TRUE,dimnames=list("blong",1:nx))
    
      # 方法2
      n1 <- nrow(TrnX1); n2 <- nrow(TrnX2)
      mu1 <- colMeans(TrnX1); mu2 <- colMeans(TrnX2)   # 计算均值向量
      S <- ( (n1-1)*var(TrnX1) + (n2-1)*var(TrnX2) )/(n1+n2-2)    # 混合样本方差
      mu <- (mu1 + mu2)/2   # 这是一个行向量
      W <- (mu2-mu1) %*% solve(S) %*% t(TstX - rep(1,nx) %o% mu)
      # 类别判别
      for (i in 1:nx){
        if (W[i] <= 0)
          blong[i] <- 1
        else
          blong[i] <- 2
      }
      blong
    }
    

    3.2.2 借助MASS

    MASS包是R语言自带的一个包,不用安装。

    classX1$target <- 1
    classX2$target <- 2
    data <- rbind(classX1,classX2)   # 按行合并数据
    library(MASS)
    attach(data)
    ld <- lda(target ~ x2+x2+x3+x4+x5+x6+x7,prior=c(0.5,0.5))
    ld
    lp <- predict(ld)
    table(lp$class,data$target)
    lp$class
    detach(data)
    

    3.3 实例分析

    使用Fisher判别求解1.1.3中的实例
    代码:

    # Fisher判别
    classX1 <- data.frame(
      x1=c(6.60, 6.60, 6.10, 6.10, 8.40, 7.2, 8.40, 7.50, 7.50, 8.30, 7.80, 7.80),
      x2=c(39.00,39.00, 47.00, 47.00, 32.00, 6.0, 113.00, 52.00, 52.00,113.00,172.00,172.00),
      x3=c(1.00, 1.00, 1.00, 1.00, 2.00, 1.0, 3.50, 1.00,3.50, 0.00, 1.00, 1.50),
      x4=c(6.00, 6.00, 6.00, 6.00, 7.50, 7.0, 6.00, 6.00,7.50, 7.50, 3.50, 3.00),
      x5=c(6.00, 12.00, 6.00, 12.00, 19.00, 28.0, 18.00, 12.00,6.00, 35.00, 14.00, 15.00),
      x6=c(0.12, 0.12, 0.08, 0.08, 0.35, 0.3, 0.15, 0.16,0.16, 0.12, 0.21, 0.21),
      x7=c(20.00,20.00, 12.00, 12.00, 75.00, 30.0, 75.00, 40.00,40.00,180.00, 45.00, 45.00)
    )
    classX2<-data.frame(
      x1=c(8.40, 8.40, 8.40, 6.3, 7.00, 7.00, 7.00, 8.30,
           8.30, 7.2, 7.2, 7.2, 5.50, 8.40, 8.40, 7.50,
           7.50, 8.30, 8.30, 8.30, 8.30, 7.80, 7.80),
      x2=c(32.0 ,32.00, 32.00, 11.0, 8.00, 8.00, 8.00,161.00,
           161.0, 6.0, 6.0, 6.0, 6.00,113.00,113.00, 52.00,
           52.00, 97.00, 97.00,89.00,56.00,172.00,283.00),
      x3=c(1.00, 2.00, 2.50, 4.5, 4.50, 6.00, 1.50, 1.50,
           0.50, 3.5, 1.0, 1.0, 2.50, 3.50, 3.50, 1.00,
           1.00, 0.00, 2.50, 0.00, 1.50, 1.00, 1.00),
      x4=c(5.00, 9.00, 4.00, 7.5, 4.50, 7.50, 6.00, 4.00,
           2.50, 4.0, 3.0, 6.0, 3.00, 4.50, 4.50, 6.00,
           7.50, 6.00, 6.00, 6.00, 6.00, 3.50, 4.50),
      x5=c(4.00, 10.00, 10.00, 3.0, 9.00, 4.00, 1.00, 4.00,
           1.00, 12.0, 3.0, 5.0, 7.00, 6.00, 8.00, 6.00,
           8.00, 5.00, 5.00,10.00,13.00, 6.00, 6.00),
      x6=c(0.35, 0.35, 0.35, 0.2, 0.25, 0.25, 0.25, 0.08,
           0.08, 0.30, 0.3, 0.3, 0.18, 0.15, 0.15, 0.16,
           0.16, 0.15, 0.15, 0.16, 0.25, 0.21, 0.18),
      x7=c(75.00,75.00, 75.00, 15.0,30.00, 30.00, 30.00, 70.00,
           70.00, 30.0, 30.0, 30.0,18.00, 75.00, 75.00, 40.00,
           40.00,180.00,180.00,180.00,180.00,45.00,45.00)
    )
    
    path <- 'E:\\桌面文档\\学习文件\\大三\\多元统计\\统计建模R语言书籍'
    setwd(path)
    source('fisherDiscriminant.R')
    discriminant.fisher(classX1,classX2)
    

    结果:

          1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
    blong 1 1 1 1 1 1 1 1 1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2
          25 26 27 28 29 30 31 32 33 34 35
    blong  2  2  2  1  1  2  2  2  2  2  2
    

    将训练样本回代进行判别,有两个点判错,分别是第28,29号样本。对于多类的Fisher判别,其基本原理是相同的,最终也是可以转化为距离判别,这里就不介绍了.

    F附录——MASS包的使用

    F1 双群体

    F1.1 线性判别

    使用的数据为

    G	x1	x2	G1
    1	-1.9	3.21	-6.9	0.41	5.2	21	5	2.51	7.3	01	6.8	12.71	0.9	-5.41	-12.5	-2.51	1.5	1.31	3.8	6.82	0.2	6.22	-0.1	7.52	0.4	14.62	2.7	8.32	2.1	0.82	-4.6	4.32	-1.7	10.92	-2.6	13.12	2.6	12.82	-2.8	10

    数据可视化:

    d6.1 <- read.table("6_1.txt", header = TRUE) # 第一行作为列名
    attach(d6.1)
    plot(x1,x2)  # 绘制散点图
    text(x1,x2,G,adj=-0.5)  # 添加文本标记
    detach(d6.1)
    

    在这里插入图片描述
    判别分析代码:

    library(MASS)
    attach(d6.1)
    ld <- lda(G~x1+x2)   # 做线性判别分析
    ld
    
    Z <- predict(ld)    # 预测
    Z
    newG <- Z$class
    cbind(G,Z$x,newG)
    tab <- table(G, newG)  # 生成原始类别与预测类别的列联表
    sum(diag(prop.table(tab)))   # prop.table将列联表对应元素全部转为相应的概率
    predict(ld, data.frame(x1=8.1,x2=2.0)) # 预测一组新数据对应的输出
    detach(d6.1)
    

    结果:

    > ld
    Call:
    lda(G ~ x1 + x2)
    
    Prior probabilities of groups:
      1   2 
    0.5 0.5 
    
    Group means:
         x1   x2
    1  0.92 2.10
    2 -0.38 8.85
    
    Coefficients of linear discriminants:
              LD1
    x1 -0.1035305
    x2  0.2247957
    
    > Z
    $class
     [1] 1 1 1 1 1 2 1 1 1 1 2 2 2 2 1 2 2 2 2 2
    Levels: 1 2
    
    $posterior
                1          2
    1  0.61625867 0.38374133
    2  0.65888888 0.34111112
    3  0.89412853 0.10587147
    4  0.87143887 0.12856113
    5  0.96214822 0.03785178
    6  0.17275662 0.82724338
    7  0.98442237 0.01557763
    8  0.68514398 0.31485602
    9  0.85330562 0.14669438
    10 0.52789262 0.47210738
    11 0.43015877 0.56984123
    12 0.30676827 0.69323173
    13 0.03336323 0.96663677
    14 0.34672296 0.65327704
    15 0.88585263 0.11414737
    16 0.40213732 0.59786268
    17 0.08694507 0.91305493
    18 0.03480991 0.96519009
    19 0.08934413 0.91065587
    20 0.09926372 0.90073628
    
    $x
               LD1
    1  -0.28674901
    2  -0.39852439
    3  -1.29157053
    4  -1.15846657
    5  -1.95857603
    6   0.94809469
    7  -2.50987753
    8  -0.47066104
    9  -1.06586461
    10 -0.06760842
    11  0.17022402
    12  0.49351760
    13  2.03780185
    14  0.38346871
    15 -1.24038077
    16  0.24005867
    17  1.42347182
    18  2.01119984
    19  1.40540244
    20  1.33503926
    
    > cbind(G,Z$x,newG)
       G         LD1 newG
    1  1 -0.28674901    1
    2  1 -0.39852439    1
    3  1 -1.29157053    1
    4  1 -1.15846657    1
    5  1 -1.95857603    1
    6  1  0.94809469    2
    7  1 -2.50987753    1
    8  1 -0.47066104    1
    9  1 -1.06586461    1
    10 1 -0.06760842    1
    11 2  0.17022402    2
    12 2  0.49351760    2
    13 2  2.03780185    2
    14 2  0.38346871    2
    15 2 -1.24038077    1
    16 2  0.24005867    2
    17 2  1.42347182    2
    18 2  2.01119984    2
    19 2  1.40540244    2
    20 2  1.33503926    2
    
    > sum(diag(prop.table(tab)))   # prop.table将列联表对应元素全部转为相应的概率
    [1] 0.9
    
    > predict(ld, data.frame(x1=8.1,x2=2.0)) # 预测一组新数据对应的输出
    $class
    [1] 1
    Levels: 1 2
    
    $posterior
              1          2
    1 0.9327428 0.06725717
    
    $x
            LD1
    1 -1.591809
    

    F1.2 二次判别

    此时各群体协方差阵不同,使用的数据为:

    G	Q	C	P
    1	8.3	4	29
    1	9.5	7	68
    1	8	5	39
    1	7.4	7	50
    1	8.8	6.5	55
    1	9	7.5	58
    1	7	6	75
    1	9.2	8	82
    1	8	7	67
    1	7.6	9	90
    1	7.2	8.5	86
    1	6.4	7	53
    1	7.3	5	48
    2	6	2	20
    2	6.4	4	39
    2	6.8	5	48
    2	5.2	3	29
    2	5.8	3.5	32
    2	5.5	4	34
    2	6	4.5	36
    

    数据可视化为:

    d6.2 <- read.table("6_2.txt",header = TRUE)
    attach(d6.2)
    plot(Q,C)
    text(Q,C,G,adj=-0.5)
    
    plot(Q,P)
    text(Q,P,G,adj=-0.5)
    
    plot(C,P)
    text(C,P,G,adj=-0.5)
    

    在这里插入图片描述

    判别分析代码:

    library(MASS)
    qd <- qda(G~Q+C+P)
    qd
    
    Z <- predict(qd) # 预测
    newG <- Z$class
    cbind(G,newG)  # 按列合并
    
    predict(qd, data.frame(Q=8,C=7.5,P=65))   # 预测一组新数据
    
    > qd
    Call:
    qda(G ~ Q + C + P)
    
    Prior probabilities of groups:
       1    2 
    0.65 0.35 
    
    Group means:
             Q        C        P
    1 7.976923 6.730769 61.53846
    2 5.957143 3.714286 34.00000
    
    > cbind(G,newG)  # 按列合并
          G newG
     [1,] 1    1
     [2,] 1    1
     [3,] 1    1
     [4,] 1    1
     [5,] 1    1
     [6,] 1    1
     [7,] 1    1
     [8,] 1    1
     [9,] 1    1
    [10,] 1    1
    [11,] 1    1
    [12,] 1    1
    [13,] 1    1
    [14,] 2    2
    [15,] 2    2
    [16,] 2    2
    [17,] 2    2
    [18,] 2    2
    [19,] 2    2
    [20,] 2    2
    
    > predict(qd, data.frame(Q=8,C=7.5,P=65))   # 预测一组新数据
    $class
    [1] 1
    Levels: 1 2
    
    $posterior
              1            2
    1 0.9998462 0.0001537705
    

    F2 多总体

    F2.1 线性判别

    此时协方差阵相等,使用的数据为

    G	Q	C	P
    1	8.3	4	29
    1	9.5	7	68
    1	8	5	39
    1	7.4	7	50
    1	8.8	6.5	55
    2	9	7.5	58
    2	7	6	75
    2	9.2	8	82
    2	8	7	67
    2	7.6	9	90
    2	7.2	8.5	86
    2	6.4	7	53
    2	7.3	5	48
    3	6	2	20
    3	6.4	4	39
    3	6.8	5	48
    3	5.2	3	29
    3	5.8	3.5	32
    3	5.5	4	34
    3	6	4.5	36
    

    数据可视化为

    d6.3 <- read.table("6_3.txt",header = TRUE)
    attach(d6.3)
    plot(Q,C)
    text(Q,C,G,adj=-0.5)
    
    plot(Q,P)
    text(Q,P,G,adj=-0.5)
    
    plot(C,P)
    text(C,P,G,adj=-0.5)
    

    在这里插入图片描述

    线性判别分析为

    ld <- lda(G~Q+C+P)
    ld
    
    Z <- predict(ld)
    newG <- Z$class
    cbind(G,Z$x,newG)
    
    tab <- table(G,newG)
    diag(prop.table(tab,1))
    
    sum(diag(prop.table(tab)))
    
    plot(Z$x, main="陈炜")
    text(Z$x[,1], Z$x[,2],G, adj=-0.5,cex=0.75)
    
    predict(ld,data.frame(Q=8,C=7.5,P=65))
    

    结果:

    > ld
    Call:
    lda(G ~ Q + C + P)
    
    Prior probabilities of groups:
       1    2    3 
    0.25 0.40 0.35 
    
    Group means:
             Q        C      P
    1 8.400000 5.900000 48.200
    2 7.712500 7.250000 69.875
    3 5.957143 3.714286 34.000
    
    Coefficients of linear discriminants:
              LD1         LD2
    Q -0.81173396  0.88406311
    C -0.63090549  0.20134565
    P  0.01579385 -0.08775636
    
    Proportion of trace:
       LD1    LD2 
    0.7403 0.2597 
    
    > cbind(G,Z$x,newG)
       G        LD1          LD2 newG
    1  1 -0.1409984  2.582951755    1
    2  1 -2.3918356  0.825366275    1
    3  1 -0.3704452  1.641514840    1
    4  1 -0.9714835  0.548448277    1
    5  1 -1.7134891  1.246681993    1
    6  2 -2.4593598  1.361571174    1
    7  2  0.3789617 -2.200431689    2
    8  2 -2.5581070 -0.467096091    2
    9  2 -1.1900285 -0.412972027    2
    10 2 -1.7638874 -2.382302324    2
    11 2 -1.1869165 -2.485574940    2
    12 2 -0.1123680 -0.598883922    2
    13 2  0.3399132  0.232863397    3
    14 3  2.8456561  0.936722573    3
    15 3  1.5592346  0.025668216    3
    16 3  0.7457802 -0.209168159    3
    17 3  3.0062824 -0.358989534    3
    18 3  2.2511708  0.008852067    3
    19 3  2.2108260 -0.331206768    3
    20 3  1.5210939  0.035984885    3
    
    > diag(prop.table(tab,1))
       1    2    3 
    1.00 0.75 1.00 
    > 
    > sum(diag(prop.table(tab)))
    [1] 0.9
    
    > predict(ld,data.frame(Q=8,C=7.5,P=65))
    $class
    [1] 2
    Levels: 1 2 3
    
    $posterior
              1        2           3
    1 0.2114514 0.786773 0.001775594
    
    $x
            LD1        LD2
    1 -1.537069 -0.1367865
    

    结果的可视化:

    plot(Z$x)
    text(Z$x[,1], Z$x[,2],G, adj=-0.5,cex=0.75)
    

    在这里插入图片描述

    F2.2 二次判别

    使用F2.1中的数据,判别分析的代码为:

    qd <- qda(G~Q+C+P)
    Z <- predict(qd)
    newG <- Z$class
    cbind(G,newG)
    sum(diag(prop.table(tab)))
    predict(qd,data.frame(Q=8,C=7.5,P=65))
    

    结果:

    > cbind(G,newG)
          G newG
     [1,] 1    1
     [2,] 1    1
     [3,] 1    1
     [4,] 1    1
     [5,] 1    1
     [6,] 2    2
     [7,] 2    2
     [8,] 2    2
     [9,] 2    2
    [10,] 2    2
    [11,] 2    2
    [12,] 2    2
    [13,] 2    3
    [14,] 3    3
    [15,] 3    3
    [16,] 3    3
    [17,] 3    3
    [18,] 3    3
    [19,] 3    3
    [20,] 3    3
    > sum(diag(prop.table(tab)))
    [1] 0.9
    > predict(qd,data.frame(Q=8,C=7.5,P=65))
    $class
    [1] 2
    Levels: 1 2 3
    
    $posterior
                1         2            3
    1 0.008221165 0.9915392 0.0002396287
    

    F2.3 贝叶斯判别

    使用F2.1中的数据,由于在协方差阵相同的情况下,贝叶斯判别的判别函数具有与距离判别判别函数相似的形式,差别在于多了先验概率和误判代价,因此仍然可以使用线性判别分析的函数做贝叶斯判别。不考虑误判代价,则判别分析的代码为:

    ld1 <- lda(G~Q+C+P,prior=rep(1,3)/3)   # 先验概率相等的贝叶斯模型
    ld1
    
    ld2 <- lda(G~Q+C+P,prior=c(5,8,7)/20)   # 先验概率不等的贝叶斯模型
    ld2
    
    Z1 <- predict(ld1)
    cbind(G,Z2$x,newG=Z$class)
    table(G,Z1$class)
    round(Z1$post,3)   # 保留三位有效数字
    
    predict(ld1,data.frame(Q=8,C=7.5,P=65))
    predict(ld2,data.frame(Q=8,C=7.5,P=65))
    

    结果:

    > ld1
    Call:
    lda(G ~ Q + C + P, prior = rep(1, 3)/3)
    
    Prior probabilities of groups:
            1         2         3 
    0.3333333 0.3333333 0.3333333 
    
    Group means:
             Q        C      P
    1 8.400000 5.900000 48.200
    2 7.712500 7.250000 69.875
    3 5.957143 3.714286 34.000
    
    Coefficients of linear discriminants:
              LD1         LD2
    Q -0.92307369  0.76708185
    C -0.65222524  0.11482179
    P  0.02743244 -0.08484154
    
    Proportion of trace:
       LD1    LD2 
    0.7259 0.2741 
    
    > ld2
    Call:
    lda(G ~ Q + C + P, prior = c(5, 8, 7)/20)
    
    Prior probabilities of groups:
       1    2    3 
    0.25 0.40 0.35 
    
    Group means:
             Q        C      P
    1 8.400000 5.900000 48.200
    2 7.712500 7.250000 69.875
    3 5.957143 3.714286 34.000
    
    Coefficients of linear discriminants:
              LD1         LD2
    Q -0.81173396  0.88406311
    C -0.63090549  0.20134565
    P  0.01579385 -0.08775636
    
    Proportion of trace:
       LD1    LD2 
    0.7403 0.2597 
    
    > table(G,Z1$class)
       
    G   1 2 3
      1 5 0 0
      2 1 6 1
      3 0 0 7
      
    > round(Z1$post,3)   # 保留三位有效数字
           1     2     3
    1  0.983 0.006 0.012
    2  0.794 0.206 0.000
    3  0.937 0.043 0.020
    4  0.654 0.337 0.009
    5  0.905 0.094 0.000
    6  0.928 0.072 0.000
    7  0.003 0.863 0.133
    8  0.177 0.822 0.000
    9  0.185 0.811 0.005
    10 0.003 0.997 0.000
    11 0.002 0.997 0.001
    12 0.111 0.780 0.109
    13 0.292 0.325 0.383
    14 0.001 0.000 0.999
    15 0.012 0.023 0.965
    16 0.079 0.243 0.678
    17 0.000 0.000 1.000
    18 0.001 0.003 0.996
    19 0.001 0.004 0.995
    20 0.014 0.025 0.961
    
    > predict(ld1,data.frame(Q=8,C=7.5,P=65))
    $class
    [1] 2
    Levels: 1 2 3
    
    $posterior
             1         2           3
    1 0.300164 0.6980356 0.001800378
    
    $x
            LD1        LD2
    1 -1.426693 -0.5046594
    
    > predict(ld2,data.frame(Q=8,C=7.5,P=65))
    $class
    [1] 2
    Levels: 1 2 3
    
    $posterior
              1        2           3
    1 0.2114514 0.786773 0.001775594
    
    $x
            LD1        LD2
    1 -1.537069 -0.1367865
    

    案例分析

    # ==================判别分析===================== #
    # 两群体同协方差
    path <- "E:\\桌面文档\\学习文件\\大三\\多元统计\\实验\\第一次上机\\第一次上机课的数据\\第一次上机课的数据\\判别分析\\陈炜\\"
    setwd(path)  # 设置工作路径
    
    # ==========================================案例1========================================#
    # =======线性判别分析============= #
    d6.1 <- read.table("6_1.txt", header = TRUE) # 第一行作为列名
    attach(d6.1)
    plot(x1,x2)  # 绘制散点图
    text(x1,x2,G,adj=-0.5)  # 添加文本标记
    detach(d6.1)
    
    library(MASS)
    attach(d6.1)
    ld <- lda(G~x1+x2)   # 做线性判别分析
    ld
    
    Z <- predict(ld)    # 预测
    Z
    newG <- Z$class
    cbind(G,Z$x,newG)
    tab <- table(G, newG)  # 生成原始类别与预测类别的列联表
    sum(diag(prop.table(tab)))   # prop.table将列联表对应元素全部转为相应的概率
    predict(ld, data.frame(x1=8.1,x2=2.0)) # 预测一组新数据对应的输出
    detach(d6.1)
    
    # =========================================案例2===============================#
    # ========二次判别分析============ #
    d6.2 <- read.table("6_2.txt",header = TRUE)
    attach(d6.2)
    plot(Q,C)
    text(Q,C,G,adj=-0.5)
    
    plot(Q,P)
    text(Q,P,G,adj=-0.5)
    
    plot(C,P)
    text(C,P,G,adj=-0.5)
    
    library(MASS)
    qd <- qda(G~Q+C+P)
    qd
    
    Z <- predict(qd) # 预测
    newG <- Z$class
    cbind(G,newG)  # 按列合并
    
    predict(qd, data.frame(Q=8,C=7.5,P=65))   # 预测一组新数据
    
    # =========一次线性判别==================#
    ld <- lda(G~.,data = d6.2)   # 假定协方差相等
    ld
    
    W <- predict(ld)  # 预测值
    cbind(G,W$x,newG=W$class)
    predict(ld,data.frame(Q=8,C=7.5,P=65))
    
    detach(d6.2)
    
    # =========================================案例3===============================#
    ### 多总体距离判别
    ## (等协方差的时候仍然是线性判别)
    d6.3 <- read.table("6_3.txt",header = TRUE)
    attach(d6.3)
    plot(Q,C)
    text(Q,C,G,adj=-0.5)
    
    plot(Q,P)
    text(Q,P,G,adj=-0.5)
    
    plot(C,P)
    text(C,P,G,adj=-0.5)
    
    # -------------------------------------以下部分还没放进word
    ld <- lda(G~Q+C+P)
    ld
    
    Z <- predict(ld)
    newG <- Z$class
    cbind(G,Z$x,newG)
    
    tab <- table(G,newG)
    diag(prop.table(tab,1))
    
    sum(diag(prop.table(tab)))
    
    plot(Z$x)
    text(Z$x[,1], Z$x[,2],G, adj=-0.5,cex=0.75)
    
    predict(ld,data.frame(Q=8,C=7.5,P=65))
    
    ## 二次判别,异方差
    qd <- qda(G~Q+C+P)
    Z <- predict(qd)
    newG <- Z$class
    cbind(G,newG)
    sum(diag(prop.table(tab)))
    predict(qd,data.frame(Q=8,C=7.5,P=65))
    
    ## 贝叶斯判别
    ld1 <- lda(G~Q+C+P,prior=rep(1,3)/3)   # 先验概率相等的贝叶斯模型
    ld1
    
    ld2 <- lda(G~Q+C+P,prior=c(5,8,7)/20)   # 先验概率不等的贝叶斯模型
    ld2
    
    Z1 <- predict(ld1)
    cbind(G,Z2$x,newG=Z$class)
    table(G,Z1$class)
    round(Z1$post,3)   # 保留三位有效数字
    
    predict(ld1,data.frame(Q=8,C=7.5,P=65))
    predict(ld2,data.frame(Q=8,C=7.5,P=65))
    
    detach(d6.3)
    
    ### 案例分析
    Case5 <- read.table('6_case.txt',header = T)
    head(Case5,30)    # 显示前30行数据
    
    attach(Case5)
    plot(Case5[,2:5], gap=0, main="陈炜")    # 绘制两两散点图
    ld <- lda(G~.,data=Case5)
    ld
    
    plot(ld,main="陈炜")
    
    Zld <- predict(ld)
    data.frame(Case5$G,Zld$class,round(Zld$x,3))
    addmargins(table(Case5$G,Zld$class))
    
    qd <- qda(G~.,data=Case5)  # 二次判别
    qd
    
    Zqd <- predict(qd)
    addmargins(table(Case5$G,Zqd$class))
    
    detach(Case5)
    
    ###########################################主成分分析
    path <- "E:\\桌面文档\\学习文件\\大三\\多元统计\\实验\\第一次上机\\第一次上机课的数据\\第一次上机课的数据\\主成分分析\\陈炜\\"
    setwd(path)
    # (1)计算相关矩阵
    X <- read.table("d7_2.txt",header=T)
    cor(X)
    
    # (2)求相关矩阵的特征值和主成分负荷
    PCA <- princomp(X,cor=T)  # cor=T的意思是使用相关矩阵
    PCA
    
    PCA$loadings   # 载荷矩阵
    summary(PCA,loading=TRUE) # 显示载荷矩阵
    
    # (3)确定主成分
    screeplot(PCA,type='lines',main="陈炜")   # 使用直线绘制碎石图
    
    # (4)主成分得分
    PCA$scores[,1:2]   # 其实就是每一个原始变量对应的主成分值
    
    
    ## 案例
    Case8 <- read.table('case8.txt', header = T)
    PC <- princomp(Case8,cor=T)
    summary(PC)
    m <- 2   # 选两个主成分(只有前两个特征值大于1)
    PC$loadings[1:m]
    princomp.rank(PC,m)  # 主成分排序
    princomp.rank(PC,m,plot=T)
    
    展开全文
  • 贝叶斯判别分析

    千次阅读 2020-04-27 09:46:17
    Bayes判别:假定对研究对象已经有一定的认识,但...贝叶斯判别函数如下: discriminiant.bayes<-function (TrnX1, TrnX2, rate=1, TstX = NULL, var.equal = FALSE) #输入X1类样本集,输入X2类样本集,TstX为待测...
  • 贝叶斯判别算法matlab的实现,详细地介绍了贝叶斯算法
  • 判别分析中,至少有一个已经明确知道类别的“训练样本”,利用这个数据,就可以建立判别准则,并通过预测变量来为未知类别的观测值进行判别了。 所谓Fisher判别法,就是一种先投影的方法。 考虑只有两个(预测)...
  • 贝叶斯判别法(SAS实现)

    千次阅读 2020-09-03 17:00:28
    贝叶斯判别法常用于分类问题中 代码: data ex; input g x1-x3 @@; cards; 1 76 99 5374 1 79.5 99 5359 1 78 99 5372 1 72.1 95.9 5242 1 73.8 77.7 5370 2 71.2 93 4250 2 75.3 94.9 3412 2 70 91.2 3390 2 72.8...
  • 贝叶斯判别分析matlab程序

    热门讨论 2008-11-19 16:28:22
    给出总体为正态分布,损失矩阵为0,1的贝叶斯判别分析的matlab程序。
  • 贝叶斯判别

    万次阅读 2017-08-21 21:14:40
    贝叶斯判别法  Naive Bayes 算法基本原理 假设 为一任意样本,它的特征为 其中 表示该样本中出现的第i个特征项。预定义的样本类别为 。假设在给定的条件下,特征项之间都是相互独立的,不存在任何依赖关系。...
  • 作为智能配电网的重要设备之一,电力变压器运行状态的准确判别和...建立基于模糊综合评价和贝叶斯判别的电力变压器状态判别和预警模型。实例证明该模型能准确高效地判别变压器的工作状态并提出相应的预警和检修计划。
  • 从贝叶斯公式到贝叶斯判别准则

    千次阅读 2018-07-26 17:27:20
    原来线性判别分析、平方判别分析、朴素贝叶斯这么简单直白。 前方将出现大量数学公式推导证明,为防止烦躁不适,先复习一下几个重要概念。 1.1一维高斯变量X~N(μ,),则概率密度函数   1.2多维高斯变量 ,X~N(μ...
  • 提出了复杂采空区危险程度辨识的贝叶斯判别方法。基于多元判别分析理论,将贝叶斯判别方法应用于金属矿山采空区危险程度的预测判别问题中,建立了相应的贝叶斯判别分析模型。该模型选用岩石单轴抗压强度、岩石弹性...
  • 机器学习(一)贝叶斯判别式 2018/2/13 by Chenjing Ding 符号 含义 CkCkC_k 第k类 p 概率密度 P(Ck)P(Ck)P(C_k) 第k类的概率。本文中的概率密度和概率在公式推导时已严格区分 x 输入...
  • 煤与瓦斯突出影响因素多,难以为其建立合适的多指标非线性预测模型,为提高突出预测的准确性和增强预测预报方法的实用性,文章选用瓦斯压力、瓦斯放散速度等多项影响判别因子,建立了基于贝叶斯判别分析煤与瓦斯突出的...
  • 利用统计软件SPSS的逐步判别分析功能筛选出判别无烟煤、烟煤、褐煤的主要指标——氢含量和氧含量,以该指标为变量建立贝叶斯逐步线性判别函数,并采用该函数分别对建模样本和测试样本进行识别。识别结果显示:基于SPSS...
  • 针对目前深水油藏分类评价研究现状的不足,基于模糊C均值聚类算法和贝叶斯判别函数,建立了深水油藏指标选择标准和分类评价体系。优选世界三大深水油气区19例油田的特征属性参数作为典型样品集,采用模糊聚类分析对深水...
  • 利用贝叶斯判别函数设计分类器

    千次阅读 2017-03-19 20:48:23
    (1)假设P(w1)=P(w2)=0.5,P(w3)=0,仅利用x1特征值为这两类判别设计一个分类器,并确定样本的经验训练误差,即误分点的百分比。 (2)假设P(w1)=P(w2)=0.5,P(w3)=0,利用x1,x2两个特征值为这两类判别...
  • 介绍了贝叶斯判别和逐步判别法的基本原理,分析了目前出现的一些特征变量优化方法,以油气解释评价中的贝叶斯判别应用为例,对于逐步贝叶斯判别中的变量优化方法进行了研究和总结,提出了变量的多步优化策略和分步多...
  • 包括示例案例
  • 结合主成分分析和贝叶斯(Bayes)判别简化构建突水水源识别模型,水样变量因子选取Ca2+、Na++K+、Mg2+、HCO3-、Cl-、SO42-六个指标。采用潘二矿新生界松散层、煤系砂岩以及太原组灰岩中的水质分析资料作为训练样本和...
  • 首先引入稀疏贝叶斯学习机对上下文模型中图像特 征的后验概率进行建模, 然后通过贝叶斯原理将稀疏贝叶斯模型与隐马尔可夫模型结合, 提出一种能够实现上下文 场景识别模型的判别学习方法. 在真实场景数据库上...
  • 朴素贝叶斯判别法--Matlab

    千次阅读 2016-09-10 16:01:25
    load fisheriris %这是matlab自带的数据,朴素贝叶斯判别的例子 objbayes = NaiveBayes.fit(meas,species);%meas是训练样本 species是分类结果 pred = objbayes.predict(meas);%objbayes.predict预测类别归属 ...
  • 贝叶斯判别函数和绝策面,通过对不同情况的讨论和分析,进一步理解应用贝叶斯理论进行分类器设计的理论知识。
  • 作为统计判别问题的模式分类  模式识别的目的就是要确定某一个给定的模式样本属于哪一类。 可以通过对被识别对象的多次观察和测量,构成特征向量,并将其作为某一个判决规则的输入,按此规则来对样本进行分类。在...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 15,923
精华内容 6,369
关键字:

贝叶斯判别