15. R的bnlearn包中的per.node.score函数的学习探究

最新推荐文章于 2023-05-18 16:18:45 发布

Rilkean heart

最新推荐文章于 2023-05-18 16:18:45 发布

阅读量721

点赞数

分类专栏： bnlearn R语言

本文链接：https://blog.csdn.net/qq_42147816/article/details/105000794

版权

R语言同时被 2 个专栏收录

18 篇文章 2 订阅

订阅专栏

bnlearn

15 篇文章 4 订阅

订阅专栏

0. 写在前面

原来只知道per.node.score是对贝叶斯网络单个节点进行打分，具体实现从来没有触及过，最近终于努力摁着头学习了一下。诸多不足，欢迎交流。
此文比较适合有一定使用bnlearn包基础以及C语言基础并且也为这个函数的机理头疼的人观看
本文使用的数据是一种n*n的基因表型数据，暂把它理解为列名和行名都表示基因，不理解数据不影响阅读。
csdn竟然不支持R代码的高亮，故用python的高亮代替，不是很美貌QAQ
这大概是我写过的最长的blog了，研究了快一周，写博客用了将近10h QAQ

1.per.node.score()与score.R

这一切都要从R看起，下面是我的R代码per.node.score的一段调用，部分参数已说明。
当然这里想计算“单节点分数”有两种方法，一是用score函数外加参数by.node=TRUE

node.score = score(start,x,type = "bic-g",by.node=TRUE)

第二种方法，也就是上面这个方法的本质，就是直接用per.node.score()。这里还是推荐第一种。

# network:bn类型的网络；score：打分类型，这里是bic-g也就是bic的连续数据打分类型；
# targets：传进来的数据的列名，在这里是各种基因；data: 原数据集，dataframe类型(应该是)
node.score = per.node.score(network = start, score = score,
            targets = nodes, extra.args = extra.args, data = x)

这个函数扒开它第一层外皮，来到score.R文件中，发现十hao分wu严yong谨chu，只是告诉你这里调用了per_node_score这个C函数。

# compute individual node contributions to the network score.
per.node.score = function(network, data, score, targets, extra.args,
    debug = FALSE) {

  .Call("per_node_score",
        network = network,
        data = data,
        score = score,
        targets = targets,
        extra.args = extra.args,
        debug = debug)

}#PER.NODE.SCORE

1.5 per_node_score()与score.h

这里第1.5层皮是从这个头文件中可以找到如下代码，SEXP是指针，属于高级用法，具体入门可以看看这个链接：

http://adv-r.had.co.nz/

（这种“高级”的用法就是促使我那个学姐把这一大串函数一怒之下改成python的原因。）

在这里per_node_score的具体实现要看向c_per_node_score这个函数

/* from per.node.score.c */
SEXP per_node_score(SEXP network, SEXP data, SEXP score, SEXP targets,
    SEXP extra_args, SEXP debug);
void c_per_node_score(SEXP network, SEXP data, SEXP score, SEXP targets,
    SEXP extra_args, int debuglevel, double *res);

2. c_per_node_score与per.node.score.c

第二层：c_per_node_score()在per.node.score.c中。
这个函数代码乍一看很长，但是都是switch分支，只要看到我们这里需要的那个case就好，这里的case是根据传进来的参数score打分类型进行区分的。此样例中是bic-g也就是bic的连续数据打分类型；
从下面的代码中可以看到我们下一步是看看loglik_gnode的实现。
有部分代码解释以注释方式给出： // wyj：XXXX

 /* AIC and BIC scores, Gaussian data. */
 // wyj： 可见aic和bic在这里是一样的实现逻辑
    case AIC_G:
    case BIC_G:

      k = REAL(getListElement(extra_args, "k"));

      for (i = 0; i < ntargets; i++) {

        SET_STRING_ELT(cur, 0, STRING_ELT(targets, i));
        DEBUG_BEFORE();
		
		//wyj: 上面都看不懂没关系，这一步调用是核心操作
        res[i] = loglik_gnode(cur, network, data, &nparams, debuglevel);
        //wyj: 可以看到res存放了返回值，，下面是惩罚值的减去
        res[i] -= (*k) * nparams;

        if (debuglevel > 0)
          Rprintf("  > penalty is %lf x %.0lf = %lf.\n", *k, nparams, (*k) * nparams);

      }/*FOR*/
      break;

3. loglik_gnode()与gaussian.loglikelihood.c

第三层皮： loglik_gnode()在gaussian.loglikelihood.c中。
此函数是在计算（连续）对数似然估计值，我浅薄的高数基础只会求个导看是否是极值点but在这里毫无用处。因此可以看一看这些博客补课：

https://blog.csdn.net/wydbyxr/article/details/83212703

依旧定位到我们需要的函数实现部分

double loglik_gnode(SEXP target, SEXP x, SEXP data, double *nparams, int debuglevel) {

double loglik = 0;
char *t = (char *)CHAR(STRING_ELT(target, 0));
SEXP nodes, node_t, parents, data_t;

  /* get the node cached information. */
  nodes = getListElement(x, "nodes");
  node_t = getListElement(nodes, t);
  /* get the parents of the node. */
  parents = getListElement(node_t, "parents");
  /* extract the node's column from the data frame. */
  data_t = c_dataframe_column(data, target, TRUE, FALSE);

  /* compute the log-likelihood. */
  // wyj: 这里才是我们重点要看的
  if (length(parents) == 0)
    loglik = glik(data_t, nparams);
  else
    loglik = cglik(data_t, data, parents, nparams);

  if (debuglevel > 0)
    Rprintf("  > loglikelihood is %lf.\n", loglik);

  return loglik;

}/*LOGLIK_GNODE*/

这么一大段也就这两句是我们需要的吧：
这里可以看到，计算节点分数分为了两种情况，一是无父节点，一是有父节点。
这里先看parents == 0的情况

/* compute the log-likelihood. */
  if (length(parents) == 0)
    loglik = glik(data_t, nparams);
  else
    loglik = cglik(data_t, data, parents, nparams);

3.1.无父节点：glik()

第四层皮：该节点无父节点的情况下，glik()仍然在gaussian.loglikelihood.c中：
代码结构是：先算了这一列基因数据的平均值xm，又算了标准差sd。sd如果<0认为res返回值，也就是最后的“打分分数”为负无穷；否则计算dnorm。

double glik(SEXP x, double *nparams) {

int i = 0, num = length(x);
double *xx = REAL(x), xm = 0, res = 0;
long double sd = 0;

  /* compute the mean, */
  // wyj:平均值
  for (i = 0; i < num; i++)
    xm += xx[i];
  xm /= num;

  /* compute the standard deviation. */
  // wyj:标准差
  SD_GUARD(num, 1, sd,
    for (i = 0; i < num; i++)
      sd += (xx[i] - xm) * (xx[i] - xm);
    sd = sqrt(sd / (num - 1));
  )

  /* compute the log-likelihood (singular models haze zero density). */
  if (sd < MACHINE_TOL)  //wyj: 我这里理解(浮点数)sd小于0
    res = R_NegInf;      /*wyj：负无穷*/
  else
    for (i = 0; i < num; i++)
      res += dnorm(xx[i], xm, (double)sd, TRUE);

  /* we may want to store the number of parameters. */
  if (nparams)
    *nparams = 1;

  return res;

}/*GLIK*/

这里dnorm我没有找到C实现，我认为是调用的R函数dnorm，虽然知道C可以再去调用R但是对于具体的实现毫无头绪（希望有大神可以甩链接甩教程）。不过这里用python30行代码就可以轻松实现上面这“四层皮”，对比了结果验证了一样之后，认为这里是（分别）计算了我这个数据集的这一列中每一个基因的标准正态的概率密度值。贴上这一句改写的python代码：

       for i in range(0, num):   # 依次累积该基因对其他基因的res值 
            # norm是正态化函数，样本总体的均值和标准差是两个需要的参数，
            # pdf：概率密度函数，对于连续函数，pdf是变量值为x的概率。
            # 这里log的返回值是标准正态概率密度再取log值，最大似然函数的公式就是这样写的
            res = res + log(scipy.stats.norm(xm, sd).pdf(node[i])) 
        res = res - nparams       # 最后的累积值-nparams即为单节点打分分数

