Commit 96757e50 authored by Gilles Kratzer's avatar Gilles Kratzer
Browse files

update: print method summary, examples and dist.asia

parent b0fd386f
Pipeline #2546 passed with stage
in 3 seconds
...@@ -2,7 +2,7 @@ export(mcmcabn, ...@@ -2,7 +2,7 @@ export(mcmcabn,
query, query,
plot.mcmcabn, plot.mcmcabn,
summary.mcmcabn, summary.mcmcabn,
print.mcmcabn print.summary.mcmcabn
) )
importFrom(gRbase, topoSortMAT) importFrom(gRbase, topoSortMAT)
...@@ -19,4 +19,6 @@ print.mcmcabn ...@@ -19,4 +19,6 @@ print.mcmcabn
#methods for class mcmcabn #methods for class mcmcabn
S3method(plot, mcmcabn) S3method(plot, mcmcabn)
S3method(summary, mcmcabn) S3method(summary, mcmcabn)
S3method(print, mcmcabn)
#method for class summary.mcmcabn
S3method(print, summary.mcmcabn)
...@@ -9,3 +9,4 @@ ...@@ -9,3 +9,4 @@
## mcmcabn 0.3: ## mcmcabn 0.3:
* update mcmcabn to make it compatible with constraints imported from the cache of scores. Heating parameter to increase or decrease acceptance probability. * update mcmcabn to make it compatible with constraints imported from the cache of scores. Heating parameter to increase or decrease acceptance probability.
* new article
############################################################################### summary.mcmcabn.R --- Author : Gilles Kratzer Last modified : 19.02.2019 :
print.summary.mcmcabn <- function(x, ...) {
cat("MCMC summary:\n", sep = "")
cat("Number of burn-in steps: ", (x$object$burnin), "\n", sep = "")
cat("Number of MCMC steps: ", (x$object$iterations), "\n", sep = "")
cat("Thinning: ", (x$object$thinning), "\n\n", sep = "")
cat("Maximum score: ", max(x$object$scores), "\n", sep = "")
cat("Empirical mean: ", mean(x$object$scores), "\n", sep = "")
cat("Empirical standard deviation: ", format(sd(x$object$scores),...), "\n", sep = "")
cat("Quantiles of the posterior network score:\n")
print(x$quant, ...)
cat("\n\nGlobal acceptance rate: ", format( 1 - mean(x$object$rejection), ...), "\n", sep = "")
print(x$AR, ...)
cat("\n\nSample size adjusted for autocorrelation: ",
format( unname(effectiveSize(x$object$scores)), ...), "\n", sep = "")
cat("\nAutocorrelations by lag:\n")
print(x$acf, ...)
} #EOF
...@@ -2,41 +2,17 @@ ...@@ -2,41 +2,17 @@
summary.mcmcabn <- function(object, quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975), lag.max = 10, ...) { summary.mcmcabn <- function(object, quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975), lag.max = 10, ...) {
cat("MCMC summary:\n", sep = "")
cat("Number of burn-in steps: ", (object$burnin), "\n", sep = "")
cat("Number of MCMC steps: ", (object$iterations), "\n", sep = "")
cat("Thinning: ", (object$thinning), "\n\n", sep = "")
cat("Maximum score: ", max(object$scores), "\n", sep = "")
cat("Empirical mean: ", mean(object$scores), "\n", sep = "")
cat("Empirical standard deviation: ", format(sd(object$scores),...), "\n", sep = "")
cat("Quantiles of the posterior network score:\n")
out1 <- matrix(data = quantile(x = object$scores, probs = quantiles), nrow = 1, ncol = length(quantiles), dimnames = list("BN score", out1 <- matrix(data = quantile(x = object$scores, probs = quantiles), nrow = 1, ncol = length(quantiles), dimnames = list("BN score",
quantiles)) quantiles))
print(out1, ...)
cat("\n\nGlobal acceptance rate: ", format( 1 - mean(object$rejection), ...), "\n", sep = "")
out2 <- matrix(data = table(object$method, object$rejection), ncol = 2, dimnames = list(levels(factor(object$method)), out2 <- matrix(data = table(object$method, object$rejection), ncol = 2, dimnames = list(levels(factor(object$method)),
c("Accepted", "Rejected"))) c("Accepted", "Rejected")))
print(out2, ...) out3 <- matrix(data = acf(object$scores, lag.max = 10, plot = FALSE)$acf, nrow = 1, ncol = (lag.max + 1), dimnames = list("acf",
cat("\n\nSample size adjusted for autocorrelation: ",
format( unname(effectiveSize(object$scores)), ...), "\n", sep = "")
# unname(acf(object$scores, lag.max = 10, plot = FALSE))
cat("\nAutocorrelations by lag:\n")
out2 <- matrix(data = acf(object$scores, lag.max = 10, plot = FALSE)$acf, nrow = 1, ncol = (lag.max + 1), dimnames = list("acf",
acf(object$scores, lag.max = 10, plot = FALSE)$lag)) acf(object$scores, lag.max = 10, plot = FALSE)$lag))
print(out2, ...) ans <- list( object=object, quant=out1, AR=out2, acf=out3)
class( ans) <- "summary.mcmcabn"
return( ans)
} #EOF } #EOF
...@@ -29,16 +29,6 @@ colnames(asia) <- c("Asia", ...@@ -29,16 +29,6 @@ colnames(asia) <- c("Asia",
"XRay", "XRay",
"Dyspnea") "Dyspnea")
## Defining the distribution list
dist.asia <- list(Asia = "binomial",
Smoking = "binomial",
Tuberculosis = "binomial",
LungCancer = "binomial",
Bronchitis = "binomial",
Either = "binomial",
XRay = "binomial",
Dyspnea = "binomial")
bsc.compute.asia <- buildscorecache(data.df = asia, bsc.compute.asia <- buildscorecache(data.df = asia,
data.dists = dist.asia, data.dists = dist.asia,
max.parents = 2) max.parents = 2)
......
...@@ -20,30 +20,6 @@ ...@@ -20,30 +20,6 @@
library(bnlearn) #for the dataset library(bnlearn) #for the dataset
library(abn) #for the cache of scores computing function library(abn) #for the cache of scores computing function
#renaming columns of the dataset
colnames(asia) <- c("Asia",
"Smoking",
"Tuberculosis",
"LungCancer",
"Bronchitis",
"Either",
"XRay",
"Dyspnea")
#lets define the distribution list
dist.asia <- list(Asia = "binomial",
Smoking = "binomial",
Tuberculosis = "binomial",
LungCancer = "binomial",
Bronchitis = "binomial",
Either = "binomial",
XRay = "binomial",
Dyspnea = "binomial")
bsc.compute.asia <- buildscorecache(data.df = asia,
data.dists = dist.asia,
max.parents = 2)
mcmc.out.asia <- mcmcabn(score.cache = bsc.compute.asia, mcmc.out.asia <- mcmcabn(score.cache = bsc.compute.asia,
score = "mlik", score = "mlik",
data.dists = dist.asia, data.dists = dist.asia,
......
...@@ -197,8 +197,8 @@ mcmc.out.asia.small <- mcmcabn(score.cache = bsc.compute.asia, ...@@ -197,8 +197,8 @@ mcmc.out.asia.small <- mcmcabn(score.cache = bsc.compute.asia,
score = "mlik", score = "mlik",
data.dists = dist.asia, data.dists = dist.asia,
max.parents = 2, max.parents = 2,
mcmc.scheme = c(100,0,0), mcmc.scheme = c(150,0,0),
seed = 321, seed = 41242,
verbose = FALSE, verbose = FALSE,
start.dag = "random", start.dag = "random",
prob.rev = 0.03, prob.rev = 0.03,
......
% print-summar.mcmcabn.Rd ---
% Author : Gilles Kratzer
% Created on : 31.10.2019
% Last modification :
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\name{print.summary.mcmcabn}
\alias{print.summary.mcmcabn}
\title{Methods for printing the summary of mcmcabn objects}
\description{
Method for printing the summary of \code{mcmcabn} objects.
}
\usage{
\method{print}{summary.mcmcabn}(x, \dots)
}
\arguments{
\item{x}{an object of class \code{summary.mcmcabn}.}
\item{\dots}{additional arguments passed to \code{print}.}
}
\details{
There exists a \code{\link{summary}} S3 function that displays more details.
}
\author{Gilles Kratzer}
\examples{
## Example from the asia dataset from Lauritzen and Spiegelhalter (1988)
## provided by Scutari (2010)
summary(mcmc.out.asia)
}
No preview for this file type
...@@ -172,6 +172,8 @@ tab <- apply(X = u.list.dag, MARGIN = 3, FUN = function(x){ ...@@ -172,6 +172,8 @@ tab <- apply(X = u.list.dag, MARGIN = 3, FUN = function(x){
nb.dag <- 10 nb.dag <- 10
scores.dags <- vector(length = nb.dag) scores.dags <- vector(length = nb.dag)
num.arcs <- vector(length = nb.dag) num.arcs <- vector(length = nb.dag)
shd <- vector(length = nb.dag)
for(i in 1:nb.dag){ for(i in 1:nb.dag){
dag <- u.list.dag[,,order(tab,decreasing = TRUE)[i]] dag <- u.list.dag[,,order(tab,decreasing = TRUE)[i]]
...@@ -179,6 +181,7 @@ for(i in 1:nb.dag){ ...@@ -179,6 +181,7 @@ for(i in 1:nb.dag){
fabn <- fitabn(dag.m = (dag),data.df = asia,data.dists = dist.asia,method = "bayes") fabn <- fitabn(dag.m = (dag),data.df = asia,data.dists = dist.asia,method = "bayes")
scores.dags[i] <- fabn$mlik scores.dags[i] <- fabn$mlik
num.arcs[i] <- sum(dag) num.arcs[i] <- sum(dag)
shd[i] <- compareDag(ref = u.list.dag[,,order(tab,decreasing = TRUE)[1]],u.list.dag[,,order(tab,decreasing = TRUE)[i]])$`Hamming-distance`
} }
``` ```
...@@ -198,4 +201,6 @@ plot(x = 1:nb.dag,y = scores.dags,col="red", type = 'b', lwd=2, axes = FALSE, xl ...@@ -198,4 +201,6 @@ plot(x = 1:nb.dag,y = scores.dags,col="red", type = 'b', lwd=2, axes = FALSE, xl
axis(4, col.axis = 'red') axis(4, col.axis = 'red')
mtext("DAGs scores", side=4, line=1.5, col="red") mtext("DAGs scores", side=4, line=1.5, col="red")
axis(1, col.axis = 'black',at = 1:nb.dag,labels = num.arcs) axis(1, col.axis = 'black',at = 1:nb.dag,labels = num.arcs)
axis(3, col.axis = 'orange',at = 1:nb.dag,labels = shd)
mtext("Structural Hamming distances", side=3, line=2, col="orange")
``` ```
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