1

I am trying to reproduce a network i have previously created in GeNie, but i get an error that some conditional probability distributions of node probability do not sum to one.

   # create empty BN
    HazardNet <- empty.graph(nodes = c("scenario", "probability", "intensity", "hazard"))
    HazardNet

# dependencies
HazardNet <- set.arc(HazardNet, from = "scenario", to = "probability")
HazardNet <- set.arc(HazardNet, from = "scenario", to = "intensity")
HazardNet <- set.arc(HazardNet, from = "intensity", to = "hazard")
HazardNet <- set.arc(HazardNet, from = "probability", to = "hazard")
HazardNet

# define states - global distribution
lmh <- c("low", "moderate", "high") 
scenario.lv <- c("30 year", "100 year", "300year", "none")
intensity.lv <- lmh
probability.lv <- lmh
hazard.lv <- lmh

# create conditional probability tables
scenario.prob <- array(c(0.01, 0.03, 0.05, 0.91), dim = 4,
                        dimnames = list(scenario = scenario.lv))
probability.prob <- array(c(rep(1/3, 12)), dim = c(4,3),
                          dimnames = list(scenario =scenario.lv, probability = probability.lv))
intensity.prob <- array(c(0, 0.5, 0.5, 0, 0.3, 0.7, 0, 0.1, 0.9, 0.98, 0.019, 0.001), dim = c(3, 4),
                          dimnames = list(intensity = intensity.lv, scenario = scenario.lv ))
hazard.prob <- array(c(0, 1, 0, 0.5, 0.5, 0, 1, 0, 0, 0.5, 0.5, 0, 0, 1, 0, 0, 0.5, 0.5, 0, 0, 1, 0, 0, 1, 0, 0, 1 ), 
                     dim = c(3, 3, 4), 
                     dimnames = list( intensity = intensity.lv, hazard= hazard.lv, scenario = scenario.lv))

cpt <- list(scenario = scenario.prob, probability = probability.prob, intensity = intensity.prob, hazard = hazard.lv)
bn <- custom.fit(HazardNet, cpt)

I assume the problem is the node scenario? If I see right all my other cpt are fine.

How can i code equal probabilities for a node with 3 states? 

> scenario.prob
scenario
 30 year 100 year  300year     none 
    0.01     0.03     0.05     0.91 
> probability.prob
          probability
scenario         low  moderate      high
  30 year  0.3333333 0.3333333 0.3333333
  100 year 0.3333333 0.3333333 0.3333333
  300year  0.3333333 0.3333333 0.3333333
  none     0.3333333 0.3333333 0.3333333
> intensity.prob
          scenario
intensity  30 year 100 year 300year  none
  low          0.0      0.0     0.0 0.980
  moderate     0.5      0.3     0.1 0.019
  high         0.5      0.7     0.9 0.001
> hazard.prob
, , scenario = 30 year

          hazard
intensity  low moderate high
  low        0      0.5    1
  moderate   1      0.5    0
  high       0      0.0    0

, , scenario = 100 year

          hazard
intensity  low moderate high
  low      0.5        0  0.0
  moderate 0.5        1  0.5
  high     0.0        0  0.5

, , scenario = 300year

          hazard
intensity  low moderate high
  low        0        0    0
  moderate   0        0    0
  high       1        1    1

, , scenario = none

          hazard
intensity  low moderate high
  low        0      0.5    1
  moderate   1      0.5    0
  high       0      0.0    0
user1607
  • 531
  • 7
  • 28
  • there are a couple of errors. `hazard.prob` should have dims (3,3,3) as its parents are `intensity` and `probability` (not `scenario`). For the tables, IIRC the child node should be on the first dimension, so `probability.prob` will also need updated. – user20650 May 01 '18 at 10:29
  • thanks, that indeed worked... however when i set evidence using `jintensity <- setEvidence(junction, nodes = "intensity", states = "low" )` the conditional probability remains the same as the marginal one for node `intensity`... am i doing a mistake somewhere? – user1607 May 02 '18 at 12:53
  • Is `setEvidence` in `bnlearn`? - if not, please update your question with all code you have used. But if you set the state in a variable you would expect it to be one in the state of the marginal of the same node. (ps ways to get marginals in `bnlearn`: for prior marginal of `intensity`, `x = cpdist(bn,nodes="intensity" , evidence = TRUE, method="ls", n=1e5) ; prop.table(table(x))` or in this simple case you can calculate manually `tcrossprod(bn$scenario$prob, bn$intensity$prob)` (there may be simpler ways) – user20650 May 02 '18 at 15:09

0 Answers0