Commit 411ff1ec authored by Gilles Kratzer's avatar Gilles Kratzer
Browse files

bug corrected with R3.6

parent 25bded24
Pipeline #1697 passed with stage
in 3 seconds
......@@ -142,6 +142,7 @@ mc3 <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, prior.choic
# alpha <- min(exp((score.Gprime + log(prior.Gprime) + log(n.Gprime)) - ((score.G + log(prior.G)) + log(n.G))),1)
alpha <- min(exp(score.Gprime - score.G) * (n.G/n.Gprime) * (prior.Gprime/prior.G), 1)
if(!is.numeric(alpha)) alpha <- 0
......
......@@ -30,13 +30,17 @@
if (is.numeric(prob.mbr) && prob.mbr > 1 && prob.mbr < 0)
stop("prob.mbr should be a probability.")
if (is.matrix(start.dag) && dim(start.dag) != c(length(data.dists), length(data.dists)))
if (is.matrix(start.dag) && (dim(start.dag)[1] != length(data.dists)))
stop("start.dag should be a squared matrix with dimension equal to the number of variables.")
if (is.matrix(start.dag) && is.null(topoSortMAT((start.dag)))) {
stop("start.dag should be a named DAG.")
}
if (is.matrix(start.dag) && max(rowSums(start.dag))>max.parents) {
stop("start.dag should not have more parent than max.parent argument.")
}
if (length(data.dists) != max(score.cache$children))
stop("data.dists should be a named list of all variables used to build the cache of precomputed score.")
......@@ -48,7 +52,7 @@
stop("prior.choice should be either 1,2 or 3.")
}
if (is.matrix(prior.dag) && dim(prior.dag) != c(length(data.dists), length(data.dists)))
if (is.matrix(prior.dag) && dim(prior.dag)[1] != length(data.dists))
stop("prior.dag should be a squared matrix with dimension equal to the number of variables.")
}
......@@ -186,7 +190,7 @@ formula.mcmcabn <- function(f, name) {
descendents <- function(nodes, dag) {
if (!identical(topoSortMAT(amat = dag, index = FALSE), character(0))) {
if (!is.integer0(which(x = dag[, nodes] == 1, arr.ind = TRUE)) && diag(dag) == 0) {
if (!is.integer0(which(x = dag[, nodes] == 1, arr.ind = TRUE)) && sum(diag(dag)) == 0) {
if (length(nodes) == 1) {
tmp <- unname(which(x = dag[, nodes] == 1, arr.ind = TRUE))
} else {
......
......@@ -142,11 +142,16 @@
<h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2>
<pre class="examples"><span class='co'># NOT RUN {</span>
<span class='co'>## This data set was generated using the following code:</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>library</a></span>(<span class='no'>bnlearn</span>) <span class='co'>#for the dataset</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>library</a></span>(<span class='no'>abn</span>) <span class='co'>#for the cache of score function</span>
<pre class="examples"><div class='input'><span class='co'>## This data set was generated using the following code:</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>library</a></span>(<span class='no'>bnlearn</span>) <span class='co'>#for the dataset</span></div><div class='output co'>#&gt; <span class='message'></span>
#&gt; <span class='message'>Attaching package: ‘bnlearn’</span></div><div class='output co'>#&gt; <span class='message'>The following object is masked from ‘package:testthat’:</span>
#&gt; <span class='message'></span>
#&gt; <span class='message'> compare</span></div><div class='output co'>#&gt; <span class='message'>The following object is masked from ‘package:stats’:</span>
#&gt; <span class='message'></span>
#&gt; <span class='message'> sigma</span></div><div class='input'><span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>library</a></span>(<span class='no'>abn</span>) <span class='co'>#for the cache of score function</span></div><div class='output co'>#&gt; <span class='message'>Loading required package: nnet</span></div><div class='output co'>#&gt; <span class='message'>Loading required package: Cairo</span></div><div class='output co'>#&gt; <span class='message'>Loading required package: MASS</span></div><div class='output co'>#&gt; <span class='message'>Loading required package: lme4</span></div><div class='output co'>#&gt; <span class='warning'>Warning: package ‘lme4’ was built under R version 3.5.2</span></div><div class='output co'>#&gt; <span class='message'>Loading required package: Matrix</span></div><div class='output co'>#&gt; <span class='message'></span>
#&gt; <span class='message'>Attaching package: ‘abn’</span></div><div class='output co'>#&gt; <span class='message'>The following object is masked from ‘package:bnlearn’:</span>
#&gt; <span class='message'></span>
#&gt; <span class='message'> mb</span></div><div class='input'>
<span class='co'>#renaming columns of the dataset</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/colnames'>colnames</a></span>(<span class='no'>asia</span>) <span class='kw'>&lt;-</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/c'>c</a></span>(<span class='st'>"Asia"</span>,
<span class='st'>"Smoking"</span>,
......@@ -169,20 +174,7 @@
<span class='no'>bsc.compute.asia</span> <span class='kw'>&lt;-</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/abn/topics/build_score_cache'>buildscorecache</a></span>(<span class='kw'>data.df</span> <span class='kw'>=</span> <span class='no'>asia</span>,
<span class='kw'>data.dists</span> <span class='kw'>=</span> <span class='no'>dist.asia</span>,
<span class='kw'>max.parents</span> <span class='kw'>=</span> <span class='fl'>2</span>)
<span class='no'>mcmc.out.asia</span> <span class='kw'>&lt;-</span> <span class='fu'><a href='mcmc.html'>mcmcabn</a></span>(<span class='kw'>score.cache</span> <span class='kw'>=</span> <span class='no'>bsc.compute.asia</span>,
<span class='kw'>score</span> <span class='kw'>=</span> <span class='st'>"mlik"</span>,
<span class='kw'>data.dists</span> <span class='kw'>=</span> <span class='no'>dist.asia</span>,
<span class='kw'>max.parents</span> <span class='kw'>=</span> <span class='fl'>2</span>,
<span class='kw'>mcmc.scheme</span> <span class='kw'>=</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/c'>c</a></span>(<span class='fl'>1000</span>,<span class='fl'>99</span>,<span class='fl'>1000</span>),
<span class='kw'>seed</span> <span class='kw'>=</span> <span class='fl'>42</span>,
<span class='kw'>verbose</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>,
<span class='kw'>start.dag</span> <span class='kw'>=</span> <span class='st'>"random"</span>,
<span class='kw'>prob.rev</span> <span class='kw'>=</span> <span class='fl'>0.03</span>,
<span class='kw'>prob.mbr</span> <span class='kw'>=</span> <span class='fl'>0.03</span>,
<span class='kw'>prior.choice</span> <span class='kw'>=</span> <span class='fl'>2</span>)
<span class='co'># }</span></pre>
<span class='kw'>max.parents</span> <span class='kw'>=</span> <span class='fl'>2</span>)</div></pre>
</div>
<div class="col-md-3 hidden-xs hidden-sm" id="sidebar">
<h2>Contents</h2>
......
......@@ -142,8 +142,7 @@
<h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2>
<pre class="examples"><span class='co'># NOT RUN {</span>
<span class='co'>## This data set was generated using the following code:</span>
<pre class="examples"><div class='input'><span class='co'>## This data set was generated using the following code:</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>library</a></span>(<span class='no'>bnlearn</span>) <span class='co'>#for the dataset</span>
<span class='co'>#renaming columns of the dataset</span>
......@@ -164,8 +163,7 @@
<span class='kw'>Bronchitis</span> <span class='kw'>=</span> <span class='st'>"binomial"</span>,
<span class='kw'>Either</span> <span class='kw'>=</span> <span class='st'>"binomial"</span>,
<span class='kw'>XRay</span> <span class='kw'>=</span> <span class='st'>"binomial"</span>,
<span class='kw'>Dyspnea</span> <span class='kw'>=</span> <span class='st'>"binomial"</span>)
<span class='co'># }</span></pre>
<span class='kw'>Dyspnea</span> <span class='kw'>=</span> <span class='st'>"binomial"</span>)</div></pre>
</div>
<div class="col-md-3 hidden-xs hidden-sm" id="sidebar">
<h2>Contents</h2>
......
......@@ -142,8 +142,7 @@
<h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2>
<pre class="examples"><span class='co'># NOT RUN {</span>
<span class='co'>## This data set was generated using the following code:</span>
<pre class="examples"><div class='input'><span class='co'>## This data set was generated using the following code:</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>library</a></span>(<span class='no'>bnlearn</span>) <span class='co'>#for the dataset</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>library</a></span>(<span class='no'>abn</span>) <span class='co'>#for the cache of scores computing function</span>
......@@ -181,8 +180,7 @@
<span class='kw'>start.dag</span> <span class='kw'>=</span> <span class='st'>"random"</span>,
<span class='kw'>prob.rev</span> <span class='kw'>=</span> <span class='fl'>0.03</span>,
<span class='kw'>prob.mbr</span> <span class='kw'>=</span> <span class='fl'>0.03</span>,
<span class='kw'>prior.choice</span> <span class='kw'>=</span> <span class='fl'>2</span>)
<span class='co'># }</span></pre>
<span class='kw'>prior.choice</span> <span class='kw'>=</span> <span class='fl'>2</span>)</div></pre>
</div>
<div class="col-md-3 hidden-xs hidden-sm" id="sidebar">
<h2>Contents</h2>
......
......@@ -234,25 +234,194 @@
<h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2>
<pre class="examples"><span class='co'># NOT RUN {</span>
<span class='co'>## Example from the asia dataset from Lauritzen and Spiegelhalter (1988) provided by Scutari (2010)</span>
<pre class="examples"><div class='input'><span class='co'>## Example from the asia dataset from Lauritzen and Spiegelhalter (1988) provided by Scutari (2010)</span>
<span class='co'># the number of MCMC run is delibaretelly chosen too small (computing time)</span>
<span class='co'># no thinning (usually not recommended)</span>
<span class='co'># no burn-in (usually not recommended,</span>
<span class='co'># even if not supported by any theoretical arguments)</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/utils/topics/data'>data</a></span>(<span class='st'>"mcmc_run_asia"</span>)
<span class='no'>mcmc.out.asia</span> <span class='kw'>&lt;-</span> <span class='fu'>mcmcabn</span>(<span class='kw'>score.cache</span> <span class='kw'>=</span> <span class='no'>bsc.compute.asia</span>,
<span class='co'># let us run: 0.03 REV, 0.03 MBR, 0.94 MC3 MCMC jumps</span>
<span class='co'># with a random DAG as starting point</span>
<span class='no'>mcmc.out.asia.small</span> <span class='kw'>&lt;-</span> <span class='fu'>mcmcabn</span>(<span class='kw'>score.cache</span> <span class='kw'>=</span> <span class='no'>bsc.compute.asia</span>,
<span class='kw'>score</span> <span class='kw'>=</span> <span class='st'>"mlik"</span>,
<span class='kw'>data.dists</span> <span class='kw'>=</span> <span class='no'>dist.asia</span>,
<span class='kw'>max.parents</span> <span class='kw'>=</span> <span class='fl'>2</span>,
<span class='kw'>mcmc.scheme</span> <span class='kw'>=</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/c'>c</a></span>(<span class='fl'>1000</span>,<span class='fl'>99</span>,<span class='fl'>10000</span>),
<span class='kw'>seed</span> <span class='kw'>=</span> <span class='fl'>42</span>,
<span class='kw'>mcmc.scheme</span> <span class='kw'>=</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/c'>c</a></span>(<span class='fl'>100</span>,<span class='fl'>0</span>,<span class='fl'>0</span>),
<span class='kw'>seed</span> <span class='kw'>=</span> <span class='fl'>321</span>,
<span class='kw'>verbose</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>,
<span class='kw'>start.dag</span> <span class='kw'>=</span> <span class='st'>"random"</span>,
<span class='kw'>prob.rev</span> <span class='kw'>=</span> <span class='fl'>0.03</span>,
<span class='kw'>prob.mbr</span> <span class='kw'>=</span> <span class='fl'>0.03</span>,
<span class='kw'>prior.choice</span> <span class='kw'>=</span> <span class='fl'>2</span>)
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/summary'>summary</a></span>(<span class='no'>mcmc.out.asia</span>)
<span class='co'># }</span></pre>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/summary'>summary</a></span>(<span class='no'>mcmc.out.asia.small</span>)</div><div class='output co'>#&gt; MCMC summary:
#&gt; Number of Burn in steps: 0
#&gt; Number of MCMC steps: 100
#&gt; Thinning: 0
#&gt;
#&gt; Maximum score: -11370.47
#&gt; Empirical mean: -12614.7
#&gt; Empirical standard deviation: 1392.967
#&gt; Quantiles of the posterior network score:
#&gt; 0.025 0.25 0.5 0.75 0.975
#&gt; BN score -15127.73 -13458.25 -12044.84 -11371.57 -11370.47
#&gt;
#&gt;
#&gt; Global acceptance rate: 0.3168317
#&gt; Accepted Rejected
#&gt; MBR 0 2
#&gt; MC3 30 67
#&gt; REV 2 0
#&gt;
#&gt;
#&gt; Sample size adjusted for autocorrelation: 1.406361
#&gt;
#&gt; Autocorrelations by lag:
#&gt; 0 1 2 3 4 5 6 7
#&gt; acf 1 0.9722601 0.9437878 0.9152083 0.8843368 0.8535088 0.8225977 0.7915496
#&gt; 8 9 10
#&gt; acf 0.7605535 0.7286543 0.6975366</div><div class='input'>
<span class='co'># Uniquelly with MC3 moves</span>
<span class='no'>mcmc.out.asia.small</span> <span class='kw'>&lt;-</span> <span class='fu'>mcmcabn</span>(<span class='kw'>score.cache</span> <span class='kw'>=</span> <span class='no'>bsc.compute.asia</span>,
<span class='kw'>score</span> <span class='kw'>=</span> <span class='st'>"mlik"</span>,
<span class='kw'>data.dists</span> <span class='kw'>=</span> <span class='no'>dist.asia</span>,
<span class='kw'>max.parents</span> <span class='kw'>=</span> <span class='fl'>2</span>,
<span class='kw'>mcmc.scheme</span> <span class='kw'>=</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/c'>c</a></span>(<span class='fl'>100</span>,<span class='fl'>0</span>,<span class='fl'>0</span>),
<span class='kw'>seed</span> <span class='kw'>=</span> <span class='fl'>42</span>,
<span class='kw'>verbose</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>,
<span class='kw'>start.dag</span> <span class='kw'>=</span> <span class='st'>"random"</span>,
<span class='kw'>prob.rev</span> <span class='kw'>=</span> <span class='fl'>0</span>,
<span class='kw'>prob.mbr</span> <span class='kw'>=</span> <span class='fl'>0</span>,
<span class='kw'>prior.choice</span> <span class='kw'>=</span> <span class='fl'>2</span>)
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/summary'>summary</a></span>(<span class='no'>mcmc.out.asia.small</span>)</div><div class='output co'>#&gt; MCMC summary:
#&gt; Number of Burn in steps: 0
#&gt; Number of MCMC steps: 100
#&gt; Thinning: 0
#&gt;
#&gt; Maximum score: -13046.48
#&gt; Empirical mean: -13667.21
#&gt; Empirical standard deviation: 777.3341
#&gt; Quantiles of the posterior network score:
#&gt; 0.025 0.25 0.5 0.75 0.975
#&gt; BN score -15147.68 -14640.66 -13136.8 -13136.28 -13046.48
#&gt;
#&gt;
#&gt; Global acceptance rate: 0.3960396
#&gt; Accepted Rejected
#&gt; MC3 40 61
#&gt;
#&gt;
#&gt; Sample size adjusted for autocorrelation: 1.732075
#&gt;
#&gt; Autocorrelations by lag:
#&gt; 0 1 2 3 4 5 6 7
#&gt; acf 1 0.965945 0.9316219 0.8967652 0.8616907 0.8313675 0.8032237 0.7761743
#&gt; 8 9 10
#&gt; acf 0.7491203 0.7197729 0.6902472</div><div class='input'>
<span class='co'>#let us define a starting DAG</span>
<span class='no'>startDag</span> <span class='kw'>&lt;-</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/matrix'>matrix</a></span>(<span class='kw'>data</span> <span class='kw'>=</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/c'>c</a></span>(<span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>1</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>,
<span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>1</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>,
<span class='fl'>1</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>,
<span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>1</span>, <span class='fl'>0</span>, <span class='fl'>0</span>,
<span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>,
<span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>,
<span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>,
<span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>, <span class='fl'>0</span>),<span class='kw'>nrow</span> <span class='kw'>=</span> <span class='fl'>8</span>,<span class='kw'>ncol</span> <span class='kw'>=</span> <span class='fl'>8</span>, <span class='kw'>byrow</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>)
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/colnames'>colnames</a></span>(<span class='no'>startDag</span>) <span class='kw'>&lt;-</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/colnames'>rownames</a></span>(<span class='no'>startDag</span>) <span class='kw'>&lt;-</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/names'>names</a></span>(<span class='no'>dist.asia</span>)
<span class='co'># Additionally, let us use the non informative prior</span>
<span class='no'>mcmc.out.asia.small</span> <span class='kw'>&lt;-</span> <span class='fu'>mcmcabn</span>(<span class='kw'>score.cache</span> <span class='kw'>=</span> <span class='no'>bsc.compute.asia</span>,
<span class='kw'>score</span> <span class='kw'>=</span> <span class='st'>"mlik"</span>,
<span class='kw'>data.dists</span> <span class='kw'>=</span> <span class='no'>dist.asia</span>,
<span class='kw'>max.parents</span> <span class='kw'>=</span> <span class='fl'>2</span>,
<span class='kw'>mcmc.scheme</span> <span class='kw'>=</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/c'>c</a></span>(<span class='fl'>100</span>,<span class='fl'>0</span>,<span class='fl'>0</span>),
<span class='kw'>seed</span> <span class='kw'>=</span> <span class='fl'>42</span>,
<span class='kw'>verbose</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>,
<span class='kw'>start.dag</span> <span class='kw'>=</span> <span class='no'>startDag</span>,
<span class='kw'>prob.rev</span> <span class='kw'>=</span> <span class='fl'>0</span>,
<span class='kw'>prob.mbr</span> <span class='kw'>=</span> <span class='fl'>0</span>,
<span class='kw'>prior.choice</span> <span class='kw'>=</span> <span class='fl'>1</span>)
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/summary'>summary</a></span>(<span class='no'>mcmc.out.asia.small</span>)</div><div class='output co'>#&gt; MCMC summary:
#&gt; Number of Burn in steps: 0
#&gt; Number of MCMC steps: 100
#&gt; Thinning: 0
#&gt;
#&gt; Maximum score: -12351.24
#&gt; Empirical mean: -12726.16
#&gt; Empirical standard deviation: 574.7013
#&gt; Quantiles of the posterior network score:
#&gt; 0.025 0.25 0.5 0.75 0.975
#&gt; BN score -14160.15 -12492.22 -12461.74 -12449.72 -12351.25
#&gt;
#&gt;
#&gt; Global acceptance rate: 0.3960396
#&gt; Accepted Rejected
#&gt; MC3 40 61
#&gt;
#&gt;
#&gt; Sample size adjusted for autocorrelation: 2.557271
#&gt;
#&gt; Autocorrelations by lag:
#&gt; 0 1 2 3 4 5 6 7
#&gt; acf 1 0.950125 0.8971036 0.8441016 0.7908641 0.7417565 0.6928941 0.6438907
#&gt; 8 9 10
#&gt; acf 0.5983815 0.5505099 0.5010223</div><div class='input'>
<span class='co'># let us define our very own prior</span>
<span class='co'># we know that there should be a link between Smoking and LungCancer nodes</span>
<span class='co'># uninformative prior matrix</span>
<span class='no'>priorDag</span> <span class='kw'>&lt;-</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/matrix'>matrix</a></span>(<span class='kw'>data</span> <span class='kw'>=</span> <span class='fl'>0.5</span>,<span class='kw'>nrow</span> <span class='kw'>=</span> <span class='fl'>8</span>,<span class='kw'>ncol</span> <span class='kw'>=</span> <span class='fl'>8</span>)
<span class='co'># name it</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/colnames'>colnames</a></span>(<span class='no'>priorDag</span>) <span class='kw'>&lt;-</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/colnames'>rownames</a></span>(<span class='no'>priorDag</span>) <span class='kw'>&lt;-</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/names'>names</a></span>(<span class='no'>dist.asia</span>)
<span class='co'># parent = smoking; child = LungCancer</span>
<span class='no'>priorDag</span>[<span class='st'>"LungCancer"</span>,<span class='st'>"Smoking"</span>] <span class='kw'>&lt;-</span> <span class='fl'>1</span>
<span class='no'>mcmc.out.asia.small</span> <span class='kw'>&lt;-</span> <span class='fu'>mcmcabn</span>(<span class='kw'>score.cache</span> <span class='kw'>=</span> <span class='no'>bsc.compute.asia</span>,
<span class='kw'>score</span> <span class='kw'>=</span> <span class='st'>"mlik"</span>,
<span class='kw'>data.dists</span> <span class='kw'>=</span> <span class='no'>dist.asia</span>,
<span class='kw'>max.parents</span> <span class='kw'>=</span> <span class='fl'>2</span>,
<span class='kw'>mcmc.scheme</span> <span class='kw'>=</span> <span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/c'>c</a></span>(<span class='fl'>100</span>,<span class='fl'>0</span>,<span class='fl'>0</span>),
<span class='kw'>seed</span> <span class='kw'>=</span> <span class='fl'>42</span>,
<span class='kw'>verbose</span> <span class='kw'>=</span> <span class='fl'>FALSE</span>,
<span class='kw'>start.dag</span> <span class='kw'>=</span> <span class='no'>startDag</span>,
<span class='kw'>prob.rev</span> <span class='kw'>=</span> <span class='fl'>0</span>,
<span class='kw'>prob.mbr</span> <span class='kw'>=</span> <span class='fl'>0</span>,
<span class='kw'>prior.choice</span> <span class='kw'>=</span> <span class='fl'>3</span>,
<span class='kw'>prior.dag</span> <span class='kw'>=</span> <span class='no'>priorDag</span>)
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/summary'>summary</a></span>(<span class='no'>mcmc.out.asia.small</span>)</div><div class='output co'>#&gt; MCMC summary:
#&gt; Number of Burn in steps: 0
#&gt; Number of MCMC steps: 100
#&gt; Thinning: 0
#&gt;
#&gt; Maximum score: -12351.24
#&gt; Empirical mean: -12726.16
#&gt; Empirical standard deviation: 574.7013
#&gt; Quantiles of the posterior network score:
#&gt; 0.025 0.25 0.5 0.75 0.975
#&gt; BN score -14160.15 -12492.22 -12461.74 -12449.72 -12351.25
#&gt;
#&gt;
#&gt; Global acceptance rate: 0.3960396
#&gt; Accepted Rejected
#&gt; MC3 40 61
#&gt;
#&gt;
#&gt; Sample size adjusted for autocorrelation: 2.557271
#&gt;
#&gt; Autocorrelations by lag:
#&gt; 0 1 2 3 4 5 6 7
#&gt; acf 1 0.950125 0.8971036 0.8441016 0.7908641 0.7417565 0.6928941 0.6438907
#&gt; 8 9 10
#&gt; acf 0.5983815 0.5505099 0.5010223</div></pre>
</div>
<div class="col-md-3 hidden-xs hidden-sm" id="sidebar">
<h2>Contents</h2>
......
......@@ -145,8 +145,7 @@
<h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2>
<pre class="examples"><span class='co'># NOT RUN {</span>
<span class='co'>## This data set was generated using the following code:</span>
<pre class="examples"><div class='input'><span class='co'>## This data set was generated using the following code:</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>library</a></span>(<span class='no'>bnlearn</span>) <span class='co'>#for the dataset</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/library'>library</a></span>(<span class='no'>abn</span>) <span class='co'>#for the cache of score function</span>
......@@ -184,8 +183,7 @@
<span class='kw'>start.dag</span> <span class='kw'>=</span> <span class='st'>"random"</span>,
<span class='kw'>prob.rev</span> <span class='kw'>=</span> <span class='fl'>0.03</span>,
<span class='kw'>prob.mbr</span> <span class='kw'>=</span> <span class='fl'>0.03</span>,
<span class='kw'>prior.choice</span> <span class='kw'>=</span> <span class='fl'>2</span>)
<span class='co'># }</span></pre>
<span class='kw'>prior.choice</span> <span class='kw'>=</span> <span class='fl'>2</span>)</div></pre>
</div>
<div class="col-md-3 hidden-xs hidden-sm" id="sidebar">
<h2>Contents</h2>
......
......@@ -170,13 +170,13 @@ Scutari, M. (2010). Learning Bayesian Networks with the bnlearn R Package. Journ
<h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2>
<pre class="examples"><span class='co'># NOT RUN {</span>
<span class='co'>## Example from the asia dataset from Lauritzen and Spiegelhalter (1988) provided by Scutari (2010)</span>
<pre class="examples"><div class='input'><span class='co'>## Example from the asia dataset from Lauritzen and Spiegelhalter (1988) provided by Scutari (2010)</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/utils/topics/data'>data</a></span>(<span class='st'>"mcmc_run_asia"</span>)
<span class='co'>#plot the mcmc run</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/graphics/topics/plot'>plot</a></span>(<span class='no'>mcmc.out.asia</span>)
<span class='co'># }</span></pre>
<span class='fu'><a href='https://www.rdocumentation.org/packages/graphics/topics/plot'>plot</a></span>(<span class='no'>mcmc.out.asia</span>)</div><div class='img'><img src='plot-1.png' alt='' width='700' height='433' /></div><div class='input'>
<span class='co'>#plot cumulative max score</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/graphics/topics/plot'>plot</a></span>(<span class='no'>mcmc.out.asia</span>, <span class='kw'>max.score</span> <span class='kw'>=</span> <span class='fl'>TRUE</span>)</div><div class='img'><img src='plot-2.png' alt='' width='700' height='433' /></div></pre>
</div>
<div class="col-md-3 hidden-xs hidden-sm" id="sidebar">
<h2>Contents</h2>
......
......@@ -158,13 +158,12 @@ print(x, &#8230;)</pre>
<h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2>
<pre class="examples"><span class='co'># NOT RUN {</span>
<span class='co'>## Example from the asia dataset from Lauritzen and Spiegelhalter (1988) provided by Scutari (2010)</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/utils/topics/data'>data</a></span>(<span class='st'>"mcmc_run_asia"</span>)
<pre class="examples"><div class='input'><span class='co'>## Example from the asia dataset from Lauritzen and Spiegelhalter (1988) provided by Scutari (2010)</span>
<span class='co'>#print the MCMC run</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/print'>print</a></span>(<span class='no'>mcmc.out.asia</span>)
<span class='co'># }</span></pre>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/print'>print</a></span>(<span class='no'>mcmc.out.asia</span>)</div><div class='output co'>#&gt; Posterior Bayesian network score estimated using MCMC:Number of Burn in steps: 10000
#&gt; Number of MCMC steps: 1e+05
#&gt; Thinning: 99
#&gt; </div></pre>
</div>
<div class="col-md-3 hidden-xs hidden-sm" id="sidebar">
<h2>Contents</h2>
......
......@@ -173,13 +173,35 @@ summary(object,
<h2 class="hasAnchor" id="examples"><a class="anchor" href="#examples"></a>Examples</h2>
<pre class="examples"><span class='co'># NOT RUN {</span>
<span class='co'>## Example from the asia dataset from Lauritzen and Spiegelhalter (1988) provided by Scutari (2010)</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/utils/topics/data'>data</a></span>(<span class='st'>"mcmc_run_asia"</span>)
<pre class="examples"><div class='input'><span class='co'>## Example from the asia dataset from Lauritzen and Spiegelhalter (1988) provided by Scutari (2010)</span>
<span class='co'>#summary the MCMC run</span>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/summary'>summary</a></span>(<span class='no'>mcmc.out.asia</span>)
<span class='co'># }</span></pre>
<span class='fu'><a href='https://www.rdocumentation.org/packages/base/topics/summary'>summary</a></span>(<span class='no'>mcmc.out.asia</span>)</div><div class='output co'>#&gt; MCMC summary:
#&gt; Number of Burn in steps: 10000
#&gt; Number of MCMC steps: 1e+05
#&gt; Thinning: 99
#&gt;
#&gt; Maximum score: -11151.14
#&gt; Empirical mean: -11696.06
#&gt; Empirical standard deviation: 544.0571
#&gt; Quantiles of the posterior network score:
#&gt; 0.025 0.25 0.5 0.75 0.975
#&gt; BN score -13198.51 -12149.93 -11499.1 -11271.92 -11156.35
#&gt;
#&gt;
#&gt; Global acceptance rate: 0.2577423
#&gt; Accepted Rejected
#&gt; MBR 0 32
#&gt; MC3 247 698
#&gt; REV 11 13
#&gt;
#&gt;
#&gt; Sample size adjusted for autocorrelation: 380.8894
#&gt;
#&gt; Autocorrelations by lag:
#&gt; 0 1 2 3 4 5 6
#&gt; acf 1 0.3996997 0.2091346 0.07279787 0.04835775 0.03308019 0.00552735
#&gt; 7 8 9 10
#&gt; acf 0.0004674749 -0.01328015 -0.02254146 -0.01459581</div></pre>
</div>
<div class="col-md-3 hidden-xs hidden-sm" id="sidebar">
<h2>Contents</h2>
......
......@@ -108,7 +108,7 @@ mcmc.out.asia.small <- mcmcabn(score.cache = bsc.compute.asia,
data.dists = dist.asia,
max.parents = 2,
mcmc.scheme = c(100,0,0),
seed = 42,
seed = 321,
verbose = FALSE,
start.dag = "random",
prob.rev = 0.03,
......@@ -133,12 +133,19 @@ mcmc.out.asia.small <- mcmcabn(score.cache = bsc.compute.asia,
summary(mcmc.out.asia.small)
#let us define a starting DAG (empty matrix = no arcs)
startDag <- matrix(data = 0,nrow = 8,ncol = 8)
#name it
#let us define a starting DAG
startDag <- matrix(data = c(0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0,
1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0),nrow = 8,ncol = 8, byrow = TRUE)
colnames(startDag) <- rownames(startDag) <- names(dist.asia)
# Additionally, let us use the non informative prior
# Additionally, let us use the non informative prior
mcmc.out.asia.small <- mcmcabn(score.cache = bsc.compute.asia,
score = "mlik",
data.dists = dist.asia,
......@@ -156,8 +163,8 @@ summary(mcmc.out.asia.small)
# let us define our very own prior
# we know that there should be a link between Smoking and LungCancer nodes
# empty matrix
priorDag <- matrix(data = 0,nrow = 8,ncol = 8)
# uninformative prior matrix
priorDag <- matrix(data = 0.5,nrow = 8,ncol = 8)
# name it
colnames(priorDag) <- rownames(priorDag) <- names(dist.asia)
# parent = smoking; child = LungCancer
......
No preview for this file type
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment