1

我是 R 的相对初学者，并试图弄清楚如何将 cpquery 函数用于 DAG 的所有边缘的 bnlearn 包。

首先，我创建了一个bn对象、一个bn网络和一个具有所有优势的表。

library(bnlearn)
data(learning.test)
baynet = hc(learning.test)
fit = bn.fit(baynet, learning.test)
sttbl = arc.strength(x = baynet, data = learning.test)

然后我尝试在sttbl数据集中创建一个新变量，这是 cpquery 函数的结果。

sttbl = sttbl %>% mutate(prob = NA) %>% arrange(strength)
sttbl[1,4] = cpquery(fit, `A` == 1, `D` == 1)

它看起来相当不错（尤其是在更大的数据上），但是当我试图以某种方式自动化这个过程时，我正在努力解决错误，例如：

采样错误（拟合 = 拟合，事件 = 事件，证据 = 证据，：证据的逻辑向量长度为​​ 1 而不是 10000。

在完美的情况下，我需要创建一个函数来填充sttbl数据集的概率生成变量，而不管它的大小。我试图用for循环来做，但一次又一次地偶然发现上面的错误。不幸的是，我正在删除失败的尝试，但它们是这样的：

for (i in 1:nrow(sttbl)) {
     j = sttbl[i,1]
     k = sttbl[i,2]
     sttbl[i,4]=cpquery(fit, fit$j %in% sttbl[i,1]==1, fit$k %in% sttbl[i,2]==1)
}

或这个：

for (i in 1:nrow(sttbl)) {
     sttbl[i,4]=cpquery(fit, sttbl[i,1] == 1, sttbl[i,2] == 1)
}

现在我想我误解了 R 或 bnlearn 包中的一些东西。

您能否告诉我如何通过多个 cpqueries 填充列来实现此任务？这对我的研究有很大帮助！

4

1 回答 1

1

cpquery很难以编程方式使用。如果您查看帮助页面中的示例，您可以看到作者用于eval(parse(...))构建查询。我在下面添加了两种方法，一种使用帮助页面中的方法，另一种cpdist用于绘制样本并重新加权以获得概率。

你的例子

library(bnlearn); library(dplyr)
data(learning.test)
baynet = hc(learning.test)
fit = bn.fit(baynet, learning.test)
sttbl = arc.strength(x = baynet, data = learning.test)
sttbl = sttbl %>% mutate(prob = NA) %>% arrange(strength)

这使用cpquery和备受诟病的eval(parse(...))- 这是bnlearn作者在?cpquery示例中以编程方式执行此操作的方法。反正，

# You want the evidence and event to be the same; in your question it is `1`
# but for example using learning.test data we use 'a'
state = "\'a\'" # note if the states are character then these need to be quoted
event = paste(sttbl$from, "==", state)
evidence = paste(sttbl$to, "==", state)

# loop through using code similar to that found in `cpquery`
set.seed(1) # to make sampling reproducible
for(i in 1:nrow(sttbl)) {
  qtxt = paste("cpquery(fit, ", event[i], ", ", evidence[i], ",n=1e6", ")")
  sttbl$prob[i] = eval(parse(text=qtxt))
}

我发现最好使用cpdistwhich 来根据某些证据生成随机样本。然后，您可以使用这些示例来构建查询。如果您使用似然加权 ( method="lw")，以编程方式执行此操作会稍微容易一些（并且不使用evil(parse(...))）。证据被添加到命名列表中，即list(A='a').

# The following just gives a quick way to assign the same
# evidence state to all the evidence nodes.  
evidence = setNames(replicate(nrow(sttbl), "a", simplify = FALSE), sttbl$to)

# Now loop though the queries
# As we are using likelihood weighting we need to reweight to get the probabilities
# (cpquery does this under the hood)
# Also note with this method that you could simulate from more than
# one variable (event) at a time if the evidence was the same.
for(i in 1:nrow(sttbl)) {
  temp = cpdist(fit, sttbl$from[i], evidence[i], method="lw")
  w = attr(temp, "weights")  
  sttbl$prob2[i] = sum(w[temp=='a'])/ sum(w)
}

sttbl
#   from to   strength      prob     prob2
# 1    A  D -1938.9499 0.6186238 0.6233387
# 2    A  B -1153.8796 0.6050552 0.6133448
# 3    C  D  -823.7605 0.7027782 0.7067417
# 4    B  E  -720.8266 0.7332107 0.7328657
# 5    F  E  -549.2300 0.5850828 0.5895373
于 2020-06-01T15:33:58.350 回答