精华内容
下载资源
问答
  • 贝叶斯判别分析

    2020-04-27 09:46:17
    Bayes判别:假定对研究对象已经有一定的认识,但...贝叶斯判别函数如下: discriminiant.bayes<-function (TrnX1, TrnX2, rate=1, TstX = NULL, var.equal = FALSE) #输入X1类样本集,输入X2类样本集,TstX为待测...

    代码部分借鉴于:https://blog.csdn.net/TOMOCAT/article/details/78723277?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522158795131619725256730493%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fall.%2522%257D&request_id=158795131619725256730493&biz_id=0&utm_source=distribute.pc_search_result.none-task-blog-2allfirst_rank_v2~rank_v25-1

    Bayes判别:假定对研究对象已经有一定的认识,但这种认识常用先验概率来描述,取得样本后,就可以用样本修正已有的先验概率,得到后验概率。

    R语言实现
    回代函数如下:
    discriminiant.bayes<-function
    (TrnX1, TrnX2, rate=1, TstX = NULL, var.equal = FALSE)
    #输入X1类样本集,输入X2类样本集,TstX为待测样本,如果不输入则计算实验样本的回代情况,var.equal=T
    #认为两总体的协方差矩阵相同,.rate= c(1|2) / c(2|1)*p2/p1
    {

    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)
    x1<-nrow(TrnX1)
    x2<-nrow(TrnX2)
    nx<-nrow(TstX)
    mu1<-colMeans(TrnX1)
    mu2<-colMeans(TrnX2)
    S1<-var(TrnX1)
    S2<-var(TrnX2)
    blong<-matrix(rep(0, nx), nrow=1, byrow=TRUE,
    dimnames=list(“blong”, 1:nx))
    if (var.equal == TRUE || var.equal == T)
    {
    S<-((x1-1)S1+(x2-1)S2)/(x1+x2-2)
    #S<-cov(rbind(TrnX1,TrnX2))
    #S<-var(rbind(TrnX1,TrnX2)) 两个S的值相等
    #var(rbind(dat[1:12,2:8],dat[12:35,2:8]))
    beta<-2
    log(rate)
    w<-mahalanobis(TstX, mu2, S)-mahalanobis(TstX, mu1, S)
    }
    else
    {
    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
    }
    交叉确认估计函数:
    Bayes.jiaocha<-function(TrnX1,TrnX2,rate=1,var.equal = FALSE)
    {
    if (is.matrix(TrnX1) != TRUE) TrnX1<-as.matrix(TrnX1)
    if (is.matrix(TrnX2) != TRUE) TrnX2<-as.matrix(TrnX2)
    x<-nrow(TrnX1)+nrow(TrnX2)
    blong1<-matrix(rep(0,x),nrow = 1,byrow = T,dimnames = list(‘blong’,1:x))
    blong2<-matrix(rep(0,x),nrow = 1,byrow = T,dimnames = list(‘blong’,1:x))
    if (rate==1)#先验概率相同
    {
    rate1<-1
    rate2<-1
    }
    else#按比例分配
    {
    rate1<-23/11
    rate2<-22/12
    }
    if (var.equal == TRUE || var.equal == T)#假设协方差相等的时
    {
    for (i in 1:x)
    {
    if(i<=nrow(TrnX1))#需要去除TrnX1中的一???,nrowTrnX1=12
    {
    blong1[i]<-discriminiant.bayes(TrnX1[-i,],TrnX2,TstX = TrnX1[i,],rate=rate1,var.equal = T)
    }
    else
    {
    blong1[i]<-discriminiant.bayes(TrnX1,TrnX2[-(i-12),],TstX = TrnX2[i-12,],rate=rate2,var.equal = T)
    }
    }
    blong1
    }
    else#假设协方差不相等的时
    {
    for (j in 1:x)
    {
    if(i<=nrow(TrnX1))#需要去除TrnX1中的一???,nrowTrnX1=12
    {
    blong2[j]<-discriminiant.bayes(TrnX1[-j,],TrnX2,TstX = TrnX1[j,],rate=rate1,var.equal = F)
    }
    else
    {
    blong2[j]<-discriminiant.bayes(TrnX1,TrnX2[-(j-12),],TstX = TrnX2[j-12,],rate=rate2,var.equal = F)
    }
    blong2
    }
    }

    }

    展开全文
  • 煤与瓦斯突出影响因素,难以为其建立合适的指标非线性预测模型,为提高突出预测的准确性和增强预测预报方法的实用性,文章选用瓦斯压力、瓦斯放散速度等项影响判别因子,建立了基于贝叶斯判别分析煤与瓦斯突出的...
  • 贝叶斯判别分析matlab程序

    热门讨论 2008-11-19 16:28:22
    给出总体为正态分布,损失矩阵为0,1的贝叶斯判别分析的matlab程序。
  • 贝叶斯判别法是基于平均错判损失最小的前提下,定义的判别法。

    十四、贝叶斯判别法

    1.贝叶斯判别的定义

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

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

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

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

    先验概率代表了出现类别的概率分布,这是在没有任何样本信息时能做出的关于类的直接判断。假设来自第ii类的先验概率是qiq_i,那么此时的错判平均损失,实际上是一种关于先验概率的加权平均。现在,我们可以定义判别法DD的错判平均损失为
    g(D)=i=1kqij=1kP(ji)L(ji)=defi=1kqiri(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).
    这样,贝叶斯判别准则就可以被视为这样的最优化问题:找到一个DD^*,使得g(D)=minDg(D)g(D^*)=\min_Dg(D)

    2.贝叶斯判别的解

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

    据此,我们可以先给出错判概率为P(ji)P(j|i),它表示样品XX本身来自密度函数fi(X)f_i(X),但落在区域DjD_j内:
    P(ji;D)=Djfi(X)dX. P(j|i;D)=\int_{D_j}f_i(X){\rm d}X.
    所以此时的错判平均损失是
    g(D)=i=1kqij=1kL(ji)Djfi(X)dX=j=1kDji=1kqiL(ji)fi(X)dX=dj=1kDjhj(X)dX. \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}
    这里定义
    hj(X)=defi=1kqiL(ji)fi(X), h_j(X)\stackrel {\rm def}= \sum_{i=1}^k q_iL(j|i)f_i(X),
    它表示把样品XX归到GjG_j类的平均损失,注意到hj(X)h_j(X)DD无关,对hj(X)h_j(X)求和,就得到了错判平均损失。对于贝叶斯判别的解DD^*,要使得g(D)g(D^*)是所有DD中最小的,所以
    g(D)g(D)=i=1kDihi(X)dXj=1kDjhj(X)dX=i=1kj=1kDiDj[hi(X)hj(X)]dX0. \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}
    由此,我们能够得到贝叶斯判别的解是:在所有hi(X)h_i(X)中,如果hj(X)h_j(X)最小,就将XX判别为第jj类,在这个判别法的条件下,能够让g(D)g(D)0g(D^*)-g(D)\le 0恒成立。

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

    3.广义马氏距离

    对于正态总体,在错判损失都相等的情况下,
    qifi(X)=qi(2π)m/2Σ1/2exp{12(Xμ(i))Σ1(Xμ(i))},lnqifi(X)=C+lnqi12lnΣ12(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)}).
    因此,我们定义样本XX到总体GiG_i的广义马氏距离为
    Di2(X)=di2(X)+lnS2lnqi. D_i^2(X)=d^2_i(X)+\ln |S|-2\ln q_i.
    可以看到,当样本XX到总体GiG_i的广义马氏距离最小的时候,它会被归类到GiG_i。因此,在每一类都是多元正态总体,且错判损失相等的情况下,用广义马氏距离替代马氏距离,贝叶斯判别的解与直接判别法是一样的。

    回顾总结

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

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

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

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

    展开全文
  • Spark 分布式实现贝叶斯判别 贝叶斯公式 假设事件B1,B2...Bn是样本空间Ω的一个分割,且他们各自的概率为P(B1),P(B2),P(Bn)假设事件 B_1,B_2...B_n 是样本空间Ω的一个分割,且他们各自的概率为P(B_1),P(B_2),P(B_n...

    Spark 分布式实现贝叶斯判别

    贝叶斯公式

    贝叶斯
    B1,B2...BnΩP(B1),P(B2),P(Bn)假设事件 B_1,B_2...B_n 是样本空间Ω的一个分割,且他们各自的概率为P(B_1),P(B_2),P(B_n)
    A是事件Ω中的一个事件,则在A给定的条件下,事件Bi的条件概率如下:
    贝叶斯公式
    Bi 通常视为A发生的”原因“,P(Bi)称为先验概率(主观概率),表示各种原因发生的可能性大小;P(Bi|A)(i=1,2…)则反映当出现结果A之后,再对各种“原因”概率的新认识,故称后验概率。

    一个小案例

    大家都知道狼来了的故事,我们就有贝叶斯的思想来解释一下这个故事:
    在最开始的时候,大家对于放羊娃的认识不够深刻,主观上认为放羊娃说真话(记为事件B1)和说假话(记为事件B2)的概率相同。即:
    P(B1)=P(B2)=0.5P(B_1) = P(B_2) = 0.5
    再假设狼来了(记为事件A),说谎话喊狼来了时,狼来的概率为1/3,说真话喊狼来了时,狼来的概率是2/3,即:
    P(AB1)=1/3;P(AB2)=2/3P(A|B_1) = 1/3 ; P(A|B_2) =2/3

    (A),第一次村民上山打狼,狼没来(记为事件\overline{A}),此时村民们对放羊娃就有了新的认识:
    狼没来的情况下小孩说谎了(在村民们的主观印象上,小孩说谎的概率增加了):
    P(B1A)=P(AB1)P(B1)k=1NP(ABi)P(Bi)=811P(B_1|\overline{A})=\frac{P(\overline{A}|B_1)P(B_1)}{\sum_{k=1}^N P(\overline{A}|B_i)P(B_i)} = \frac{8}{11}

    随着小孩说谎的次数增加,村民们对于小孩说谎的主观概率也不断增加,当这个概率增加到一定程度时计算小孩说真话,村民们就不会再相信他。

    贝叶斯判别

    之前提到的两种判别分析方法都非常简单,实用,但是也存在着一定的缺点:一是判别方法与各个总体出现的概率大小无关,而是与错判后造成的损失无关。贝叶斯判别则考虑了这两种情况:贝叶斯判别假定对样本有一定的认知(先验概率),然后计算得出后验概率并进行统计推断的判别方法。

    贝叶斯判别求解过程

    由于过程较长,公式比较多,直接上书
    在这里插入图片描述在这里插入图片描述
    在这里插入图片描述

    代码实现过程

    本案例使用的数据为鸢尾花数据集

    def main(args: Array[String]): Unit = {
    
        val spark = SparkSession
          .builder()
          .appName(s"${this.getClass.getSimpleName}")
          .master("local[*]")
          .getOrCreate()
    
        import spark.implicits._
    
        val sc = spark.sparkContext
    
        val irisData = spark.read
          .option("header", true)
          .option("inferSchema", true)
          .csv("F:\\DataSource\\iris.csv")
    
        val schema = irisData.schema
        val fts = schema.filterNot(_.name == "class").map(_.name).toArray
    
        val amountVectorAssembler: VectorAssembler = new VectorAssembler()
          .setInputCols(fts)
          .setOutputCol("features")
    
        val vec2Array = udf((vec: DenseVector) => vec.toArray)
    
        val irisFeatrus = amountVectorAssembler
          .transform(irisData)
          .select($"class", vec2Array($"features") as "features")
    
        val p: Long = irisFeatrus.count()
    
       // 计算均值向量的自定义聚合函数(请参考之前的两篇文章)
        val ui = spark.udf.register("udafMedian", new meanVector(fts.length))
    
        // 计算样本均值向量
        val uiGroup = irisFeatrus
          .groupBy($"class")
          .agg(ui($"features") as "ui", count($"class") as "len")
    
        val covMatrix = irisFeatrus
          .join(uiGroup, "class")
          .rdd
          .map(row => {
            val lable = row.getAs[String]("class")
            val len = row.getAs[Long]("len")
            val u = densevec(row.getAs[Seq[Double]]("ui").toArray)
            val x = densevec(row.getAs[Seq[Double]]("features").toArray)
            val denseMatrix = (x - u).toDenseMatrix
            lable -> (denseMatrix, u, len)
          })
          .reduceByKey((d1, d2) => {
            // 矩阵合并,均值向量,样本大小
            (DenseMatrix.vertcat(d1._1, d2._1), d1._2, d1._3)
          })
          .mapValues(tp => {
            val covm: DenseMatrix[Double] =
              (tp._1.t * tp._1).map(_ / (tp._3 - 1)) //协方差矩阵
            val qi = math.log(tp._3.toDouble / p) // 先验概率,在此默认为各类样本的频率
              (covm, tp._2.toDenseMatrix, qi)
          })
    
    
        val covBroad = sc.broadcast(covMatrix.collect())
        val predictudf = udf((seq: Seq[Double]) => {
          val dist: Array[(String, Double)] = covBroad.value
            .map(tp => {
            /**
            * 计算判别函数的相关指标
            **/
              val xi = densevec(seq.toArray).toDenseMatrix
              val inCov = inv(tp._2._1)
              val lnCov = math.log(det(tp._2._1)) / 2
              val xdiff = (xi * inCov * xi.t).data.head / 2
              val mdist = (tp._2._2 * inv(tp._2._1) * tp._2._2.t).data.head / 2
              val xu = (xi * inCov * tp._2._2.t).data.head
              val d = tp._2._3 - lnCov - xdiff - mdist + xu
              (tp._1, d)
            })
    
          val pm = dist.map(x => math.exp(x._2)).sum
    
         // 计算后验概率
          dist.map(tp => {
            tp._1 -> math.exp(tp._2) / pm
          })
          
        })
    
        irisFeatrus
          .withColumn("prediction", predictudf($"features"))
          .show(truncate = false)
    
        spark.stop()
      }
    

    结果查看:从结果看出,分类效果还是很好的

    |class      |features            |prediction |
    +-----------+--------------------+-----------+
    |Iris-setosa|[5.1, 3.5, 1.4, 0.2]|Iris-setosa|
    |Iris-setosa|[4.9, 3.0, 1.4, 0.2]|Iris-setosa|
    |Iris-setosa|[4.7, 3.2, 1.3, 0.2]|Iris-setosa|
    

    由于作者水平有限,在介绍及实现过程中难免有纰漏之处,欢迎细心的朋友指正

    参考资料:

    《多元统计分析及R语言建模》–王斌会
    《概率论与数理统计》 --茆师松

    展开全文
  • 判别分析是用以判别个体所属群体的一种统计方法,它产生于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类训练样本,其输入格式是数据框,或矩阵(样本按行输入). rate=L(12)L(21)p2p1rate=\frac{L(1|2)}{L(2|1)}\cdot \frac{p_2}{p_1},缺省值为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)
    
    展开全文
  • 贝叶斯判别:对象(总体)在抽样前已有一定的认识,常用先验分布来描述这种认识,然后给予抽取的样本再对先验认识作修正,得到后验分布,而各种统计推断均基于后验分布进行。将Bayes 统计的思想用于判别分析,就得到...
  • 断它来自哪个总体贝叶斯(BAYES) 每个人脑中都有一个贝叶斯 判别思想是根据先验概率求出后验概率,并依据后验概率分布作出统计推断。 所谓先验概率,就是用概率来描述人们对所研究的对象的认识的程度;
  • 两类(M=2)情况的贝叶斯最小风险判别 选M=2,即全部的模式样本只有ω1和ω2两类,要求分类器将模式样本分到ω1和ω2两类中,则平均风险可写成: 当分类器将x判别为ω1时: 当分类器将x判别为ω2时: 若r1(x...
  • 贝叶斯和全概率的区分和做题判别

    千次阅读 2019-06-23 14:44:31
    前面是问总体看来被偷的概率是多少,现在是知道了总体被偷了这件事,概率并不知道,问你个更有意思的问题,像是侦探断案:是哪个小偷的偷的,计算每个小偷偷的概率。 这个特性用在机器学习,人工智能领域相当好用。...
  • 到目前为止,我们主要...Part4 生成模型、高斯判别分析、朴素贝叶斯1.判别学习算法和生成学习算法① 判别学习算法(discriminative learning algorithm):训练出一个总模型,把新来的样本放到这个总模型中,直接...
  • 多元高斯密度下的判别函数:线性判别函数LDF、二次判别函数QDF (三)贝叶斯错误率 (四)生成式模型的参数估计:贝叶斯学派与频率学派;极大似然估计、最大后验概率估计、贝叶斯估计;多元高斯密度下的参数估计 ...
  • R语言 判别分析

    千次阅读 2016-12-27 10:48:29
    #贝叶斯判别 贝叶斯判别式假定对研究对象已有一定的认识 这种认识常用先验概率来描述 #当取得样本后 就可以用样本来修正已经有的先验概率分布 得出后验概率分布 #然后通过后验概率分布 进行各种统计推断 #实际上...
  • 贝叶斯

    2009-03-31 20:52:00
    贝叶斯【理论概述】 贝叶斯 Thomas Bayes,英国数学家.1702年出生于伦敦,做过神甫。1742年成为英国皇家学会会员。1763年4月7日逝世。贝叶斯在数学方面主要研究概率论。他首先将归纳推理法用于概率论基础理论,并...
  • 朴素贝叶斯

    2018-09-10 14:23:37
    本文章参考: 周志华 《机器学习》 李航《统计学习方法》 很多人谈朴素贝叶斯都喜欢从它每个属性独立的特点来...贝叶斯决策论的目的是最小化总体风险。我们面对一个分类问题,假设可能的类别有n种,针对一个...
  • 贝叶斯 - 《贝叶斯统计》笔记

    万次阅读 多人点赞 2017-04-21 17:13:28
    贝叶斯统计 - 茆诗松》 茆诗松《贝叶斯统计》目前看过的讲贝叶斯方法最通俗易懂的书了 下载了在这里 第一章 先验分布和后验分布 1.1 三种信息  统计学的两个主要学派:频率学派,贝叶斯学派  统计推断的...
  • 多元统计:判别分析

    2020-06-09 13:11:26
    贝叶斯判别法1. 贝叶斯判别法原理2. 贝叶斯判别法则3. 贝叶斯判别法例题四. 费希尔判别法1. 费希尔判别法原理五. 距离判别法、贝叶斯判别法和费希尔判别法的异同 一. 判别分析介绍 判别分析,是一种统计判别和分组...
  • 最小错误率贝叶斯决策和最小风险贝叶斯决策matlab代码,贝叶斯决策理论方法,是统计模式识别中的一个基本方法。贝叶斯决策判据,既考虑了各类参考总体出现的概率大小,又考虑了因误判造成的损失大小,判别能力强。
  • 贝叶斯贝叶斯公式

    2019-08-04 11:22:51
    人机与认知实验室人机与认知实验室 Table of Contents 贝叶斯公式 理论分析 ...【为了证明上帝的存在,他发明了概率统计学原理,遗憾的是,他的这一美好愿望至死也未能实现】 ...贝叶斯(约1701-1761) Thomas Baye...
  • Part IV Generative Learning Algorithms 回顾上一部分的内容,我们解决问题的出发点在于直接对p(y|x;...这样建模的算法,称为判别学习算法(discriminative learning algorithms)。除此之外,有另...
  • 贝叶斯分类

    2018-09-05 19:37:59
    前言 ...贝叶斯判定准则的大意是说,对于一个分类问题,我们的目的是要找到一个判定准则hhh,使得总体风险最小化(这里的风险是指一个属于i的样本被误分到j类中)。而为使总体风险最小化,...
  • 贝叶斯算法

    2018-05-25 00:04:26
    谈到贝叶斯,会提到概率论的两大学派,频率学派和贝叶斯学派,也是机器学习中的判别方法和生成方法,对于大多数的分类算法,如决策树,SVM,逻辑回归,KNN等,这些都是判别方法,也就是直接学习出特征输出Y和特征X...
  • 贝叶斯决策论

    2021-05-20 05:00:34
    文章目录贝叶斯决策论先验概率似然概率最小错误率贝叶斯贝叶斯公式:最大后验分类规则贝叶斯判定准则最小风险贝叶斯最小风险判决步骤最小风险贝叶斯的目标贝叶斯判定准则困难与策略判别式模型生成式模型Example极大...
  • 1. 贝叶斯分类 贝叶斯策论是概率框架下实施决策的基本方法,适用于各种分类任务; 简单聊一下 先验 和 后验 概率: 条件: 今天风大/小(用A表示该事件) 结果:是否可以骑行(用B表示该事件) 那么: 先验概率就是 P(A)P(A)...
  • 贝叶斯统计

    2019-09-13 11:47:50
    20世纪仍然是频率统计占领上风,但很学者预测甚至迹象表明,21世纪将是贝叶斯统计的天下[1]。 1. 贝叶斯定理 不像频率统计只根据当前抽样的结果进行推断,贝叶斯统计还需要考虑经验(pri...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 2,797
精华内容 1,118
关键字:

多总体贝叶斯判别