3.2 有父节点：cglik

第四层皮：该节点有父节点的情况下，计算稍微复杂了一些，而且又多了几层调用而且传入的参数的名字会变来变去，已在代码注释中做了说明。
有父节点情况下使用了cglik()函数，该函数仍然在gaussian.loglikelihood.c中：
首先是做了一些数据的处理与空间的分配：

/*x：target目标打分节点那一列的数据 data：数据集*/
double cglik(SEXP x, SEXP data, SEXP parents, double *nparams) {
int i = 0, nrow = length(x);
int ncol = length(parents); // wyj: 注意ncol是parents个数，即父节点个数
double *xx = REAL(x), **dd = NULL, res = 0, *fitted = NULL, sd = 0;
SEXP data_x;

  /* dereference the data and prepare it for the QR decomposition. */
  // wyj：取消引用数据并为QR分解做好准备
  // wyj：data_x也就是下面的dd，存放了parents列的data值
  PROTECT(data_x = c_dataframe_column(data, parents, FALSE, FALSE));
  dd = Calloc1D(ncol, sizeof(double *));
  for (i = 0; i < ncol; i++)
    dd[i] = REAL(VECTOR_ELT(data_x, i));
  /* allocate the fitted values. */
  fitted = Calloc1D(nrow, sizeof(double));

然后比较关键的一步： fitted和sd将在这里被计算

 // wyj: 尽可能利用特殊情况(根据ncol分类，即父节点个数)有效地计算最小二乘
  // wyj: fitted和sd在这里被计算
  /*dd：parents列的数据值 
    xx：target目标打分节点那一列的数据 ncol：parents数，即父节点个数 */
  c_ols(dd, xx, nrow, ncol, fitted, NULL, NULL, &sd);

然后这一部分和无父节点的情况基本一样，都调用了dnorm函数去计算res，只是dnorm传进去的参数不一样：

//wyj :计算对数似然(奇异模型密度为零)
  /* compute the log-likelihood (singular models have zero density). */
  if (sd < MACHINE_TOL)
    res = R_NegInf;
  else
    for (i = 0; i < nrow; i++)
    // wyj: 和没有父节点相比，这里把平均数xm换成了fitted
      res += dnorm(xx[i], fitted[i], (double)sd, TRUE);

  /* we may want to store the number of parameters. */
  if (nparams)
    *nparams = ncol + 1;

由此可以看出，这里关键是c_ols函数中的最小二乘法的实现过程。
关于最小二乘法可以看看这个博客以及这个知乎问题，只需要一点点线代的知识：

https://blog.csdn.net/fu6543210/article/details/84380508
https://www.zhihu.com/question/37031188

3.2.1 c_ols()与least.squares.c

第五层皮：c_ols()函数在least.squares.c中。注意这里传进去的参数的含义。
这个函数的功能是：尽可能利用特殊情况(根据ncol分类，即父节点个数)有效地计算最小二乘
首先是传参和数据处理，然后可以看到分为3种情况分别处理：没有父节点/1个父节点/多个父节点。注意ncol表示父节点个数。

/* compute least squares efficiently by special-casing whenever possible. */
/* x；parents列的数据值 y: target目标打分节点那一列的数据 ncol: 父节点个数
   fitted ：适应值，类型为nrow大小的double*数组
*/
void c_ols(double **x, double *y, int nrow, int ncol, double *fitted,
    double *resid, double *beta, double *sd) {
  double *qr = NULL;

	if (ncol == 0) {

    // wyj: 没有父节点
    /* special case: null model. */
    c_ols0(y, nrow, fitted, resid, beta, sd);

  }/*THEN*/
  
  else if (ncol == 1) {

    // wyj: 父节点为1————简单回归
    /* special case: simple regression. */
    c_ols1(x, y, nrow, fitted, resid, beta, sd);

  }/*THEN*/
  else {
		//有些长，这里先不贴上代码
  }/*ELSE*/

3.2.1.1 没有父节点：c_ols0()

首先看看没有父节点的情况，此时调用c_ols0函数，也在least.squares.c中

// wyj: 快速实现最小二乘的截距
/* fast implementation of least squares with just the intercept. */
/*y: target目标打分节点那一列的数据 */
static void c_ols0(double *y, int nrow, double *fitted, double *resid,
    double *beta, double *sd) {
    
	int i = 0;
	double mean = c_mean(y, nrow);

   // wyj: 唯一的系数是截距，也就是均值
   /* the only coefficient is the intercept, which is the mean of the response. */
  if (beta)
    *beta = mean;
  if (sd)
    c_sd(y, nrow, 1, mean, FALSE, sd);
  if (fitted)
    for (i = 0; i < nrow; i++)
      fitted[i] = mean;
  if (resid)
    for (i  = 0; i < nrow; i++)
      resid[i] = y[i] - mean;
}/*C_OLS0*/

3.2.1.2 1个父节点：c_ols1()

其次是有1个父节点的情况，此时调用c_ols1函数，也在least.squares.c中。
这个函数的框架是：先算出回归系数a与b，再用a与b去计算其他的东西。
首先是一堆变量声明：

// wyj : 快速实现最小二乘与一个回归系数。
/* fast implementation of least squares with one regression coefficient. */
/*x；parents列的数据值 y: target目标打分节点那一列的数据*/
static void c_ols1(double **x, double *y, int nrow, double *fitted, double *resid,
    double *beta, double *sd) {

int i = 0, singular = FALSE;
double *m[2] = {y, *x}, mean[2] = {0, 0}, cov[4] = {0, 0, 0, 0}, a = 0, b = 0;
long double sse = 0;

然后利用闭式估计或者闭式解计算回归系数a与b：调用的这两个函数的具体实现在covariance.c中，比较简单不再深究。

/* the regression coefficients are computed using the closed form estimators
   * for simple regression. */
  /*利用闭式估计计算回归系数*/
  //wyj: 下标0和1分别表示target和parents
  c_meanvec(m, mean, nrow, 2, 0);      //wyj: 从一组变量中计算均值向量
  //wyj：cov是一个压缩矩阵（一维数组）
  c_covmat(m, mean, nrow, 2, cov, 0);  //wyj: 填写协方差矩阵

什么是闭式解？

闭式解closed form solution）也叫解析解(analytical solution)，就是一些严格的公式,给出任意的自变量就可以求出其因变量,也就是问题的解, 他人可以利用这些公式计算各自的问题。比如我们最常见的二元一次方程的解的公式。
https://www.cnblogs.com/peanutk/p/10207051.html

然后接下来判断奇异性。MACHINE_TOL我们上面提到过，可以看做机器0，实际上是一个比0大那么一丢丢丢丢的数。

/* the only way for a simple regression to be singular is if the regressor
   * is constant, so it is collinear with the intercept. */
  /*简单回归方程显示出奇异性的唯一情况是回归方程为常数，也就是它与应该与截距共线。*/
  //wyj: cov[3]为cov[parent,parent]也就是D(parent)方差
  //wyj: 当父节点方差基本上为0时，说明父节点那一列的值不发生什么变化，
  //wyj: 此时图像是一条平行x轴的直线，singular为1作为奇异性的标志
  singular = (fabs(cov[3]) < MACHINE_TOL);

  if (singular) {

    b = NA_REAL;
    a = mean[0];

  }/*THEN*/

什么是奇异性呢？

奇异型分布(singular distribution)一类理论上很有价值的概率分布函数，实际应用中也经常遇到。如果分布函数F (x)连续，但它的导函数几乎处处(关于勒贝格测度而言)等于0，则称分布函数F(x)为奇异型的.如果随机变量X的分布函数是奇异型的，则称该随机变量及其概率分布P为奇异型的.
对奇异型分布函数F（x），不存在密度函数f（x）。
https://baike.baidu.com/item/%E5%A5%87%E5%BC%82%E5%9E%8B%E5%88%86%E5%B8%83/19096761
这里是判断奇异性的代码

如果是非奇异性的话，b和a分别这样计算

else {

    b = cov[1] / cov[3];       //wyj: b=cov[target,parent]/D(parent)
    a = mean[0] - mean[1] * b; //wyj: a=mean(target)=mean(parent)*b

  }/*ELSE*/

实际上从这里可以看到，奇异性情况下b为NA，因为b的公式分母cov[3]也就是D(parent)为0，而分母是8可以为0的。这两种情况的公式是相通的。
关于回归系数为什么要这样算呢？数学依据在这里：

https://www.zhihu.com/question/58223009
https://wenku.baidu.com/view/cf5a4f8f84254b35eefd34b7.html

然后我们就可以继续用求出来的一元线性回归表达式去求其他的项。在这里我只用到了fitted预测值，和sd标准差。至于什么是resid这里不再探究。

 //wyj: 计算适应值，或者说预测值
  if (fitted) {
   //Rprintf("  > wyj-test2-3.26-fitted < \n");

    if (singular)
      for (i = 0; i < nrow; i++)
        fitted[i] = a;
    else
      for (i = 0; i < nrow; i++)
        fitted[i] = a + b * x[0][i];  //wyj: 一元线性回归

  }/*THEN*/

  //wyj: what is resid?
  if (resid) {
    //wyj: 下面这句并没有输出
    //Rprintf("  > wyj-test2-3.26-resid < \n");
    if (fitted) 
      for (i  = 0; i < nrow; i++)
        resid[i] = y[i] - fitted[i];
    else {

      if (singular)
        for (i  = 0; i < nrow; i++)
          resid[i] = y[i] - a + b * x[0][i];
      else
        for (i  = 0; i < nrow; i++)
          resid[i] = y[i] - a;

    }/*ELSE*/

  }/*THEN*/

  if (sd) {
     //Rprintf("  > wyj-test2-3.26-sd < \n");
    SD_GUARD(nrow, 2, *sd,
      if (resid)
        for (i = 0; i < nrow; i++)
          sse += resid[i] * resid[i];
      else if (fitted){                                    //wyj: 这里是根据fitted算出sd
        //Rprintf("  > wyj-test2-3.27-fitted < \n");
        for (i = 0; i < nrow; i++)
          sse += (y[i] - fitted[i]) * (y[i] - fitted[i]);  //wyj: sse：和方差；y: target
      } else {
        //wyj: 下面的语句并未输出
        //Rprintf("  > wyj-test2-3.27-sd-3 < \n");
        if (singular)
          sse = c_sse(y, a, nrow);
        else{                                               
          //wyj: 下面的语句并未输出
          //Rprintf("  > wyj-test2-3.27-4 < \n");
          for (i = 0; i < nrow; i++)
            sse += (y[i] - a - b * x[0][i]) * (y[i] - a - b * x[0][i]);
        }
      }/*ELSE*/
	  //wyj: sse/(n-p-1) 是误差项的方差的无偏估计。p表示自变量个数
      *sd = sqrt(sse / (nrow - 2));
    )

3.2.1.3 多个父节点：c_qr_matrix()和c_qr()

从c_ols中可以看到，有多个父节点需要qr分解，并且涉及到2个函数c_qr_matrix()和c_qr()

/* prepare the data for the QR decomposition. */
    qr = Calloc1D(nrow * (ncol + 1), sizeof(double));
    c_qr_matrix(qr, x, nrow, ncol);
    /* general case: multiple regression. */
    c_qr(qr, y, nrow, ncol + 1, fitted, resid, beta, sd);

关于qr分解与最小二乘法的结合可以看这些资料：

https://baike.baidu.com/item/QR%E5%88%86%E8%A7%A3/8918473?fr=aladdin
https://www.cnblogs.com/AndyJee/p/3846455.html

3.2.1.3.1 c_qr_matrix()与common.c

c_qr_matrix()在common.c中，这个函数很简单，目的就是对qr这个结构初始化为1，并且将父节点数据复制到qr矩阵的最右列。这种结构是qr分解需要的特殊结构：第一列为全1，后面几列为父节点的数据。

void c_qr_matrix(double *qr, double **x, int nrow, int ncol) {

int i = 0;

  //wyj: 填充截距列
  /* fill the intercept column. */
  for (i = 0; i < nrow; i++)
    qr[i] = 1;

  //wyj: 将父节点数据复制到qr矩阵的最右列。
  /* copy the data to the right comun(column) of the matrix. */
  for (i = 0; i < ncol; i++)
    memcpy(qr + (i + 1) * nrow, x[i], nrow * sizeof(double));

}/*C_QR_MATRIX*/

注意下memcpy的用法

void *memcpy(void *str1, const void *str2, size_t n)
参数：
str1 – 指向用于存储复制内容的目标数组，类型强制转换为 void*指针。
str2 – 指向要复制的数据源，类型强制转换为 void* 指针。
n – 要被复制的字节数。

3.2.1.3.2 c_qr()与linear.algebra.c

c_qr这个函数是一个C级函数接口，通过QR分解执行OLS（普通最小二乘法）。位于linear.algebra.c中。不再展开具体，调用就完事了。

/* general case: multiple regression. */
    //wyj: 多元回归
    c_qr(qr, y, nrow, ncol + 1, fitted, resid, beta, sd);

4. 总结（思维导图）

我把bnlearn中使用bic打分制成了一个思维导图：
在这里插入图片描述
查看更清晰的版本（xmind、png、pdf）可以点此下载

5. 暂留的疑问

1. bnlearn的单节点打分？

bnlearn代码中的src\gaussian.loglikelihood.c中有：

/* compute the log-likelihood. */
  if (length(parents) == 0)
    loglik = glik(data_t, nparams);
  else
    loglik = cglik(data_t, data, parents, nparams);

显然这里是把parent==0和parent！=0分开算了节点分数。后者调用了cglik()，但是cglik()中却有这样的一段

if (ncol == 0) {

    // wyj: ncol表示父节点数量，此时没有父节点
    /* special case: null model. */
    c_ols0(y, nrow, fitted, resid, beta, sd);

  }/*THEN*/

为什么在parent！=0中又有一段parent==0的逻辑？尽管我做了一下测试，这一个代码块并不会被调用，但是它的存在还是很迷惑。ncol和parent究竟是不是一个意思？为什么这里又会有这样一段逻辑？

2. 惩罚值？

bnlearn中惩罚值是这样的：

res[i] -= (*k) * nparams;

k是什么？
我现在知道的信息是：因为res减去了一个惩罚值，也就是bic，我从手边的一个论文A Transformational Characterization of Equivalent Bayesian Network Structures里看到了bic打分的公式就是这个样子的：
在这里插入图片描述
至于bic的公式1/2到底乘不乘先按下不表，bnlearn代码中的那个k一直没有交代清楚来源，只知道k是由参数extra.arg决定的，但是代码注释也没有给出k或者extra.arg的意义。疑惑。

2021.6.26，评论区补充：k: the penalty coefficient to be used by the AIC and BIC scores. The default value is 1 for AIC and log(nrow(data))/2 for BIC.他的R文档里是这么写的，文档地址https://www.bnlearn.com/documentation/man/score.html
感谢这位朋友

Rilkean heart

关注

0
点赞
踩
3

收藏

觉得还不错? 一键收藏
7
评论
15. R的bnlearn包中的per.node.score函数的学习探究

0. 写在前面原来只知道per.node.score是对贝叶斯网络单个节点进行打分，具体实现从来没有触及过，最近终于努力摁着头学习了一下。诸多不足，欢迎交流。此文比较适合有一定使用bnlearn包基础以及C语言基础并且也为这个函数的机理头疼的人观看本文使用的数据是一种n*n的基因表型数据，暂把它理解为列名和行名都表示基因，不理解数据不影响阅读。csdn竟然不支持R代码的高亮，故用pyth...
复制链接

扫一扫

专栏目录