# R/f02.chain.code.r In rbmn: Handling Linear Gaussian Bayesian Networks

#### Documented in chain2correlationchain2gemachain2mnchain2nbnchain2prechain4chaincheck8chaingenerate8chaininout4chainis8nbn8chainmarginal4chainnbn2chainorder4chainprint8chainreverse8chainstate4chain

```#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
generate8chain <- function(rnn=c(3,7),proo=0.5,rcor=c(-1,1),
rmu=c(0,0),rsig=c(0,1),
nona=r.form3names(max(rnn)))
#TITLE generation of a /chain/ /nbn/
#DESCRIPTION [randomly] generates a /chain/ /nbn/.
#DETAILS
# Proposed ranges can be a unique value, implying no randomness
# in the value.\cr
# Roots are placed according to \code{proo} probabilities, then collider
# are placed in between with uniform probability on the possibles nodes.
#KEYWORDS
#INPUTS
#[INPUTS]
#{rnn} <<Range of the number of nodes.>>
#{proo} <<Probabilit[y|ies] that the successive and acceptable nodes be colliders. Can be a vector.>>
#{rcor} <<Range of the correlations between neighbour nodes.>>
#{rmu}  <<Range of the expectations.>>
#{rsig} <<Range of the standard deviations.>>
#{nona} <<Proposed names for the maximum number of nodes, only the
#         necessary first ones will be used.>>
#VALUE
# A /chain/ coding list is returned.
#EXAMPLE
# set.seed(1234);
# print8chain(generate8chain());
# print8chain(generate8chain());
# print8chain(generate8chain(rnn=10,rcor=0.5));
# print8chain(generate8chain(rnn=10,rcor=0.5));
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 13_07_03
#REVISED 13_07_03
#--------------------------------------------
{
# checking
# < to be done >
# dealing with unique values
if (length(rnn)==1) { rnn <- c(rnn,rnn);}
if (length(rcor)==1) { rcor <- c(rcor,rcor);}
if (length(rmu)==1) { rmu <- c(rmu,rmu);}
if (length(rsig)==1) { rsig <- c(rsig,rsig);}
# getting the node number.
nn <- floor(runif(1,rnn[1],rnn[2]+1))
# getting the node names
nona <- nona[r.bc(nn)]
# getting the probabilities of downstream
pp <- proo
while (length(pp) < nn) { pp <- c(pp,proo);}
pp <- pp[1:nn]
# building the /chain/
# miscellaneous part
res <- vector("list",0)
res\$names <- nona
res\$mu <- runif(nn,rmu[1],rmu[2])
res\$sigma <- runif(nn,rsig[1],rsig[2])
res\$corre <- runif(nn-1,rcor[1],rcor[2])
# placing roots first
root <- rbinom(nn,1,pp)
for (ii in r.bc(nn-1)) {
if (root[ii]==1) { root[ii+1] <- 0;}
}
# at least one root
if (sum(root)==0) {
root[sample.int(nn,1)] <- 1
}
res\$roots <- nona[root==1]
root <- which(root==1)
# placing colliders second
res\$colliders <- character(0)
for (ii in r.bc(length(root)-1)) {
r1 <- root[ii]; r2 <- root[ii+1]
if (r2-r1<2) { stop("Erreur in generate8chain");}
qui <- sample.int(r2-r1-1,1)+r1
res\$colliders <- c(res\$colliders,nona[qui])
}
# returning
res
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
print8chain <- function(chain,digits=3)
#TITLE prints a /chain/ object
#DESCRIPTION prints a /chain/ object.
#DETAILS
# See \code{nbn2chain} code for some details about the
# definition of a /chain/.
#KEYWORDS
#INPUTS
#{chain} << The \code{chain} object to print.>>
#[INPUTS]
#{digits} << when not null, the number of digits for rounding the
# numerical values.>>
#VALUE
# nothing but something is printed
#EXAMPLE
# print8chain(rbmn0chain.01);
# print8chain(rbmn0chain.02);
# print8chain(rbmn0chain.03);
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 12_01_18
#REVISED 12_01_30
#--------------------------------------------
{
# checking
che <- check8chain(chain)
if (is.character(che)) {
print(che)
stop("The provided 'chain' is not valid!")
}
# node number
nn <- length(chain\$names)
orien <- aux2(chain)/2+1.5
cat("#-------------------------\n")
for (node in r.bc(nn)) {
nna <- chain\$names[node]
nona <- paste("(",nna,")",sep="")
if (nna %in% chain\$roots) {
nona <- paste("<--(",nna,")-->",sep="")
}
if (nna %in% chain\$colliders) {
nona <- paste("-->)",nna,"(<--",sep="")
}
cat(r.form3justifie(nona,nbc=10,format=2))
cat("  ",round(chain\$mu[node],digits),"              (",
round(chain\$sigma[node],digits),")",sep="")
cat("\n")
if (node < nn) {
cat(r.form3justifie(c("^","v")[orien[node]],nbc=10,format=2))
cat("     <",round(chain\$corre[node],digits),">",sep="")
cat("\n")
}
}
cat("#-------------------------\n")
# returning
invisible()
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
check8chain <- function(chain)
#TITLE checks a /chain/ object
#DESCRIPTION checks the consistency of \code{chain} as a /chain/ object
# issues a fatal error with some clues if inconsistent.
#DETAILS
# Looking a the code of this function provides a way to know which
# are the requirements of a /chain/ object.
#KEYWORDS
#INPUTS
#{chain} << The \code{chain} object to check.>>
#[INPUTS]
#VALUE
# \code{TRUE} or a \code{character} containing some clue
#EXAMPLE
# check8chain(rbmn0chain.01);
# if (is.na(as.logical(res))) { print(res);}
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 13_07_26
#REVISED 13_07_26
#--------------------------------------------
{
# checking the component names
compo <- names(rbmn0chain.01)
if (!setequal(compo,names(chain))) {
return("The names of the chain are not 'names(rbmn0chain.01)'")
}
# checking the lengths
nbno <- length(chain\$names)
res <- character(0)
if (nbno != length(chain\$mu)) { res <- c(res,"'\$mu' has got a bad length");}
if (nbno != length(chain\$mu)) { res <- c(res,"'\$sigma' has got a bad length");}
if ((nbno-1) != length(chain\$corre)) { res <- c(res,"'\$corre' has got a bad length");}
if ((nbno-1)/2 < length(chain\$colliders)) { res <- c(res,"'\$colliders' has got a bad length");}
if ((nbno+1)/2 < length(chain\$roots)) { res <- c(res,"'\$roots' has got a bad length");}
if (length(res) > 0) { return(res);}
# checking the correlations
if (any(abs(chain\$corre)>1)) { return("Some correlation doesn't make sense");}
# returning
TRUE
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
order4chain <- function(chain,ord=NULL)
#TITLE returns a topological order of a /chain/ or
# checks a proposed order.
#DESCRIPTION From a \code{chain} object
# returns one of the possible topological orders,
# through a permutation when \code{is.null(ord)}.
# If not \code{ord} must be a proposed order to be
# checked given as a permutation if \code{is.numeric(ord)}
# or a vector of ordered names if \code{is.character(ord)}.
#DETAILS
# For the moment the \code{ord} option is
# bad and an error message is returned when used.
#KEYWORDS
#INPUTS
#{chain}<< the \code{chain} object to be considered.>>
#[INPUTS]
#{ord} << Indicates what must be done. \code{NULL} to get a topological
# order associated to the chain otherwise a permutation to be checked as
# one of the possible topological orders of the chain.>>
#VALUE
# a permutation vector of the nodes of the /nbn/
#        or a named character with the nodes not having
#        their parents before them; when it is of
#        length zero this means that the check
#        was successful.
#EXAMPLE
# order4chain(rbmn0chain.02);
# order4chain(rbmn0chain.02,order4chain(rbmn0chain.02));
#REFERENCE
#CALLING
#COMMENT
#FUTURE Correct the 'ord' option!
#AUTHOR J.-B. Denis
#CREATED 12_01_31
#REVISED 12_01_31
#--------------------------------------------
{
# checking
# < to be done >
# getting the permutation
ooo <- aux2(chain,TRUE)
if (is.null(ord)) {
res <- ooo
} else {
return("Sorry, the 'ord' option is wrong and must be corrected!")
nam <- chain\$names
nn <- length(nam)
if (is.numeric(ord)) {
# a numeric permutation is expected
if (length(union(ord,r.bc(nn)))!=nn) {
r.erreur(ord,message="'ord' is not a permutation.")
}
} else {
# a list of names is expected
if (length(union(ord,nam))!=nn) {
r.erreur(list(ord,nam),
message="'ord' is not a permutation of 'nam'")
}
ord <- r.numero(ord,nam)
}
# checking and storing inconsistencies
res <- numeric(0); nm <- character(0)
names(res) <- nm
etat <- state4chain(chain)
for (uu in ord) {
if (etat[uu] == "r") {
if (uu > 1) {
if (etat[uu-1] == "t") { etat[uu-1] <- "r";}
if (etat[uu-1] == "c") { etat[uu-1] <- "t";}
}
if (uu < nn) {
if (etat[uu+1] == "t") { etat[uu+1] <- "r";}
if (etat[uu+1] == "c") { etat[uu+1] <- "t";}
}
} else {
res <- c(res,uu)
nm <- c(nm,nam[uu])
names(res) <- nm
}
etat[uu] <- "f"
}
}
# returning
res
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
state4chain <- function(chain)
#TITLE returns the states of each node of a chain
#DESCRIPTION From a \code{chain} object
# returns a named character precising the role of each node:
# "r" for root, "c" for collider, "t" for transmitter and
# "l" for leaf.
#DETAILS
#KEYWORDS
#INPUTS
#{chain}<< the \code{chain} object to be considered.>>
#[INPUTS]
#VALUE
# a character of the states named with node names.
#EXAMPLE
# state4chain(rbmn0chain.01);
# state4chain(rbmn0chain.03);
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 12_01_31
#REVISED 12_01_31
#--------------------------------------------
{
# checking
# < to be done >
# identifying
nn <- length(chain\$names)
etat <- rep("t",nn)
etat[c(1,nn)] <- "l"
if (length(chain\$colliders) > 0) {
etat[r.numero(chain\$colliders,chain\$names)] <- "c"
}
etat[r.numero(chain\$roots,chain\$names)] <- "r"
# naming
names(etat) <- chain\$names
# returning
etat
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
marginal4chain <- function(chain)
#TITLE returns marginal expectations and standard deviations of a chain
#DESCRIPTION From a \code{chain} object
# returns a list with two components: \code{\$mu} and \code{\$sigma}
# vectors of marginal expectations and standard deviations.\cr
#DETAILS
#KEYWORDS
#INPUTS
#{chain}<< the \code{chain} object to be considered.>>
#[INPUTS]
#VALUE
# a list with the two components \code{\$mu} and \code{\$sigma}.
#EXAMPLE
# marginal4chain(rbmn0chain.02)
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 12_01_31
#REVISED 12_01_31
#--------------------------------------------
{
# checking
# < to be done >
# getting a topological order
# and the status of each node
nn <- length(chain\$names)
ooo <- order4chain(chain)
etat <- state4chain(chain)
# following this order computing
# progressively expectations and
# standard deviations
mu <- sigma <- rep(NA,nn)
orien <- aux2(chain)
for (uu in ooo) {
if (etat[uu]=="r") {
# root, no modification
mu[uu] <- chain\$mu[uu]
sigma[uu] <- chain\$sigma[uu]
} else {
if (etat[uu]=="c") {
# collider: the two neighbours are parents
sigma[uu] <- chain\$sigma[uu]/aux0(chain\$corre[(uu-1):uu]);
mu[uu] <- chain\$mu[uu] +
chain\$corre[uu-1]*sigma[uu]/sigma[uu-1]*mu[uu-1] +
chain\$corre[uu]  *sigma[uu]/sigma[uu+1]*mu[uu+1]
} else {
# other cases: only one parent
if (uu>1) { if (orien[uu-1]==1) {
# parent upward
sigma[uu] <- chain\$sigma[uu]/aux0(chain\$corre[uu-1]);
mu[uu] <- chain\$mu[uu] +
chain\$corre[uu-1]*sigma[uu]/sigma[uu-1]*mu[uu-1]

}}
if (uu<nn) { if (orien[uu]==-1) {
# parent downward
sigma[uu] <- chain\$sigma[uu]/aux0(chain\$corre[uu]);
mu[uu] <- chain\$mu[uu] +
chain\$corre[uu]*sigma[uu]/sigma[uu+1]*mu[uu+1]
}}
}
}
}
if (any(is.na(mu))) {
r.erreur(chain,message="either the chain or the algo")
}
# returning
list(mu=mu,sigma=sigma)
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
chain2nbn <- function(chain)
#TITLE transforms a /chain/ to a /nbn/
#DESCRIPTION From a \code{chain} object
# returns the \code{nbn} translation.
#DETAILS
#KEYWORDS
#INPUTS
#{chain}<< the \code{chain} object to be transformed.>>
#[INPUTS]
#VALUE
# The corresponding \code{nbn} object.
#EXAMPLE
# print8nbn(chain2nbn(rbmn0chain.02),ordering=names(rbmn0nbn.02))
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 12_01_19
#REVISED 13_07_04
#--------------------------------------------
{
# checking
# < to be done >
# constants
nna <- chain\$names
nn <- length(nna)
ooo <- aux2(chain,TRUE); # topological order
rho <- aux2(chain); # arc orientations
avan <- c(0,rho)== 1; # parent before
apre <- c(rho,0)==-1; # parent after
aaa <- sort(r.numero(chain\$roots,chain\$names))
# transforming
res <- vector("list",nn)
sigma <- marginal4chain(chain)\$sigma
mu <- chain\$mu
for (nnu in ooo) {
papa <- character(0); # parent specification
coef <- numeric(0);   # regression coefficients
# only non roots must be filled
if (!(nnu %in% aaa)) {
corri <- 1
if (avan[nnu] & apre[nnu]) {
corri <- aux0(chain\$corre[(nnu-1):nnu])
} else {
if (avan[nnu]) { corri <- aux0(chain\$corre[nnu-1]);}
if (apre[nnu]) { corri <- aux0(chain\$corre[nnu  ]);}
}
if (avan[nnu]) {
# the node before is a parent
papa <- nna[nnu-1]
coef <- chain\$sigma[nnu] / sigma[nnu-1] /
corri * chain\$corre[nnu-1]
}
if (apre[nnu]) {
# the node after is a parent
papa <- c(papa,nna[nnu+1])
coef <- c(coef,chain\$sigma[nnu] / sigma[nnu+1] /
corri * chain\$corre[nnu])
}
}
res[[nnu]]\$parents <- papa
res[[nnu]]\$regcoef <- coef
res[[nnu]]\$mu <- chain\$mu[nnu]
res[[nnu]]\$sigma <- chain\$sigma[nnu]
}
names(res) <- nna
# reordering the nodes
ooo <- order4chain(chain)
res <- res[ooo]
# testing possible NaNs
for (ii in r.bf(res)) {
if (any(is.na(res[[ii]]\$regcoef))) {
stop("Dectected Numerical Inaccuracy in 'chain2nbn'")
}
}
# returning
res
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
chain2gema <- function(chain)
#TITLE transforms a /chain/ to a /gema/
#DESCRIPTION From a \code{chain} object
# returns the \code{gema} using a direct formulae.\cr
# Much precised than to use the /nbn/ way.
#DETAILS
#KEYWORDS
#INPUTS
#{chain}<< the \code{chain} object to be transformed.>>
#[INPUTS]
#VALUE
# The corresponding \code{gema} object.
#EXAMPLE
# identical(chain2gema(rbmn0chain.02)\$mu,rbmn0gema.02\$mu);
# print(chain2gema(rbmn0chain.02)\$li-rbmn0gema.02\$li);
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 12_02_06
#REVISED 12_02_06
#--------------------------------------------
{
# checking
# < to be done >
# useful constants
nna <- chain\$names
nn <- length(nna)
pare <- aux2(chain,parents=TRUE)
# initializing
res <- list(mu=marginal4chain(chain)\$mu,
li=diag(nrow=nn))
names(res\$mu) <- nna
dimnames(res\$li) <- list(nna,NULL)
# computing linear combination
## for correlations colliding structure
for (no in r.bc(nn)) {
if (pare[no] == -1) {
# a colliding node
res\$li[no,no-1] <- chain\$corre[no]
res\$li[no,no] <- aux0(chain\$corre[c(no,no+1)])
res\$li[no,no+1] <- chain\$corre[no+1]
} else {
if (pare[no] > 0) {
# a standard node
if (pare[no] < no) { rho <- chain\$corre[no-1]
} else { rho <- chain\$corre[no];}
res\$li[no,pare[no]] <- rho
res\$li[no,no] <- aux0(rho)
}
}
}
## for correlations sequential segments
if (length(chain\$colliders)>0) {
ccc <- r.numero(chain\$colliders,nna)
} else {
ccc <- numeric(0)
}
aaa <- r.numero(chain\$roots,nna)
ava <- aaa[-length(aaa)]
apr <- aaa[-1]
for (cc in r.bf(ccc)) {
pivo <- res\$li[ccc[cc],ccc[cc]]
ba <- ava[cc]:ccc[cc]
bp <- ccc[cc]:apr[cc]
res\$li[ba,ba] <- aux1(chain\$corre[ba[-length(ba)]],FALSE)
res\$li[bp,bp] <- aux1(chain\$corre[bp[-length(bp)]],TRUE)
res\$li[ccc[cc],ccc[cc]] <- pivo
}
## without forgetting endind sequential segments
if (aaa[1] > 1) {
aa <- 1:aaa[1]
res\$li[aa,aa] <- aux1(chain\$corre[1:(aaa[1]-1)],TRUE)
}
if (aaa[length(aaa)] < nn) {
aa <- aaa[length(aaa)]:nn
res\$li[aa,aa] <- aux1(chain\$corre[aaa[length(aaa)]:(nn-1)],FALSE)
}
#
masig <- marginal4chain(chain)\$sigma
res\$li <- res\$li * masig
# reordering the nodes
ooo <- order4chain(chain)
res\$mu <- res\$mu[ooo]
res\$li <- res\$li[ooo,ooo]
# returning
res
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
chain2correlation <- function(chain)
#TITLE computes the correlation matrix of a chain
#DESCRIPTION returns the correlation
# matrix of a /chain/ object.
#DETAILS
#KEYWORDS
#INPUTS
#{chain}<< The chain object to consider.>>
#[INPUTS]
#VALUE
# The correlation matrix. It is not
# sorted to respect a topological order
# contrary to \code{chain2mn} function.
#EXAMPLE
# chain2correlation(rbmn0chain.03);
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_02_01
#REVISED 11_02_02
#--------------------------------------------
{
# checking
# <to be done>
nn <- length(chain\$mu)
res <- matrix(0,nn,nn)
dimnames(res) <- list(chain\$names,chain\$names)
# finding the milestones
etat <- state4chain(chain)
sto <- sort(c(1,nn,which("c"==etat)))
# introducing the blocks
for (bb in r.bc(length(sto)-1)) {
deb <- sto[bb]; fin <- sto[bb+1]
res[r.bd(deb,fin),r.bd(deb,fin)] <- aux5(chain\$corre[r.bd(deb,fin-1)])
}
# returning
res
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
chain2mn <- function(chain,order=TRUE)
#TITLE computes the distribution of a chain
#DESCRIPTION returns the /mn/ object
# associated to a /chain/ object. Much better
# to use this function that the general function
# \code{nbn2mn} since exact formulae are applied.
#DETAILS
#KEYWORDS
#INPUTS
#{chain}<< The chain object to consider.>>
#[INPUTS]
#{order} << Must a topological order be imposed?>>
#VALUE
# The resulting /mn/ object. Following the
# convention of \code{mn} objects, a topological
# order is given to it. This is necessary to retrieve
# the associate /nbn/.
#EXAMPLE
# print8mn(chain2mn(rbmn0chain.01));
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_12_05
#REVISED 13_07_30
#--------------------------------------------
{
# checking
# <to be done>
# getting the marginal parameters
margi <- marginal4chain(chain)
# getting the multivariable parameters
corre <- chain2correlation(chain)
# computing the variance matrix
gamma <- outer(margi\$sigma,margi\$sigma,"*") * corre
# naming the components
names(margi\$mu) <- chain\$names
dimnames(gamma) <- list(chain\$names,chain\$names)
# getting a topological order
if (order) {
ooo <- order4chain(chain)
} else {
ooo <- r.bf(margi\$mu)
}
# returning
list(mu=margi\$mu[ooo],
gamma=gamma[ooo,ooo])
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
chain2pre <- function(chain,corre=FALSE)
#TITLE computes the precision of a chain
#DESCRIPTION returns the precision matrix
# of a chain, that is the inverse of its
# variance (correlation) matrix. Much better
# to use this function that
# \code{solve(chain2mn(chain)\$gamma)} since
# exact formulae are applied.
#DETAILS
#KEYWORDS
#INPUTS
#{chain}<< The chain object to consider.>>
#[INPUTS]
#{corre} <<To get the inverse of the correlation
#VALUE
# A dimnamed matrix
#EXAMPLE
# chain2pre(rbmn0chain.02);
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_02_02
#REVISED 11_02_06
#--------------------------------------------
{
# checking
# <to be done>
# building the inverse of the correlation matrix
# constants
nn <- length(chain\$names)
## first a long sequential chain
res <- aux6(chain\$corre)
## then the 3x3 blocks for colliders
coli <- sort(which("c"==state4chain(chain)))
for (cc in coli) {
ends <- c(cc==2,cc==nn-1)
ou <- (cc-1):(cc+1)
res[ou,ou] <- aux7(chain\$corre[ou[-3]],ends)
}
if (!corre) {
sisi <- marginal4chain(chain)\$sigma
res <- res / outer(sisi,sisi,"*")
}
# dimnaming
dimnames(res) <- list(chain\$names,chain\$names)
# returning
res
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
inout4chain <- function(chain)
#TITLE reduces a chain to its inputs and outputs
#DESCRIPTION From a \code{chain} returns the
# reduced \code{chain} comprising only inputs
# (that is root nodes) and outputs (that is
# colliders and ends which are not roots)
#DETAILS
#KEYWORDS
#INPUTS
#{chain}<< The chain object to consider.>>
#[INPUTS]
#VALUE
# The resulting chain
#EXAMPLE
# print8chain(inout4chain(rbmn0chain.02));
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_02_07
#REVISED 11_02_07
#--------------------------------------------
{
# checking
# <to be done>
# determining the node to keep
eta <- state4chain(chain)
kee <- chain\$names[which(eta!="t")]
nk <- length(kee)
# building the resulting chain
res <- vector("list",0)
res\$names <- kee
res\$roots <- chain\$roots
res\$colliders <- chain\$colliders
if (nk < 2) {
res\$mu <- chain\$mu
res\$sigma <- chain\$sigma
res\$corre <- chain\$corre
} else {
mn <- chain2mn(chain)
mn\$mu <- mn\$mu[kee]
mn\$gamma <- mn\$gamma[kee,kee,drop=FALSE]
res\$mu <- mn\$mu
res\$sigma <- sqrt(diag(mn\$gamma))
co <- cor4var(mn\$gamma)
res\$corre <- diag(co[-1,-nk,drop=FALSE])
res <- aux4(res)
}
# returning
res
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
reverse8chain <- function(chain)
#TITLE reverses the nodes of a chain
#DESCRIPTION returns the chain obtained
# after reversing its node order
#DETAILS
#KEYWORDS
#INPUTS
#{chain}<< The chain object to consider.>>
#[INPUTS]
#VALUE
# The resulting chain
#EXAMPLE
# print8chain(rbmn0chain.02);
# print8chain(reverse8chain(rbmn0chain.02));
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_02_07
#REVISED 11_02_07
#--------------------------------------------
{
# checking
# <to be done>
# building it
res <- vector("list",0)
res\$names <- rev(chain\$names)
res\$roots <- chain\$roots
res\$colliders <- chain\$colliders
res\$mu <- rev(chain\$mu)
res\$sigma <- rev(chain\$sigma)
res\$corre <- rev(chain\$corre)
# returning
res
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
is8nbn8chain <- function(nbn,order=FALSE)
#TITLE Checks if a given /nbn/ is a /chain/
#DESCRIPTION returns \code{TRUE} [the order] or
# \code{FALSE} [NULL] according that \code{nbn}
# is a chain of not [according to \code{order}].
#DETAILS
#KEYWORDS
#INPUTS
#{nbn}<< The nbn object to consider.>>
#[INPUTS]
#{order} << When \code{FALSE} the answer to the
# question is returned with \code{TRUE} or \code{FALSE}.\cr
# When \code{TRUE} the chain order
# of the nodes is returned if it is a /chain/
# else \code{NULL}.>>
#VALUE
# A \code{logical(1)} when \code{order} si \code{TRUE} if not
#  the resulting chain order versus NULL.
#EXAMPLE
#is8nbn8chain(rbmn0nbn.01);
#is8nbn8chain(rbmn0nbn.04);
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_02_14
#REVISED 11_02_14
#--------------------------------------------
{
# checking
# <to be done>
# constants
nn <- length(names(nbn))
# degenerate case
if (nn == 1) {
if (order) {
return(1)
} else {
return(TRUE)
}
}
# symmetrical relationship matrix
rela <- rela + t(rela)
# checking obvious cases
res <- TRUE
nbr <- apply(rela,1,sum)
if (sum(nbr == 1)!= 2) { res <- FALSE;}
if (sum(nbr == 2)!= nn-2) { res <- FALSE;}
if (!res) {
if (order) {
return(NULL)
} else {
return(FALSE)
}
}
# looking for a chain between nodes
anc <- which(nbr==1)[1]
orde <- anc
mauvais <- FALSE
while ((length(orde)<nn) & !mauvais) {
nou <- which(rela[anc,]==1)
if (length(nou) != 1) {
mauvais <- TRUE
} else {
rela[anc,nou] <- rela[nou,anc] <- 0
orde <- c(orde,nou)
anc <- nou
}
}
#
if (mauvais) {
if (order) {
res <- NULL
} else {
res <- FALSE
}
} else {
if (!order) {
res <- TRUE
} else {
res <- orde
}
}
# returning
res
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
nbn2chain <- function(nbn)
#TITLE transforms a /nbn/ into a /chain/
#DESCRIPTION returns the chain obtained
# from \code{nbn} which is supposed to a chain.
# If it is not a chain, an error is issued.
#DETAILS
# It is advised to use \code{is8nbn8chain} before
# calling this function.
#KEYWORDS
#INPUTS
#{nbn}<< The /nbn/ object to consider.>>
#[INPUTS]
#VALUE
# The resulting chain
#EXAMPLE
# print8chain(nbn2chain(rbmn0nbn.02));
#REFERENCE
#CALLING
#COMMENT
#FUTURE
#AUTHOR J.-B. Denis
#CREATED 11_02_14
#REVISED 11_02_15
#--------------------------------------------
{
# checking
# <to be done for /nbn/>
if (!is8nbn8chain(nbn)) {
r.erreur(nbn,message="This /nbn/ is not a /chain/")
}
# getting one of the two possible chain orders
ooo <- is8nbn8chain(nbn,TRUE)
# building the chain except for correlations
res <- vector("list",0)
res\$names <- names(nbn)[ooo]
res\$sigma <- sapply(nbn,function(a){a\$sigma;})[ooo]
res\$mu    <- sapply(nbn,function(a){a\$mu;})[ooo]
nbpa <- sapply(nbn,function(a){length(a\$parents);})
nbpa <- nbpa[ooo]
res\$roots <- res\$names[which(nbpa==0)]
res\$colliders <- res\$names[which(nbpa==2)]
res\$corre <- rep(0.5,length(res\$names)-1)
# getting a topological order
too <- order4chain(res)
res\$corre <- rep(0.5,length(res\$names)-1)
sig <- rep(NA,length(res\$names))
if (length(res\$names)>1) {
for (nnu in r.bf(res\$names)) {
nnn <- too[nnu]
nna <- res\$names[nnn]
if (!(nna %in% res\$roots)) {
papa <- nbn[[nna]]\$parents
# computing the necessary standard deviations
sb <- res\$sigma[nnn]
coco <- sasa <- numeric(0)
for (pn in r.bf(papa)) {
pp <- papa[pn]
reco <- nbn[[nna]]\$regcoef[pn]
if (nnn > 1) {if (res\$names[nnn-1]==pp) {
sa <- sig[nnn-1]
sasa <- c(sasa,sa)
coco <- c(coco,reco)
}}
if (nnn < length(res\$names)) {if (res\$names[nnn+1]==pp) {
sa <- sig[nnn+1]
sasa <- c(sasa,sa)
coco <- c(coco,reco)
}}
}
# getting the marginal variance
sb <- sqrt(sum(coco^2*sasa^2)+sb^2)
sig[nnn] <- sb
# computing the correlation(s)
for (pn in r.bf(papa)) {
pp <- papa[pn]
reco <- nbn[[nna]]\$regcoef[pn]
if (nnn > 1) {if (res\$names[nnn-1]==pp) {
sa <- sig[nnn-1]
res\$corre[nnn-1] <- reco * sa / sb
}}
if (nnn < length(res\$names)) {if (res\$names[nnn+1]==pp) {
sa <- sig[nnn+1]
res\$corre[nnn] <-  reco * sa / sb
}}
}
} else {
# conditional and marginal variances are equal
sig[nnn] <- res\$sigma[nnn]
}
}
}
# returning
res
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

#<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
chain4chain <- function(chain,nodes,
condi=numeric(0),
value=rep(0,length(condi)))
#TITLE extracts a chain from a chain
#DESCRIPTION returns the chain obtained
# from \code{chain} retaining only nodes indicated
# by \code{nodes} and conditioned with nodes
# indicated in \code{condi}.
#DETAILS
# Integration is done for nodes not belonging to
# the extracted chain nor being in the conditioning
# subset. Then the distribution of the retained nodes
# is left identical to this in the initial chain.
#KEYWORDS
#INPUTS
#{chain}<< The chain object to consider.>>
#[INPUTS]
#{nodes} << \code{numeric} (or \code{character}) vector
#         giving the numbers (or names) of the nodes
#         to be retained in the extracted chain.>>
#{condi} << \code{numeric} (or \code{character}) vector
#         giving the numbers (or names) of the
#         conditioning nodes for the extracted chain.>>
#{value} << Numerical values associated to \code{condi}.>>
#VALUE
# The resulting chain
#EXAMPLE
# chain4chain(rbmn0chain.02,c("a","d"),c("b"),12)
#REFERENCE
#CALLING
#COMMENT
#FUTURE To be rewritten
#AUTHOR J.-B. Denis
#CREATED 11_02_08
#REVISED 11_02_08
#--------------------------------------------
{
return("This function is under construction!")
# checking
# <to be done>
# constants
nn <- length(chain\$names)
# each case in turn
fait <- FALSE
# no condition
res <- chain
fait <- TRUE
}
#
# both conditioning
if (nn < 3) {
res <- NULL
} else {
res <- "'a faire"
fait <- TRUE
}
}
#
if (!fait) {
if (nn < 2) {
res <- NULL
} else {
if (is.null(tail)) { chain <- reverse8chain(chain);}
# condioning has to be done for the last node
# looking for the last root
lro <- max(which(state4chain(chain)=="r"))
ncor <- aux8(chain\$corre[r.bd(lro,nn)])
ncor <- diag(ncor[-1,-(nn-lro),drop=FALSE])
nmu <- chain\$mu[-nn]
nsi <- chain\$sigma[-nn]
res <- list(names=chain\$names[-nn],
roots=chain\$roots,
colliders=chain\$colliders,
mu=nmu,
sigma=nsi,
corre=ncor
)
fait <- TRUE
if (is.null(tail)) { res <- reverse8chain(res);}
}
}
# returning
res
}
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
```

## Try the rbmn package in your browser

Any scripts or data that you put into this service are public.

rbmn documentation built on July 9, 2023, 6:37 p.m.