In bnlearn, cpquery gives random probablities Running the example from bnlearn documentation in R
    data(gaussian.test)
    fitted = bn.fit(hc(gaussian.test), gaussian.test)
    # the result should be around 0.04.
    cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <=                 3)),evidence = (C + D < 10))

When I run the last command multiple times, I get different probabilities.
    cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)),evidence = (C+D<10))
    [1] 0.03792526
    cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)),evidence = (C+D<10))
    [1] 0.03990878
    cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)),evidence = (C+D<10))
    [1] 0.04542807
    cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)),evidence = (C+D<10))
    [1] 0.0353139
    cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)),evidence = (C+D<10))
    [1] 0.03755605
    cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)),evidence = (C+D<10))
    [1] 0.03730921
    cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)),evidence = (C+D<10))
    [1] 0.0467128
    cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)),evidence = (C+D<10))
    [1] 0.04154969
    cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)),evidence = (C+D<10))
    [1] 0.04973038

I am new to Bayesian Network, but shouldn't the probabilities be constant for an event. It is for continuous variable, so ideally its calculating the integrals form the joint probability functions. Its not like the gaussian data is generated every time a query is run. The data is generated in the first line and then reused. The manual says that the answer should be around 0.04, but I am also getting 0.049 and 0.029.
I am getting this on other queries to for my own data. So was wondering if anyone has an answer.
 A: From the help page of ?cpquery: states this : 

Note that both cpquery and cpdist are based on Monte Carlo particle filters, and therefore they may return slightly 
  different values on different runs.

You can reduce the variability in the inference runs by increasing the number of 
draws in the sampling procedure by using the tuning parameter, n.
So increase the number of draws by setting the tning parameter n. ie 
cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)),
                                        evidence = (C + D < 10), n=10000000)

To make subsequent runs reproducible use set.seed
ie 
set.seed(1)
cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)), 
                                        evidence = (C + D < 10), n=10000000)
#[1] 0.04071599

set.seed(1)
cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)), 
                                        evidence = (C + D < 10), n=10000000)
#[1] 0.04071599

So, by increasing n, the sampling procedure will run for longer, but will likely be more accurate
set.seed(1)
run1 <- replicate(100, 
          cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)), 
                                                   evidence = (C + D < 10)))
summary(run1)
#   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#0.03140 0.03751 0.04086 0.04103 0.04444 0.05387   

# Increase number of draws -> lower variation
set.seed(1)  
run2 <- replicate(100,  
          cpquery(fitted, event = ((A >= 0) & (A <= 1)) & ((B >= 0) & (B <= 3)), 
                                            evidence = (C + D < 10),  n=1000000))
summary(run2)
#    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
# 0.04006 0.04059 0.04081 0.04084 0.04103 0.04200   

Setting same axis-limits to see the lower variation

