精华内容
下载资源
问答
  • 为了评估二元经济结果之间的因果关系,我们考虑对面板数据的双变量动态概率模型进行估计,该模型具有特殊性,可以说明动态过程的初始条件。 由于似然函数的难解形式是二维积分,因此我们使用一种近似方法:自适应...
  • 本文考虑险种二项风险模型 ,对保单到达时收取的保费是一随机变量进行了研究 ,得到了其破产 概率的一般公式和lundberg不等式.
  • 为有效计及源荷侧响应对潮流分布的影响,基于系统随机注入量与响应量的概率模型,从可调度常规机组和柔性负荷的不平衡功率基准值的分配、响应随机变量的计算2个方面对基于半不变量法的概率潮流模型进行了改进,并...
  • 在前三周的作业中,我构造了概率模型并调用第三方的求解器对器进行了求解,最终获得了每个随机变量的分布(有向图),最大后验分布(双向图)。本周作业的主要内容就是自行编写概率模型的求解器。实际上,从根本...

      在前三周的作业中,我构造了概率图模型并调用第三方的求解器对器进行了求解,最终获得了每个随机变量的分布(有向图),最大后验分布(双向图)。本周作业的主要内容就是自行编写概率图模型的求解器。实际上,从根本上来说求解器并不是必要的。其作用只是求取边缘分布或者MAP,在得到联合CPD后,寻找联合CPD的最大值即可获得MAP,对每个变量进行边缘分布求取即可获得边缘分布。但是,这种简单粗暴的方法效率极其低下,对于MAP求取而言,每次得到新的evidance时都要重新搜索CPD,对于单个变量分布而言,更是对每个变量都要反复进行整体边缘化。以一个长度为6字母的单词为例,联合CPD有着多达26^6个数据,反复操作会浪费大量的计算资源。

    1、团树算法初始化

      团树算法背后的思路是分而治之。对于一组随机变量ABCDEFG,如果A和其他变量之间是独立的,那么无论是求边缘分布还是MAP都可以将A单独考虑。如果ABC联系比较紧密,CDE联系比较紧密,那么如果两个团关于C的边缘分布是相同的,则我们没有必要将ABCDE全部乘在一起再来分别求各个变量的边缘分布。因为反过来想,乘的时候也只是把对应的C乘起来,如果C的边缘分布相同,在相乘的时候其实两个团之间并没有引入其他信息,此时乘法不会对ABDE的边缘分布产生影响。团树算法的数学过程和Variable Elimination是相同的。

      PGM在计算机中的表达是factorLists,factor的var(i),var表示节点连接关系。val描述了factor中var的关系。cliqueTree其实是一种特殊的factorLists,它的var是clique,表示一堆聚类的var。它的val表示的还是var之间的关系。只不过此时var之间的连接不复存在了。所以clique由两个变量组成:1、cliqueTree 2、edges.

      团树算法的初始化可以分为两个过程:1、将变量抱团;2、获取团的初始势;

      变量抱团是一个玄学过程,因为有很多不同的抱法,而且还都是对的。比较常见的是最小边,最小割等...其实如果是人来判断很容易就能得到结果,但是使用计算机算法则要费一些功夫了。不过这不涉及我们对团树算法的理解,所以Koller教授代劳了。

      团的初始势表示团里变量之间的关系。其算法如下,需要注意的是不能重复使用factor.因为一个factor表达了一种关系,如果两个团里都有同一个factor,那么就是...这个事情。。。你帮他重复一遍。。。等于你也有责任的,晓得吧?

      

     1 %COMPUTEINITIALPOTENTIALS Sets up the cliques in the clique tree that is
     2 %passed in as a parameter.
     3 %
     4 %   P = COMPUTEINITIALPOTENTIALS(C) Takes the clique tree skeleton C which is a
     5 %   struct with three fields:
     6 %   - nodes: cell array representing the cliques in the tree.
     7 %   - edges: represents the adjacency matrix of the tree.
     8 %   - factorList: represents the list of factors that were used to build
     9 %   the tree. 
    10 %   
    11 %   It returns the standard form of a clique tree P that we will use through 
    12 %   the rest of the assigment. P is struct with two fields:
    13 %   - cliqueList: represents an array of cliques with appropriate factors 
    14 %   from factorList assigned to each clique. Where the .val of each clique
    15 %   is initialized to the initial potential of that clique.
    16 %   - edges: represents the adjacency matrix of the tree. 
    17 %
    18 % Copyright (C) Daphne Koller, Stanford University, 2012
    19 
    20 
    21 
    22 function P = ComputeInitialPotentials(C)
    23 Input = C;
    24 % number of cliques
    25 N = length(Input.nodes);
    26 
    27 % initialize cluster potentials 
    28 P.cliqueList = repmat(struct('var', [], 'card', [], 'val', []), N, 1);
    29 P.edges = zeros(N);
    30 
    31 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    32 % YOUR CODE HERE
    33 %
    34 % First, compute an assignment of factors from factorList to cliques. 
    35 % Then use that assignment to initialize the cliques in cliqueList to 
    36 % their initial potentials. 
    37 
    38 % C.nodes is a list of cliques.
    39 % So in your code, you should start with: P.cliqueList(i).var = C.nodes{i};
    40 % Print out C to get a better understanding of its structure.
    41 %
    42 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    43 % N_factors = length(C.factorList);
    44 for i = 1:N
    45     k = 1;
    46     clear clique_
    47      N_factors = length(Input.factorList);
    48     for j = 1:N_factors
    49         if min(ismember(Input.factorList(j).var,Input.nodes{i}))
    50             clique_(k) = Input.factorList(j);
    51             k = k+1;
    52             Input.factorList(j) =struct('var', [], 'card', [], 'val', []);
    53         end
    54     end
    55     Joint_Dis_cliq = ComputeJointDistribution(clique_);
    56     Joint_Dis_cliq_std = StandardizeFactors(Joint_Dis_cliq);
    57     P.cliqueList(i) = Joint_Dis_cliq_std;
    58 end
    59 P.edges = Input.edges;
    60 end
    View Code

    2、团树的校准

      继续之前的例子,ABC联系比较紧密,CDE联系比较紧密,所以抱成了两个团。如果其关于C的边缘分布相同,那么我们则可以在直接对两个团求ABDE的边缘分布,而不用乘起来了。然而令人悲伤的是现实中往往C的边缘分布是不同的。这时就需要对团树进行校准,希望经过“校准”这个操作后,两边关于C达成了一致意见。显然,一棵校准后的团树求任意一个变量的边缘分布都是方便的,只要对很小规模的联合分布进行边际化就行。

      要使得两边关于C的意见达成一致,最简单的方法就是把C在“A团”中的边缘分布乘以"E团”的势。反过来再把A在“E团”中的边缘分布乘以A团的势。那么此时C在两个团中的边缘分布就完全一样了 all = margin(C,A)*margin(C,E)。此即为团树校准的朴素想法。在数学上,团树的校准依然来自VE算法。让AB领盒饭后,C继续参加下一轮的VE。AB领盒饭剩下的C就是C在A团中的边缘分布。

      团树校准的关键是知道消息传播的顺序。消息一般先由叶向根传递,再由根向叶传递。并且,一个团在得到其所有邻团的消息之前,不能向下一个团传递消息。消息传递顺序获取算法如下:

      

     1 %COMPUTEINITIALPOTENTIALS Sets up the cliques in the clique tree that is
     2 %passed in as a parameter.
     3 %
     4 %   P = COMPUTEINITIALPOTENTIALS(C) Takes the clique tree skeleton C which is a
     5 %   struct with three fields:
     6 %   - nodes: cell array representing the cliques in the tree.
     7 %   - edges: represents the adjacency matrix of the tree.
     8 %   - factorList: represents the list of factors that were used to build
     9 %   the tree. 
    10 %   
    11 %   It returns the standard form of a clique tree P that we will use through 
    12 %   the rest of the assigment. P is struct with two fields:
    13 %   - cliqueList: represents an array of cliques with appropriate factors 
    14 %   from factorList assigned to each clique. Where the .val of each clique
    15 %   is initialized to the initial potential of that clique.
    16 %   - edges: represents the adjacency matrix of the tree. 
    17 %
    18 % Copyright (C) Daphne Koller, Stanford University, 2012
    19 
    20 
    21 
    22 function P = ComputeInitialPotentials(C)
    23 Input = C;
    24 % number of cliques
    25 N = length(Input.nodes);
    26 
    27 % initialize cluster potentials 
    28 P.cliqueList = repmat(struct('var', [], 'card', [], 'val', []), N, 1);
    29 P.edges = zeros(N);
    30 
    31 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    32 % YOUR CODE HERE
    33 %
    34 % First, compute an assignment of factors from factorList to cliques. 
    35 % Then use that assignment to initialize the cliques in cliqueList to 
    36 % their initial potentials. 
    37 
    38 % C.nodes is a list of cliques.
    39 % So in your code, you should start with: P.cliqueList(i).var = C.nodes{i};
    40 % Print out C to get a better understanding of its structure.
    41 %
    42 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    43 % N_factors = length(C.factorList);
    44 for i = 1:N
    45     k = 1;
    46     clear clique_
    47      N_factors = length(Input.factorList);
    48     for j = 1:N_factors
    49         if min(ismember(Input.factorList(j).var,Input.nodes{i}))
    50             clique_(k) = Input.factorList(j);
    51             k = k+1;
    52             Input.factorList(j) =struct('var', [], 'card', [], 'val', []);
    53         end
    54     end
    55     Joint_Dis_cliq = ComputeJointDistribution(clique_);
    56     Joint_Dis_cliq_std = StandardizeFactors(Joint_Dis_cliq);
    57     P.cliqueList(i) = Joint_Dis_cliq_std;
    58 end
    59 P.edges = Input.edges;
    60 end
    View Code

      

      在获取消息传递顺序之后,则可进一步对被传递的消息进行计算。被传递的消息应为某个团对被传播变量的“所有认知”,所有认知则包括该团本身对该消息的认知,以及该团收到的“情报”。需要注意的是,向下家报告情报的时候要对所有信息进行总结,但是不能将下家告诉你的事情重复一遍。因为。。。重复一遍你也有责任的,知道吧。。。。

      

     1     while (1)
     2 
     3       [i,j]=GetNextCliques(P,MESSAGES);
     4 
     5       if i == 0
     6         break
     7       end
     8 
     9       to_be_summed =  setdiff(P.cliqueList(i).var,P.cliqueList(j).var);
    10       to_be_propogan  =  setdiff(P.cliqueList(i).var,to_be_summed);
    11 
    12       tmp_ = 1;
    13       clear factorList
    14       for k = 1:N
    15           if P.edges(i,k)==1&&k~=j&&~isempty(MESSAGES(k,i).var)
    16               factorList(tmp_) = MESSAGES(k,i); 
    17               tmp_ = tmp_+1;
    18           end
    19       end
    20       factorList(tmp_) = P.cliqueList(i);
    21       MESSAGES(i,j) = ComputeMarginal(to_be_propogan,ComputeJointDistribution(factorList),[]);
    22     end
    View Code

      在消息完成从顶向下以及从下到上的传播后,每个团需要根据周边传来的消息进行总结。也就是把消息与本身的势相乘(消息是一种边缘分布)

     1 N = length(P.cliqueList);
     2     for i = 1:N
     3         tmp_ = 1;
     4         for k = 1:N
     5           if P.edges(i,k)==1
     6               factorList(tmp_) = MESSAGES(k,i); 
     7               tmp_ = tmp_+1;
     8           end
     9         end
    10         factorList(tmp_) = P.cliqueList(i);
    11         belief(i) = ComputeJointDistribution(factorList);
    12         clear factorList
    13     end
    View Code

      此时,团树称为已经校准。对各个团的中的变量进行marginal就可以得到每个变量的边缘分布了。

    3、基于团树的MAP估计

      在很多时候,我们可能对单个变量的分布并不感兴趣,而是对[ABCDE]这个组合取哪个值概率最大感兴趣。这个思想可以用于信号解码,OCR,图像处理等领域。很多时候我们不关心单个像素的label是啥,只关心分割出来的像素块label是啥。这类问题称为最大后验估计(MAP)。

                argmaxP(AB)  = argmaxP(A)P(B|A) = argmax_a{P(A){argmax_bP(B|A)}

      显然,从上述过程中,很容易联想到之前提到的边际。只不过这里把边际换成了argmax。P(A){argmax_bP(B|A)}的结果依旧是分布,只不过这个分布的前提是无论A取哪个值,其assignment to val都对应着argmax_b。也就是说,此时如果选择最大的val,那么assignment则对应的是argmax_ab。这种操作的意义就在于可以对一组变量的MAP分而治之,最终单个变量的MAP就是全局MAP的一部分。此时的MESSAGE计算如下:

      

     1 for i = 1:N
     2     P.cliqueList(i).val = log(P.cliqueList(i).val);
     3 end
     4 
     5     while (1)
     6 
     7       [i,j]=GetNextCliques(P,MESSAGES);
     8 
     9       if i == 0
    10         break
    11       end
    12 
    13       to_be_summed =  setdiff(P.cliqueList(i).var,P.cliqueList(j).var);
    14       to_be_propogan  =  setdiff(P.cliqueList(i).var,to_be_summed);
    15 
    16       tmp_ = 1;
    17       clear factorList
    18       for k = 1:N
    19           if P.edges(i,k)==1&&k~=j&&~isempty(MESSAGES(k,i).var)
    20               factorList(tmp_) = MESSAGES(k,i); 
    21               tmp_ = tmp_+1;
    22           end
    23       end
    24       factorList(tmp_) = P.cliqueList(i);
    25       F = factorList;
    26       Joint = F(1);
    27         for l = 2:length(F)
    28             % Iterate through factors and incorporate them into the joint distribution
    29             Joint = FactorSum(Joint, F(l));
    30         end
    31       MESSAGES(i,j) = FactorMaxMarginalization(Joint,to_be_summed);
    32     end
    View Code

        此处对val取对数是因为在map估计时,card一般都比较大。对应的val太小不便于作乘法(OCR的card是26!!!)

       消息的综合如下:

     1     
     2     for i = 1:N
     3         tmp_ = 1;
     4         for k = 1:N
     5           if P.edges(i,k)==1
     6               factorList(tmp_) = MESSAGES(k,i); 
     7               tmp_ = tmp_+1;
     8           end
     9         end    
    10      factorList(tmp_) = P.cliqueList(i);
    11      F = factorList;
    12      belief = F(1);
    13         for l = 2:length(F)
    14             % Iterate through factors and incorporate them into the joint distribution
    15             belief = FactorSum(belief, F(l));
    16         end
    17     
    18      clear factorList
    19      Belief(i) = belief;    
    20     end
    View Code

      

    4、总结

      团树算法作为一种精确推理算法在VE算法的基础上大幅减小了计算量和搜索空间。但其作为一种精确推理方法,依旧有着较大局限性。下周的Homework会以实现MCMC算法为目标~就是Alpha狗用的哪个蒙特卡罗哦~敬请期待。

          所有代码请点这里

      

     

    转载于:https://www.cnblogs.com/ironstark/p/5396791.html

    展开全文
  • 概率销售在质量差异产品寡头市场下的演化博弈研究,杨光,刘新旺,在异质产品市场的国内外寡头销售商垄断销售背景下,本文借助Hotelling模型研究了概率销售和传统销售策略作为决策变量的博弈演化过�
  • 在第二层中,将基于空间马尔可夫随机场的梯度轮廓合并到ELN惩罚中,以消除后验分布的隐性边际概率,以鼓励空间平滑性。 此外,训练样本的真实标签被固定为第二层优化模型中的附加约束,以进一步提高分类准确性。 ...
  • 使用概率模型的原因 k均值等价于假设了球对称形状的聚类。使用带权欧式距离,仍然假设了轴对齐的椭球。...双变量高斯分步,协方差矩阵的主对角线决定了展度;副对角线决定朝向 高斯混合模型 GMM...

    使用概率模型的原因

    k均值等价于假设了球对称形状的聚类。使用带权欧式距离,仍然假设了轴对齐的椭球。没有考虑聚类的形状。
    促使概率模型的原因:混合模型

    • 提供观测点到聚类的软分配soft assignment(分配包含不确定性)
    • 考虑了聚类的形状而不仅仅是中心
    • 允许从不同维度来学习权重

    高斯分布

    这里写图片描述
    双变量高斯分步,协方差矩阵的主对角线决定了展度;副对角线决定朝向
    这里写图片描述

    高斯混合模型

    这里写图片描述
    GMM估计的EM算法
    这里写图片描述

    GMM推导

    这里写图片描述
    这里写图片描述
    这里写图片描述
    这里写图片描述

    EM算法的收敛性,初始化和过拟合

    收敛性

    • EM是一种坐标下降算法,等价于交替最大化E步和M步目标函数
    • 收敛于局部最优解

    初始化

    • 许多初始化EM算法的方式
    • 对于收敛和局部最优解的质量重要
    • 距离
      • 随机选择K个中心
    • 选择类似k-means++
    • 根据k-means的解初始化
    • 通过划分(移除)簇直到形成k个簇生成混合模型

    过拟合

    不要让方差为0.
    协方差矩阵对角线添加小的常数。
    使用贝叶斯方法为参数添加先验。

    k-means与EM关系

    • k-means是两个步骤交替进行,可以分别看成E步和M步;
    • E步中将每个点分给中心距它最近的类(硬分配),可以看成是EM算法中E步(软分配)的近似。当方差无限小的时候,EM相当于k-means。
    • M步中将每类的中心更新为分给该类各点的均值,可以认为是在「各类分布均为单位方差的高斯分布」的假设下,最大化似然值;

    这里写图片描述

    参考资料

    知乎王赟 Maigo的回答
    《Machine Learning》ColumbiaX: CSMM.102x Lecture 16

    展开全文
  • 该方法是在线性回归模型中引入识别变量,借助于双层Bayesian模型和Gibbs抽样算法,给出识别变量后验概率的计算方法和变量选择的方法,通过比较这些识别变量的后验概率进行异常值定位。最后进行了大量的模拟实验,...
  • 利用多变量逻辑函数和变值原理, 在N元0-1输入/输出序对上形成变值测量四元组,建立变值路模拟模型。该模型根据多元/条件概率、同步/异步、对称/反对称等不同条件得到对应概率统计分布,形成4组16个统计直方图。...
  • 在异质产品市场的国内外寡头销售商垄断销售背景下,本文借助Hotelling模型研究了概率销售和传统销售策略作为决策变量的动态博弈演化过程,系统演化的结果表明零售商...
  • 在[0,1]线性城市模型中,两家零售商在第一阶段同时选择价格策略变量,第二阶段确定价格的大小及相应的定价概率,第三阶段由消费者选择零售商.应用演化博弈论进行分析,得到了价格促销策略为寡头零售市场的唯一演化稳定...
  • 如果概率模型变量都是观测变量,那么给定数据,可以直接用极大似然估计法,或贝叶斯估计法来估计模型参数,但是,当模型含有隐变量时,情况就复杂一些,相当于一个双层的概率模型,要估计出两层的模型参数,就需要...

    EM 算法所面对的问题跟之前的不一样,要复杂一些。

    EM 算法所用的概率模型,既含有观测变量,又含有隐变量。如果概率模型的变量都是观测变量,那么给定数据,可以直接用极大似然估计法,或贝叶斯估计法来估计模型参数,但是,当模型含有隐变量时,情况就复杂一些,相当于一个双层的概率模型,要估计出两层的模型参数,就需要换种方法求解。EM 算法是通过迭代的方法求解。

    监督学习是由训练数据 {(x(1),y(1)),(x(2),y(2)),...,(x(m),y(m))} 学习条件概率分布 P(Y|X) 或决策函数 Y=f(X) 作为模型,用于分类、回归等任务,这时训练数据中的每个样本点由输入和输出组成。但有时训练数据只有输入没有对应的输出 {(x(1),•),(x(2),•),...,(x(m),•)},从这样的数据学习模型称为非监督学习问题。EM 算法可用于生成模型的非监督学习,生成模型由联合概率分布 P(X,Y) 表示,可以认为非监督学习训练数据是联合概率分布产生的数据。X 为观测数据,Y 为未观测数据。

    我们先不管上一篇文章介绍的高斯混合模型,先来看通用的 EM 算法。

    假设有训练集 {x(1),x(2),...,x(m)},我们要寻找模型 p(x,z) 的参数来拟合这些数据,数据的似然估计为:

    要直接使用极大似然估计来求 θ 是很难的,因为 z(i) 是隐随机变量,如果 z(i) 不是隐变量,而是可以观察到的,那使用极大似然估计就简单多了。

    在这种情况下,EM 算法给出了一种高效的极大似然估计方法,直接最大化 L(θ) 很难,我们的策略是不断地构建 L 的下界(E-step),然后优化下界(M-step)。

    对于每个 i,使 Qi 为 z 上的分布(∑zQi(z)=1,Qi(z)≥0),那么

    最后一步的推导是用了 Jensen 不等式定理和 f(x)=log(x) 是凹函数的事实。

    现在,对于一些 Qi 的分布,上式给出了 L(θ) 的下界,Qi 有很多种选择,该选择哪个呢?如果已经有 θ 的猜测值,那么自然会想到让下界贴近 θ 值,也就是使不等式在 θ 处等号成立。根据 Jensen 不等式,如果要使 E[f(X)]=f(EX) 成立,应使 X 为常数。即

    也就是

    因为 ∑zQi(z(i))=1,所以可使

    所以,就简单地把 Qi 设为给定 x(i) 和 θ 后 z(i) 的后验分布即可。

    现在,有了 Q 的这个选择,就得到了 L 的下界,这是 E-step。在 M-step,对参数 θ 作极大似然估计。重复执行这两步,就是 EM 算法:

    我们怎么知道算法是不是收敛呢?假设 θ(t) 和 θ(t+1) 是 EM 算法连续迭代的两个参数, 我们现在证明 L(θ(t))≤L(θ(t+1)),以证明 EM 单调改进 log 似然。证明的关键就在于 Qi 的选择,不失一般性,我们从 EM 迭代的 θ(t) 开始,Qi(t)(z(i)) :=p(z(i))p(x(i)(t)),我们知道这使得 Jensen 不等式变恒等。

    参数 θ(t+1) 由最大化上式右边所得。

    第一个不等式来自一个事实:

    第二个不等式是因为 θ(t+1) 等于

    最后一个等式是因为 Qi 的选择使得 Jenson 不等式在 θ(t) 处等式成立。

    所以 EM 使似然单调收敛,EM 算法一直运行知道收敛,收敛测试就是看两次结果的差是不是小于一个设置的容忍值,如果 EM 改进很慢就说明收敛了。

    如果我们定义:

    从之前的推导我们知道 L(θ)≥J(Q,θ),EM 算法可以看作是 J 的坐标上升法,在 E-step,以 Q 为参数最大化,在 M-step,以 θ 为参数最大化。

     

    有了通用 EM 算法的定义,我们再来看下高斯混合模型中 Φ,µ 和 ∑ 的参数拟合。高斯混合模型应用广泛,在许多情况下,EM 算法是学习高斯混合模型的有效方法。简单起见,在 M-step 我们只推导 Φ 和 μi 的参数更新,有兴趣的可以推导下 ∑i

    假设有训练集 {x(1),x(2),...,x(m)},要把数据建模成一个联合分布 p(x(i),z(i))=p(x(i)|z(i))p(z(i)),其中 z(i)~Multinomial(Φ),x(i)|z(i)=j~N(µj,∑j),k 表示 z(i) 的可取值个数。数据 x(i) 是这样生成的:先从 {1,...,k} 中随机选择一个 z(i),再从 z(i) 所关联的高斯分布中生成 x(i)。这就是高斯混合模型,其中 z(i) 是隐变量,也就是未观测变量,正是这个变量使得问题变得复杂。

    这个模型的参数是 Φ,μ 和 ∑,我们的任务就是估计出这些参数。

    E-step 很简单,根据上面的推导,我们只需要计算

    这里 Qi(z(i)=j) 表示 z(i) 在分布 Qi 下取值 j 的概率。使用贝叶斯规则来计算 z(i) 的后验概率。

    这里 p(x(i)|z(i)=j;μ,∑) 是均值为 μj,方差为 ∑j 的高斯分布在 x(i) 处的密度,p(z(i)=j;Φ) 由 Φj 给出。

    下一步,在 M-step,我们需要最大化以 Φ,μ 和 ∑ 为参数的:

    先来推导以 μl 为参数的最大化,有

    设为 0,就得到 μl 的更新规则。

    跟上一篇文章结果一样。

    再来推导 M-step 中 Φj 的更新规则,经过观察,发现只需要最大化下式就可以了。

    还有一个附加的约束,ΣΦj=1,因为 Φj=p(z(i)=j;Φ),为了处理这个约束,我们构建拉格朗日式子:

    其中 β 是拉格朗日乘子。求导得:

    上式设为 0,解得:

    所以 Φj 正比于 Σw,由约束 ΣΦj=1,可得 -β=ΣiΣjw=m,所以

    对于 Σj 的推导也很直观。

     

    参考资料:

    [1] http://cs229.stanford.edu/notes/cs229-notes8.pdf

    [2] 李航,著.统计学习方法[M]. 清华大学出版社, 2012

    转载于:https://www.cnblogs.com/NaughtyBaby/p/5405576.html

    展开全文
  • 为使相位补偿误差影响分析结果并不依赖于哪一种特定的相位同步方案,给出了相参积累过程的一般数学模型和相位补偿误差的概率密度函数,推导了系统互模糊函数处理后的峰值输出,定义了相参积累损耗以分析相位同步误差的...
  • 然后,利用导出的TWOR-AFNC-Nodir系统中的紧闭形式下界解(theta = 2),我们研究了TWOR-AFNC-Dir系统以及两者之间的中断概率和遍历能力的整体比较系统模型。 在路径损耗模型基础上进行的比较分析表明,在城
  • 方法: 在原有部分变量的基础上,利用特征工程的方法,新增了单变量、双变量、三变量、四变量出现的频率和变量出现的条件概率等变量,利用随机森林模型,对目标变量进行预测。 结论: 一、对于训练集数据分析发现,...

    Amazon Employee Access 数据分析报告

    报告摘要

    • 目标:本分析旨在利用Amazon的员工编号相关信息,来分析和预测当员工申请访问某个编号的资源时,是否被允许访问。
    • 方法: 在原有部分变量的基础上,利用特征工程的方法,新增了单变量、双变量、三变量、四变量出现的频率和变量出现的条件概率等变量,利用随机森林模型,对目标变量进行预测。
    • 结论:
      • 一、对于训练集数据分析发现,各变量之间存在着一定的联系,其中ROLE_TITLE变量和ROLE_RODE变量存在一对一的关系,ROLE_TITLE变量和ROLE_FAMILY变量存在多对一的关系,其他变量之间也存在较强的对应关系。
      • 二、根据这种方法建模,发现模型具有一定的预测效果。

    目录

    • 问题描述
    • 数据加载
    • 数据探索
      • 描述性统计
      • 变量间的对应关系探索
      • 变量分布探索
    • 特征工程
      • 降维
      • 新增单变量频率
      • 新增双变量频率
      • 新增三变量频率
      • 新增四变量频率
      • 新增各变量出现的条件概率
    • 模型建立
    • 模型预测与评价

    一、问题描述

    利用Amazon员工的编号信息,包括员工经理的编号、员工所在分类的编号、员工所在部门编号、员工职位编号、员工类别编号等信息,来预测当员工申请访问某个 编号的资源时,是否被允许访问。

    变量名 含义
    ACTION 1代表资源被授权访问,0代表资源未被授权访问
    RESOURCE 资源编号
    MGR_ID 员工经理的编号
    ROLE_ROLLUP_1 公司员工分类1,如美国工程
    ROLE_ROLLUP_2 公司员工分类2,如美国零售
    ROLE_DEPTNAME 公司部门描述,如零售
    ROLE_TITLE 职位名称,如高级工程零售经理
    ROLE_FAMILY_DESC 公司员工类别扩展描述,如零售经理,软件工程
    ROLE_FAMILY 公司员工类别描述,如零售经理
    ROLE_CODE 员工角色编号

    二、数据加载

    加载所需的python库

    import statsmodels.api as sm
    import statsmodels.formula.api as smf
    import statsmodels.graphics.api as smg
    import patsy
    get_ipython().magic('matplotlib inline')
    import matplotlib.pyplot as plt
    import numpy as np
    import pandas as pd
    from pandas import Series,DataFrame
    from scipy import stats
    import seaborn as sns
    

    载入train数据集

    amazon = pd.read_csv("C:/Users/cs/Desktop/Amazon/train.csv")
    data =amazon
    data.head()
    
    ACTION RESOURCE MGR_ID ROLE_ROLLUP_1 ROLE_ROLLUP_2 ROLE_DEPTNAME ROLE_TITLE ROLE_FAMILY_DESC ROLE_FAMILY ROLE_CODE
    0 1 39353 85475 117961 118300 123472 117905 117906 290919 117908
    1 1 17183 1540 117961 118343 123125 118536 118536 308574 118539
    2 1 36724 14457 118219 118220 117884 117879 267952 19721 117880
    3 1 36135 5396 117961 118343 119993 118321 240983 290919 118322
    4 1 42680 5905 117929 117930 119569 119323 123932 19793 119325

    三、数据探索

    3.1 描述性统计

    train数据集共有32769个样本,不存在缺失值

    data.info()
    
    <class 'pandas.core.frame.DataFrame'>
    RangeIndex: 32769 entries, 0 to 32768
    Data columns (total 10 columns):
    ACTION              32769 non-null int64
    RESOURCE            32769 non-null int64
    MGR_ID              32769 non-null int64
    ROLE_ROLLUP_1       32769 non-null int64
    ROLE_ROLLUP_2       32769 non-null int64
    ROLE_DEPTNAME       32769 non-null int64
    ROLE_TITLE          32769 non-null int64
    ROLE_FAMILY_DESC    32769 non-null int64
    ROLE_FAMILY         32769 non-null int64
    ROLE_CODE           32769 non-null int64
    dtypes: int64(10)
    memory usage: 2.5 MB
    

    样本中,约有5.8%的员工授权申请没有通过,除ACTION外,各变量编号从个位数到六位数不等。

    data.describe()
    
    ACTION RESOURCE MGR_ID ROLE_ROLLUP_1 ROLE_ROLLUP_2 ROLE_DEPTNAME ROLE_TITLE ROLE_FAMILY_DESC ROLE_FAMILY ROLE_CODE
    count 32769.000000 32769.000000 32769.000000 32769.000000 32769.000000 32769.000000 32769.000000 32769.000000 32769.000000 32769.000000
    mean 0.942110 42923.916171 25988.957979 116952.627788 118301.823156 118912.779914 125916.152644 170178.369648 183703.408893 119789.430132
    std 0.233539 34173.892702 35928.031650 10875.563591 4551.588572 18961.322917 31036.465825 69509.462130 100488.407413 5784.275516
    min 0.000000 0.000000 25.000000 4292.000000 23779.000000 4674.000000 117879.000000 4673.000000 3130.000000 117880.000000
    25% 1.000000 20299.000000 4566.000000 117961.000000 118102.000000 118395.000000 118274.000000 117906.000000 118363.000000 118232.000000
    50% 1.000000 35376.000000 13545.000000 117961.000000 118300.000000 118921.000000 118568.000000 128696.000000 119006.000000 118570.000000
    75% 1.000000 74189.000000 42034.000000 117961.000000 118386.000000 120535.000000 120006.000000 235280.000000 290919.000000 119348.000000
    max 1.000000 312153.000000 311696.000000 311178.000000 286791.000000 286792.000000 311867.000000 311867.000000 308574.000000 270691.000000

    查看各变量上不同编号的种类数。可以发现,在30000多个样本中,RESOURCE、MGR_ID和ROLE_FAMILY上编号种类数较多,其他变量上编号种类数较少。
    值得注意的是,ROLE_TITLE和ROLE_CODE种类数一致。

    f = lambda x: x.unique().size
    data.apply(f)
    
    ACTION                 2
    RESOURCE            7518
    MGR_ID              4243
    ROLE_ROLLUP_1        128
    ROLE_ROLLUP_2        177
    ROLE_DEPTNAME        449
    ROLE_TITLE           343
    ROLE_FAMILY_DESC    2358
    ROLE_FAMILY           67
    ROLE_CODE            343
    dtype: int64
    

    3.2 变量间的对应关系探索

    3.2.1 ROLE_TITLE与ROLE_CODE

    画出ROLE_TITLE和ROLE_CODE变量的散点图,存在明显的正相关关系。

    fig,ax = plt.subplots(nrows=1,ncols=1,figsize=(8,5))
    plt.scatter(data.ROLE_TITLE,data.ROLE_CODE)
    
    <matplotlib.collections.PathCollection at 0xabb0c50>
    

    png

    将两个变量的值合并,编号的种类数目仍为343,

    TITLE_CODE = data.ROLE_TITLE*1000000+data.ROLE_CODE
    TITLE_CODE.unique().size
    
    343
    
    # 定义f2,用来计算交叉表每一行或每一列中非0值的个数
    f2 = lambda x: x[x!=0].count()
    # 画出两个变量间的交叉表
    TICO = pd.crosstab(data.ROLE_TITLE,data.ROLE_CODE)
    # 观察交叉表中ROLE_CODE变量对应的ROLE_TITLE变量个数
    TICO.apply(f2).plot()
    # 在变量ROLE_CODE上,对应的ROLE_TITLE个数为0,说明两个变量间至少存在一对多的对应关系
    TICO.apply(f2)[TICO.apply(f2)>1]
    
    Series([], dtype: int64)
    

    png

    
       观察交叉表中ROLE_TITLE变量对应的ROLE_CODE变量个数,也为0,说明两个变量间存在一一对应的关系
    TICO.apply(f2,axis=1).plot()
    TICO.apply(f2,axis=1)[TICO.apply(f2,axis=1)>1]
    
    Series([], dtype: int64)
    

    png

    3.2.2 ROLE_ROLLUP_1与ROLE_DEPTNAME

    # 将两个变量的值合并,编号的种类数目发生了较大的变化,但仍可发现,存在一定的对应关系
    RO1_DEP= data.ROLE_ROLLUP_1*10000000+data.ROLE_DEPTNAME
    data.ROLE_ROLLUP_1.unique().size, data.ROLE_DEPTNAME.unique().size, RO1_DEP.unique().size
    # ctRO1DEP = pd.crosstab(data.ROLE_ROLLUP_1,data.ROLE_DEPTNAME)
    
    (128, 449, 1185)
    

    3.2.3 ROLE_ROLLUP_2与ROLE_DEPTNAME

    # 将两个变量的值合并,编号的种类数目发生了较大的变化,但仍可发现,存在一定的对应关系
    RO2_DEP= data.ROLE_ROLLUP_2*10000000+data.ROLE_DEPTNAME
    data.ROLE_ROLLUP_2.unique().size, data.ROLE_DEPTNAME.unique().size, RO2_DEP.unique().size
    
    (177, 449, 1398)
    

    3.2.4 ROLE_ROLLUP_1与ROLE_ROLLUP_2

    # 将两个变量合并,编号的唯一值数目变化不大,说明两者之间存在很强的对应关系
    RO1_RO2= data.ROLE_ROLLUP_1*10000000+data.ROLE_ROLLUP_2
    data.ROLE_ROLLUP_1.unique().size, data.ROLE_ROLLUP_2.unique().size, RO1_RO2.unique().size
    
    (128, 177, 187)
    
    #画出两个变量间的交叉表
    ctRO12 = pd.crosstab(data.ROLE_ROLLUP_1,data.ROLE_ROLLUP_2)
    # 观察交叉表中ROLE_ROLLUP_2变量对应的ROLE_ROLLUP_1变量个数
    ctRO12.apply(f2).plot()
    # 在变量ROLE_ROLLUP_2上,只有三个值对应的ROLE_ROLLUP_1个数大于1(非一一对应关系),说明两个变量间有很强的一对多的对应关系
    ctRO12.apply(f2)[ctRO12.apply(f2)>1]
    
    ROLE_ROLLUP_2
    118164    2
    118178    2
    119256    9
    dtype: int64
    

    png

    # 统计ROLE_ROLLUP_2编号为118164、118178和119356样本的数目,样本数目的变量并不多,但总体上,未通过授权的比率比平均高
    a = data.ROLE_ROLLUP_2[(data.ROLE_ROLLUP_2==118164) | (data.ROLE_ROLLUP_2==118178)| (data.ROLE_ROLLUP_2==119256)].count()
    b = data.ACTION[(data.ROLE_ROLLUP_2==118164) | (data.ROLE_ROLLUP_2==118178)| (data.ROLE_ROLLUP_2==119256)].value_counts()
    b,a
    
    (1    380
     0     36
     Name: ACTION, dtype: int64, 416)
    
    # 观察交叉表中ROLE_ROLLUP_1变量对应的ROLE_ROLLUP_2变量个数
    # ctRO12.apply(f,axis=1).plot()
    # 在变量ROLE_ROLLUP_1上,有32个值对应的ROLE_ROLLUP_2个数大于1
    ctRO12.apply(f2,axis=1)[ctRO12.apply(f2,axis=1)>1].count()
    
    32
    

    3.2.5 ROLE_FAMILY与ROLE_FAMILY_DESC

    # 将两个变量合并,编号的唯一值数目变化不大,说明两者之间存在很强的对应关系
    FA_DESC= data.ROLE_FAMILY_DESC*1000000+data.ROLE_FAMILY   
    data.ROLE_FAMILY_DESC.unique().size,data.ROLE_FAMILY.unique().size, FA_DESC.unique().size
    
    (2358, 67, 2586)
    
    #画出两个变量间的交叉表
    ctFAFA = pd.crosstab(data.ROLE_FAMILY,data.ROLE_FAMILY_DESC)
    # 在变量ROLE_FAMILY_DESC上,有170个值对应的ROLE_FAMILY个数大于1,
    # 在变量ROLE_FAMILY上,有59个值对应的ROLE_FAMILY_DESC个数大于1,说明两个变量间有较强的一对多的对应关系
    

    3.2.6 ROLE_TITLE和ROLE_FAMILY

    # 将两个变量合并,唯一值没有发生变化,说明两者之间可能存在一对多关系
    TIFA = data.ROLE_TITLE*1000000+data.ROLE_FAMILY
    data.ROLE_TITLE.unique().size, data.ROLE_FAMILY.unique().size, TIFA.unique().size
    
    (343, 67, 343)
    
    #画出两个变量间的交叉表
    ctTIFA = pd.crosstab(data.ROLE_TITLE,data.ROLE_FAMILY)
    # 观察交叉表中ROLE_TITLE变量对应的ROLE_FAMILY变量个数
    ctTIFA.apply(f2,axis=1).plot()
    # 可以发现,ROLE_TITLE 与ROLE_FAMILY之间存在着一对多的关系,
    ctTIFA.apply(f2,axis=1)[ctTIFA.apply(f2,axis=1)>1].count()
    
    0
    

    png

    3.3 变量分布探索

    # 画出变量ACTION的条形图,大部分的申请都被授权
    fig,ax = plt.subplots(figsize=(8,5))
    data.ACTION.value_counts().plot(kind="bar",color="lightblue")
    ax.set_xticklabels(("Accessed","Not Accessed"),  rotation= "horizontal" )
    ax.set_title("Bar plot of Action")
    
    <matplotlib.text.Text at 0xd173080>
    

    png

    # 画出其余变量的分布直方图,RESOURCE和MGR_ID变量的编号大多分布在0-1000000上,且分布相对离散,其余变量分布都集中在一定的值和区域内。
    # 如变量ROLE_ROLLUP_1上,有21407个样本编号为117961;在ROLE_FAMILY上有10980个样本的编号为290919。
    # data.ROLE_ROLLUP_1.value_counts(),data.ROLE_FAMILY.value_counts()
    fig,ax = plt.subplots(nrows=4,ncols=2,figsize=(20,40))
    data.RESOURCE.hist(ax=ax[0,0],bins=100)
    ax[0,0].set_title("Hist plot of RESOURCE")
    data.MGR_ID.hist(ax=ax[0,1],bins=100)
    ax[0,1].set_title("Hist plot of MGR_ID")
    data.ROLE_ROLLUP_1.hist(ax=ax[1,0],bins=100)
    ax[1,0].set_title("Hist plot of ROLE_ROLLUP_1")
    data.ROLE_ROLLUP_2.hist(ax=ax[1,1],bins=100)
    ax[1,1].set_title("Hist plot of ROLE_ROLLUP_2")
    data.ROLE_DEPTNAME.hist(ax=ax[2,0],bins=100)
    ax[2,0].set_title("ROLE_DEPTNAME")
    data.ROLE_TITLE.hist(ax=ax[2,1],bins=100)
    ax[2,1].set_title("Hist plot of ROLE_TITLE")
    data.ROLE_FAMILY_DESC.hist(ax=ax[3,0],bins=100)
    ax[3,0].set_title("Hist plot of ROLE_FAMILY_DESC")
    data.ROLE_FAMILY.hist(ax=ax[3,1],bins=100)
    ax[3,1].set_title("Hist plot of ROLE_FAMILY")
    
    <matplotlib.text.Text at 0xdb2fa90>
    

    png

    # 画出变量间相关系数矩阵图,变量编号的值之间并没有明显的线性关系
    cm = np.corrcoef(data.values.T)
    sns.set(font_scale=1)
    cols = data.columns
    hm = sns.heatmap(cm, 
                cbar=True,
                annot=True, 
                square=True,
                fmt='.2f',
                annot_kws={'size': 10},
                yticklabels=cols,
                xticklabels=cols)
    plt.tight_layout()
    plt.show()
    

    png

    四、特征工程

    4.1 降维

    # 由于ROLE_CODE和ROLE_FAMILY与ROLE_TITLE存在一对一和一对多的关系,认为他不能包含更多的信息,删去这两个变量
    data = amazon
    del data["ROLE_CODE"]
    del data["ROLE_FAMILY"]
    
    amazon = pd.read_csv("C:/Users/cs/Desktop/Amazon/train.csv")
    

    4.2 新增单变量的频率

    # 利用循环,得到每个自变量出现的频率,赋值到新的列中。
    one= ["RESOURCE","MGR_ID","ROLE_ROLLUP_1","ROLE_ROLLUP_2","ROLE_DEPTNAME","ROLE_TITLE","ROLE_FAMILY_DESC"]
    for i in range(0,len(one)):
        a=data[one[i]]
        b=data[one[i]].value_counts()/32769
        a=a.map(b)
        data[one[i]+"_prob"]=a
    data.head()
    
    ACTION RESOURCE MGR_ID ROLE_ROLLUP_1 ROLE_ROLLUP_2 ROLE_DEPTNAME ROLE_TITLE ROLE_FAMILY_DESC RESOURCE_prob MGR_ID_prob ROLE_ROLLUP_1_prob ROLE_ROLLUP_2_prob ROLE_DEPTNAME_prob ROLE_TITLE_prob ROLE_FAMILY_DESC_prob
    0 1 39353 85475 117961 118300 123472 117905 117906 0.000092 0.001678 0.653270 0.135006 0.002197 0.109341 0.210443
    1 1 17183 1540 117961 118343 123125 118536 118536 0.000915 0.000305 0.653270 0.120388 0.004852 0.002472 0.000366
    2 1 36724 14457 118219 118220 117884 117879 267952 0.000061 0.000092 0.005615 0.005615 0.016662 0.038329 0.001007
    3 1 36135 5396 117961 118343 119993 118321 240983 0.000031 0.001892 0.653270 0.120388 0.005798 0.141872 0.037963
    4 1 42680 5905 117929 117930 119569 119323 123932 0.000244 0.000275 0.008423 0.004211 0.001373 0.002289 0.000580

    4.3 新增双变量的频率

    # 利用循环,得到每两个自变量同时出现的频率,赋值到新的列中。
    two = ["RESOURCE","MGR_ID","ROLE_ROLLUP_1","ROLE_ROLLUP_2","ROLE_DEPTNAME","ROLE_TITLE","ROLE_FAMILY_DESC"] 
    for i in range(0,len(two)):
        for j in range(i+1,len(two)):
            a=data[two[i]]+data[two[j]]*1000000
            b=a.value_counts()/32769
            a=a.map(b)
            data[two[i]+"_"+two[j]+"_prob"]=a
    data.head()
    
    ACTION RESOURCE MGR_ID ROLE_ROLLUP_1 ROLE_ROLLUP_2 ROLE_DEPTNAME ROLE_TITLE ROLE_FAMILY_DESC RESOURCE_prob MGR_ID_prob ... ROLE_ROLLUP_1_ROLE_ROLLUP_2_prob ROLE_ROLLUP_1_ROLE_DEPTNAME_prob ROLE_ROLLUP_1_ROLE_TITLE_prob ROLE_ROLLUP_1_ROLE_FAMILY_DESC_prob ROLE_ROLLUP_2_ROLE_DEPTNAME_prob ROLE_ROLLUP_2_ROLE_TITLE_prob ROLE_ROLLUP_2_ROLE_FAMILY_DESC_prob ROLE_DEPTNAME_ROLE_TITLE_prob ROLE_DEPTNAME_ROLE_FAMILY_DESC_prob ROLE_TITLE_ROLE_FAMILY_DESC_prob
    0 1 39353 85475 117961 118300 123472 117905 117906 0.000092 0.001678 ... 0.135006 0.002014 0.089200 0.180659 0.002014 0.013855 0.033233 0.000671 0.001678 0.079557
    1 1 17183 1540 117961 118343 123125 118536 118536 0.000915 0.000305 ... 0.120388 0.003815 0.002472 0.000366 0.003754 0.000580 0.000153 0.000153 0.000153 0.000366
    2 1 36724 14457 118219 118220 117884 117879 267952 0.000061 0.000092 ... 0.005615 0.000397 0.001556 0.000061 0.000397 0.001556 0.000061 0.005615 0.000061 0.000061
    3 1 36135 5396 117961 118343 119993 118321 240983 0.000031 0.001892 ... 0.120388 0.005401 0.125057 0.036956 0.005035 0.022460 0.007782 0.003052 0.001770 0.016204
    4 1 42680 5905 117929 117930 119569 119323 123932 0.000244 0.000275 ... 0.004211 0.000671 0.000488 0.000305 0.000549 0.000244 0.000244 0.000183 0.000183 0.000519

    5 rows × 36 columns

    4.4 新增三变量的频率

    # 利用循环,得到每三个自变量同时出现的频率,赋值到新的列中。
    three = ["RESOURCE","MGR_ID","ROLE_ROLLUP_1","ROLE_ROLLUP_2","ROLE_DEPTNAME","ROLE_TITLE","ROLE_FAMILY_DESC"] 
    for i in range(0,len(three)):
         for j in range(i+1,len(three)):
                for k in range(j+1,len(three)):
                    a = data[three[i]]*100000*100000+data[three[j]]*1000000+data[three[k]]
                    b = a.value_counts()/91690
                    a = a.map(b)
                    data[three[i]+"_"+three[j]+"_"+three[k]+"_"+"prob"]=a
    data.head()
    
    ACTION RESOURCE MGR_ID ROLE_ROLLUP_1 ROLE_ROLLUP_2 ROLE_DEPTNAME ROLE_TITLE ROLE_FAMILY_DESC RESOURCE_prob MGR_ID_prob ... ROLE_ROLLUP_1_ROLE_ROLLUP_2_ROLE_DEPTNAME_prob ROLE_ROLLUP_1_ROLE_ROLLUP_2_ROLE_TITLE_prob ROLE_ROLLUP_1_ROLE_ROLLUP_2_ROLE_FAMILY_DESC_prob ROLE_ROLLUP_1_ROLE_DEPTNAME_ROLE_TITLE_prob ROLE_ROLLUP_1_ROLE_DEPTNAME_ROLE_FAMILY_DESC_prob ROLE_ROLLUP_1_ROLE_TITLE_ROLE_FAMILY_DESC_prob ROLE_ROLLUP_2_ROLE_DEPTNAME_ROLE_TITLE_prob ROLE_ROLLUP_2_ROLE_DEPTNAME_ROLE_FAMILY_DESC_prob ROLE_ROLLUP_2_ROLE_TITLE_ROLE_FAMILY_DESC_prob ROLE_DEPTNAME_ROLE_TITLE_ROLE_FAMILY_DESC_prob
    0 1 39353 85475 117961 118300 123472 117905 117906 0.000092 0.001678 ... 0.000720 0.004951 0.011877 0.000185 0.000556 0.023220 0.000185 0.000556 0.003937 0.000218
    1 1 17183 1540 117961 118343 123125 118536 118536 0.000915 0.000305 ... 0.001341 0.000207 0.000055 0.000055 0.000055 0.000131 0.000055 0.000055 0.000055 0.000055
    2 1 36724 14457 118219 118220 117884 117879 267952 0.000061 0.000092 ... 0.000142 0.000556 0.000022 0.000055 0.000022 0.000022 0.000055 0.000022 0.000022 0.000022
    3 1 36135 5396 117961 118343 119993 118321 240983 0.000031 0.001892 ... 0.001800 0.008027 0.002781 0.001091 0.000633 0.005682 0.000971 0.000534 0.001451 0.000545
    4 1 42680 5905 117929 117930 119569 119323 123932 0.000244 0.000275 ... 0.000196 0.000087 0.000087 0.000044 0.000044 0.000109 0.000022 0.000022 0.000087 0.000065

    5 rows × 71 columns

    4.5 新增四变量的频率

    # 利用循环,得到每三个自变量和RESOURCE同时出现的频率,赋值到新的列中。
    four = ["RESOURCE","MGR_ID","ROLE_ROLLUP_1","ROLE_ROLLUP_2","ROLE_DEPTNAME","ROLE_TITLE","ROLE_FAMILY_DESC"] 
    for i in range(1,len(four)):
         for j in range(i+1,len(four)):
                for k in  range(j+1,len(four)):
                    a = data[four[0]]*100000*100000+data[four[i]]*1000000+data[four[j]]+data[four[k]]*0.000001
                    b = a.value_counts()/32769
                    a = a.map(b)
                    data[four[0]+"_"+four[i]+"_"+four[j]+"_"+four[k]+"_"+"prob"]=a
    data.head()
    
    ACTION RESOURCE MGR_ID ROLE_ROLLUP_1 ROLE_ROLLUP_2 ROLE_DEPTNAME ROLE_TITLE ROLE_FAMILY_DESC RESOURCE_prob MGR_ID_prob ... RESOURCE_ROLE_ROLLUP_1_ROLE_ROLLUP_2_ROLE_DEPTNAME_prob RESOURCE_ROLE_ROLLUP_1_ROLE_ROLLUP_2_ROLE_TITLE_prob RESOURCE_ROLE_ROLLUP_1_ROLE_ROLLUP_2_ROLE_FAMILY_DESC_prob RESOURCE_ROLE_ROLLUP_1_ROLE_DEPTNAME_ROLE_TITLE_prob RESOURCE_ROLE_ROLLUP_1_ROLE_DEPTNAME_ROLE_FAMILY_DESC_prob RESOURCE_ROLE_ROLLUP_1_ROLE_TITLE_ROLE_FAMILY_DESC_prob RESOURCE_ROLE_ROLLUP_2_ROLE_DEPTNAME_ROLE_TITLE_prob RESOURCE_ROLE_ROLLUP_2_ROLE_DEPTNAME_ROLE_FAMILY_DESC_prob RESOURCE_ROLE_ROLLUP_2_ROLE_TITLE_ROLE_FAMILY_DESC_prob RESOURCE_ROLE_DEPTNAME_ROLE_TITLE_ROLE_FAMILY_DESC_prob
    0 1 39353 85475 117961 118300 123472 117905 117906 0.000092 0.001678 ... 0.000092 0.000092 0.000092 0.000031 0.000031 0.000061 0.000031 0.000031 0.000061 0.000031
    1 1 17183 1540 117961 118343 123125 118536 118536 0.000915 0.000305 ... 0.000336 0.000305 0.000153 0.000031 0.000031 0.000031 0.000031 0.000031 0.000031 0.000031
    2 1 36724 14457 118219 118220 117884 117879 267952 0.000061 0.000092 ... 0.000061 0.000061 0.000031 0.000061 0.000031 0.000031 0.000061 0.000031 0.000031 0.000031
    3 1 36135 5396 117961 118343 119993 118321 240983 0.000031 0.001892 ... 0.000031 0.000031 0.000031 0.000031 0.000031 0.000031 0.000031 0.000031 0.000031 0.000031
    4 1 42680 5905 117929 117930 119569 119323 123932 0.000244 0.000275 ... 0.000061 0.000061 0.000061 0.000061 0.000061 0.000092 0.000031 0.000031 0.000061 0.000061

    5 rows × 91 columns

    4.6 新增各变量出现频率的条件概率

    # RESOURCE 确定时其他单个变量同时发生的概率
    resourcetwo = ['RESOURCE_MGR_ID_prob','RESOURCE_ROLE_ROLLUP_1_prob', 'RESOURCE_ROLE_ROLLUP_2_prob','RESOURCE_ROLE_DEPTNAME_prob', 
                   'RESOURCE_ROLE_TITLE_prob','RESOURCE_ROLE_FAMILY_DESC_prob']
    for i in range(0,len(resourcetwo)):
        a =  data[resourcetwo[i]]/data.RESOURCE_prob
        data[resourcetwo[i]+"_"+"probre"]=a
    data.head()
    
    ACTION RESOURCE MGR_ID ROLE_ROLLUP_1 ROLE_ROLLUP_2 ROLE_DEPTNAME ROLE_TITLE ROLE_FAMILY_DESC RESOURCE_prob MGR_ID_prob ... RESOURCE_ROLE_ROLLUP_2_ROLE_DEPTNAME_ROLE_TITLE_prob RESOURCE_ROLE_ROLLUP_2_ROLE_DEPTNAME_ROLE_FAMILY_DESC_prob RESOURCE_ROLE_ROLLUP_2_ROLE_TITLE_ROLE_FAMILY_DESC_prob RESOURCE_ROLE_DEPTNAME_ROLE_TITLE_ROLE_FAMILY_DESC_prob RESOURCE_MGR_ID_prob_probre RESOURCE_ROLE_ROLLUP_1_prob_probre RESOURCE_ROLE_ROLLUP_2_prob_probre RESOURCE_ROLE_DEPTNAME_prob_probre RESOURCE_ROLE_TITLE_prob_probre RESOURCE_ROLE_FAMILY_DESC_prob_probre
    0 1 39353 85475 117961 118300 123472 117905 117906 0.000092 0.001678 ... 0.000031 0.000031 0.000061 0.000031 1.000000 1.000000 1.000000 0.333333 0.666667 1.000000
    1 1 17183 1540 117961 118343 123125 118536 118536 0.000915 0.000305 ... 0.000031 0.000031 0.000031 0.000031 0.033333 0.866667 0.366667 0.033333 0.033333 0.033333
    2 1 36724 14457 118219 118220 117884 117879 267952 0.000061 0.000092 ... 0.000061 0.000031 0.000031 0.000031 0.500000 1.000000 1.000000 1.000000 1.000000 0.500000
    3 1 36135 5396 117961 118343 119993 118321 240983 0.000031 0.001892 ... 0.000031 0.000031 0.000031 0.000031 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000
    4 1 42680 5905 117929 117930 119569 119323 123932 0.000244 0.000275 ... 0.000031 0.000031 0.000061 0.000061 0.250000 0.375000 0.250000 0.250000 0.500000 0.375000

    5 rows × 97 columns

    # 其他单个变量确定时RESOURCE变量同时发生的概率
    resourcetwo = ['RESOURCE_MGR_ID_prob','RESOURCE_ROLE_ROLLUP_1_prob', 'RESOURCE_ROLE_ROLLUP_2_prob','RESOURCE_ROLE_DEPTNAME_prob', 
                   'RESOURCE_ROLE_TITLE_prob','RESOURCE_ROLE_FAMILY_DESC_prob']
    resourceone = [ 'MGR_ID_prob', 'ROLE_ROLLUP_1_prob','ROLE_ROLLUP_2_prob', 'ROLE_DEPTNAME_prob', 'ROLE_TITLE_prob','ROLE_FAMILY_DESC_prob']
    for i in range(0,len(resourcetwo)):
        a =  data[resourcetwo[i]]/data[resourceone[i]]
        data[resourcetwo[i]+"_"+"proboth"]=a
    data.head()
    
    ACTION RESOURCE MGR_ID ROLE_ROLLUP_1 ROLE_ROLLUP_2 ROLE_DEPTNAME ROLE_TITLE ROLE_FAMILY_DESC RESOURCE_prob MGR_ID_prob ... RESOURCE_ROLE_ROLLUP_2_prob_probre RESOURCE_ROLE_DEPTNAME_prob_probre RESOURCE_ROLE_TITLE_prob_probre RESOURCE_ROLE_FAMILY_DESC_prob_probre RESOURCE_MGR_ID_prob_proboth RESOURCE_ROLE_ROLLUP_1_prob_proboth RESOURCE_ROLE_ROLLUP_2_prob_proboth RESOURCE_ROLE_DEPTNAME_prob_proboth RESOURCE_ROLE_TITLE_prob_proboth RESOURCE_ROLE_FAMILY_DESC_prob_proboth
    0 1 39353 85475 117961 118300 123472 117905 117906 0.000092 0.001678 ... 1.000000 0.333333 0.666667 1.000000 0.054545 0.000140 0.000678 0.013889 0.000558 0.000435
    1 1 17183 1540 117961 118343 123125 118536 118536 0.000915 0.000305 ... 0.366667 0.033333 0.033333 0.033333 0.100000 0.001215 0.002788 0.006289 0.012346 0.083333
    2 1 36724 14457 118219 118220 117884 117879 267952 0.000061 0.000092 ... 1.000000 1.000000 1.000000 0.500000 0.333333 0.010870 0.010870 0.003663 0.001592 0.030303
    3 1 36135 5396 117961 118343 119993 118321 240983 0.000031 0.001892 ... 1.000000 1.000000 1.000000 1.000000 0.016129 0.000047 0.000253 0.005263 0.000215 0.000804
    4 1 42680 5905 117929 117930 119569 119323 123932 0.000244 0.000275 ... 0.250000 0.250000 0.500000 0.375000 0.222222 0.010870 0.014493 0.044444 0.053333 0.157895

    5 rows × 103 columns

    # RESOURCE 确定时其他两个变量同时发生的概率
    resourcethree = [ 'RESOURCE_MGR_ID_ROLE_ROLLUP_1_prob','RESOURCE_MGR_ID_ROLE_ROLLUP_2_prob', 'RESOURCE_MGR_ID_ROLE_DEPTNAME_prob',
            'RESOURCE_MGR_ID_ROLE_TITLE_prob','RESOURCE_MGR_ID_ROLE_FAMILY_DESC_prob', 'RESOURCE_ROLE_ROLLUP_1_ROLE_ROLLUP_2_prob',
            'RESOURCE_ROLE_ROLLUP_1_ROLE_DEPTNAME_prob','RESOURCE_ROLE_ROLLUP_1_ROLE_TITLE_prob','RESOURCE_ROLE_ROLLUP_1_ROLE_FAMILY_DESC_prob',
            'RESOURCE_ROLE_ROLLUP_2_ROLE_DEPTNAME_prob','RESOURCE_ROLE_ROLLUP_2_ROLE_TITLE_prob', 'RESOURCE_ROLE_ROLLUP_2_ROLE_FAMILY_DESC_prob',
            'RESOURCE_ROLE_DEPTNAME_ROLE_TITLE_prob', 'RESOURCE_ROLE_DEPTNAME_ROLE_FAMILY_DESC_prob','RESOURCE_ROLE_TITLE_ROLE_FAMILY_DESC_prob']
    for i in range(0,len(resourcethree)):
        a =  data[resourcethree[i]]/data.RESOURCE_prob
        data[resourcethree[i]+"_"+"probre"]=a
    data.head()
    
    ACTION RESOURCE MGR_ID ROLE_ROLLUP_1 ROLE_ROLLUP_2 ROLE_DEPTNAME ROLE_TITLE ROLE_FAMILY_DESC RESOURCE_prob MGR_ID_prob ... RESOURCE_ROLE_ROLLUP_1_ROLE_ROLLUP_2_prob_probre RESOURCE_ROLE_ROLLUP_1_ROLE_DEPTNAME_prob_probre RESOURCE_ROLE_ROLLUP_1_ROLE_TITLE_prob_probre RESOURCE_ROLE_ROLLUP_1_ROLE_FAMILY_DESC_prob_probre RESOURCE_ROLE_ROLLUP_2_ROLE_DEPTNAME_prob_probre RESOURCE_ROLE_ROLLUP_2_ROLE_TITLE_prob_probre RESOURCE_ROLE_ROLLUP_2_ROLE_FAMILY_DESC_prob_probre RESOURCE_ROLE_DEPTNAME_ROLE_TITLE_prob_probre RESOURCE_ROLE_DEPTNAME_ROLE_FAMILY_DESC_prob_probre RESOURCE_ROLE_TITLE_ROLE_FAMILY_DESC_prob_probre
    0 1 39353 85475 117961 118300 123472 117905 117906 0.000092 0.001678 ... 0.357389 0.119130 0.238259 0.357389 0.119130 0.238259 0.357389 0.119130 0.119130 0.238259
    1 1 17183 1540 117961 118343 123125 118536 118536 0.000915 0.000305 ... 0.131043 0.011913 0.011913 0.011913 0.011913 0.011913 0.011913 0.011913 0.011913 0.011913
    2 1 36724 14457 118219 118220 117884 117879 267952 0.000061 0.000092 ... 0.357389 0.357389 0.357389 0.178695 0.357389 0.357389 0.178695 0.357389 0.178695 0.178695
    3 1 36135 5396 117961 118343 119993 118321 240983 0.000031 0.001892 ... 0.357389 0.357389 0.357389 0.357389 0.357389 0.357389 0.357389 0.357389 0.357389 0.357389
    4 1 42680 5905 117929 117930 119569 119323 123932 0.000244 0.000275 ... 0.089347 0.089347 0.134021 0.134021 0.044674 0.089347 0.089347 0.089347 0.089347 0.134021

    5 rows × 118 columns

    # 其他两个变量确定时RESOURCE变量同时发生的概率
    resourcethree = [ 'RESOURCE_MGR_ID_ROLE_ROLLUP_1_prob','RESOURCE_MGR_ID_ROLE_ROLLUP_2_prob', 'RESOURCE_MGR_ID_ROLE_DEPTNAME_prob',
            'RESOURCE_MGR_ID_ROLE_TITLE_prob','RESOURCE_MGR_ID_ROLE_FAMILY_DESC_prob', 'RESOURCE_ROLE_ROLLUP_1_ROLE_ROLLUP_2_prob',
            'RESOURCE_ROLE_ROLLUP_1_ROLE_DEPTNAME_prob','RESOURCE_ROLE_ROLLUP_1_ROLE_TITLE_prob','RESOURCE_ROLE_ROLLUP_1_ROLE_FAMILY_DESC_prob',
            'RESOURCE_ROLE_ROLLUP_2_ROLE_DEPTNAME_prob','RESOURCE_ROLE_ROLLUP_2_ROLE_TITLE_prob', 'RESOURCE_ROLE_ROLLUP_2_ROLE_FAMILY_DESC_prob',
            'RESOURCE_ROLE_DEPTNAME_ROLE_TITLE_prob', 'RESOURCE_ROLE_DEPTNAME_ROLE_FAMILY_DESC_prob','RESOURCE_ROLE_TITLE_ROLE_FAMILY_DESC_prob']
    othertwo = ['MGR_ID_ROLE_ROLLUP_1_prob','MGR_ID_ROLE_ROLLUP_2_prob','MGR_ID_ROLE_DEPTNAME_prob', 'MGR_ID_ROLE_TITLE_prob',
           'MGR_ID_ROLE_FAMILY_DESC_prob', 'ROLE_ROLLUP_1_ROLE_ROLLUP_2_prob', 'ROLE_ROLLUP_1_ROLE_DEPTNAME_prob', 'ROLE_ROLLUP_1_ROLE_TITLE_prob',
           'ROLE_ROLLUP_1_ROLE_FAMILY_DESC_prob', 'ROLE_ROLLUP_2_ROLE_DEPTNAME_prob', 'ROLE_ROLLUP_2_ROLE_TITLE_prob',
           'ROLE_ROLLUP_2_ROLE_FAMILY_DESC_prob', 'ROLE_DEPTNAME_ROLE_TITLE_prob','ROLE_DEPTNAME_ROLE_FAMILY_DESC_prob',
           'ROLE_TITLE_ROLE_FAMILY_DESC_prob']
    for i in range(0,len(resourcethree)):
        a =  data[resourcethree[i]]/data[othertwo[i]]
        data[othertwo[i]+"_"+"proboth"]=a
    data.head()
    
    ACTION RESOURCE MGR_ID ROLE_ROLLUP_1 ROLE_ROLLUP_2 ROLE_DEPTNAME ROLE_TITLE ROLE_FAMILY_DESC RESOURCE_prob MGR_ID_prob ... ROLE_ROLLUP_1_ROLE_ROLLUP_2_prob_proboth ROLE_ROLLUP_1_ROLE_DEPTNAME_prob_proboth ROLE_ROLLUP_1_ROLE_TITLE_prob_proboth ROLE_ROLLUP_1_ROLE_FAMILY_DESC_prob_proboth ROLE_ROLLUP_2_ROLE_DEPTNAME_prob_proboth ROLE_ROLLUP_2_ROLE_TITLE_prob_proboth ROLE_ROLLUP_2_ROLE_FAMILY_DESC_prob_proboth ROLE_DEPTNAME_ROLE_TITLE_prob_proboth ROLE_DEPTNAME_ROLE_FAMILY_DESC_prob_proboth ROLE_TITLE_ROLE_FAMILY_DESC_prob_proboth
    0 1 39353 85475 117961 118300 123472 117905 117906 0.000092 0.001678 ... 0.000242 0.005415 0.000245 0.000181 0.005415 0.001574 0.000985 0.016245 0.006498 0.000274
    1 1 17183 1540 117961 118343 123125 118536 118536 0.000915 0.000305 ... 0.000997 0.002859 0.004412 0.029782 0.002906 0.018810 0.071478 0.071478 0.071478 0.029782
    2 1 36724 14457 118219 118220 117884 117879 267952 0.000061 0.000092 ... 0.003885 0.054983 0.014015 0.178695 0.054983 0.014015 0.178695 0.003885 0.178695 0.178695
    3 1 36135 5396 117961 118343 119993 118321 240983 0.000031 0.001892 ... 0.000091 0.002019 0.000087 0.000295 0.002166 0.000486 0.001402 0.003574 0.006162 0.000673
    4 1 42680 5905 117929 117930 119569 119323 123932 0.000244 0.000275 ... 0.005180 0.032490 0.067010 0.107217 0.019855 0.089347 0.089347 0.119130 0.119130 0.063069

    5 rows × 133 columns

    五、模型建立

    # 划分测试集与训练集
    from sklearn.cross_validation import train_test_split
    from sklearn.ensemble import RandomForestClassifier
    from sklearn.metrics import confusion_matrix, roc_curve,roc_auc_score,classification_report 
    y = data.ACTION
    X = data
    del X["ACTION"]
    X_train, X_test, y_train, y_test = train_test_split(
             X, y, test_size=0.3, random_state=0)
    
    # 利用以上处理所得的共133个自变量建立随机森林模型
    forest = RandomForestClassifier(criterion='entropy',
                                    n_estimators=1000, 
                                    random_state=1,
                                    n_jobs=2)
    RFfit = forest.fit(X_train , y_train)
    

    六、模型预测与评价

    # 利用模型进行预测
    preds = RFfit.predict(X_test)
    
    # 得到模型的混淆矩阵如下所示
    confusion_matrix(y_test,preds)
    
    array([[ 138,  420],
           [  59, 9214]])
    
    # 得到模型的ROC_AUC得分如下所示
    pre = RFfit.predict_proba(X_test)
    roc_auc_score(y_test,pre[:,1])
    
    0.8639483844684166
    
    # 得到摸型的ROC曲线如下所示
    fpr,tpr,thresholds = roc_curve(y_test,pre[:,1])
    fig,ax = plt.subplots(figsize=(8,5))
    plt.plot(fpr,tpr)
    ax.set_title("Roc of Logistic Randomforest")
    
    <matplotlib.text.Text at 0x26395198>
    

    png

    利用Kaggle测试集得分为0.89,说明模型具有一定的效果。

    展开全文
  • 转至:http://www.hankcs.com/ml/em-algorithm-and-its-generalization.html本文是...如果概率模型变量都是观测变量(数据中可见的变量),则可以直接用极大似然估计,或者用贝叶斯估计模型参数。但是,当模型含有...
  • EM算法及其推广

    2017-12-11 15:44:12
    本文是《统计学习方法》第九章的笔记,注解了原著的部分公式推导,补充了另一个经典的硬币模型,并且...如果概率模型变量都是观测变量(数据中可见的变量),则可以直接用极大似然估计,或者用贝叶斯估计模型
  • 文章目录小结在前由概率到分数的转换概率转换分数最终表达式变量的分值计算评分卡性能评估坏账率与通过率的审批策略双卡审批策略 将模型预测概率转化为分数,更符合人的一个直观感受 小结在前 1.把对数几率当成分数...
  • 概率 最大似然估计——用频率代替概率 ...衡量模型和真实概率之间的差距 困惑度 困惑度和交叉熵呈指数关系 互信息 指X中包含的Y的信息量 即知道了X之后,能知道多少Y 字耦合度 噪声信道模型 ...
  • 神经网络

    2019-12-20 21:03:27
    输出观测概率 初始状态概率 参数学习 观测变量一样,模型不一样,能解英文的码 解码:根据观测序列推测隐藏的模型状态 基本算法: 所有可能性: 图的性质:只跟前一状态有关系。 前项算似然 马尔科夫随机场 分部...
  • 简要介绍了印度喀拉拉邦数学... 表征-密度表征,信息度量,公理定义,矩阵自变量的伪解析函数以及正态概率定律的表征; Mathai的熵-熵优化; 方差分析; 人口问题与社会科学; 二次和双线性形式; 线性代数 概率统计。
  • RBM受限玻尔兹曼机是组成深度信念网络的一种基本单元,RBM是一种概率生成模型也是一种能量模型, RBM模块中主要包括一个隐含层和一个可见层。其中,受限玻尔兹曼机的隐含层和可见层是通过双向连接的方式进行连接的,...
  • 语音识别--gmm-hmm思考

    2020-03-01 01:29:51
    简单回顾一下今天所看的内容: gmm-hmm pdf: 概率密度函数,在这里可以由gmm来估计,同样也可以用...hmm:隐马尔科夫模型马尔科夫链的过程。关键在于理解状态。 首先需要的说的马尔科夫链。当与时间无关时,...
  • 5.4.1 顺序概率模型 207 5.4.2 分解模型 210 5.5 持续期模型 214 5.5.1 ACD模型 216 5.5.2 模拟 218 5.5.3 估计 219 5.6 非线性持续期模型 224 5.7 价格变化和持续期的二元模型 225 5.8 应用 229 附录A ...
  • 该方法首先建立字典原子及各参数的概率分布模型, 将其划分为局部变量及全局变量, 并使用固定其他参数来更新当前参数的Gibbs抽样方法对各变量赋予初始值, 然后采用随机优化方法对两类变量进行期望最大化(EM)迭代优化,...
  • 《数学之美》中的自然语言处理

    千次阅读 2017-02-09 12:30:09
    1. 信息的冗余是信息安全的保障。 2. 语言的数据,我们称之为语料,尤其是双语或者多语的对照语料对翻译至关重要...当使用条件概率的时候,每个变量的可能性就是一种语言字典的大小,而如果假定当前的词出现的概率
  • 方法:基于总共10,948个碰撞记录,建立了一个二进制logit模型,以探索各种因素对行人死亡概率的影响。 此外,一级相互作用效应被引入基本模型。 Hosmer-Lemeshow拟合优度检验用于评估模型性能。 计算分类变量的赔率...
  • 什么是统计套利

    千次阅读 2017-07-28 11:22:44
    统计套利主要是在对历史数据进行统计分析的基础上,估计相关变量概率分布,并结合基本面数据进行分析以指导套利交易,其主要分为配对/一揽子交易、多因素模型、均值回归策略、协整以及针对波动率与相关性的建模,...

空空如也

空空如也

1 2 3
收藏数 53
精华内容 21
关键字:

双变量概率模型