# NB the following codes approach the computations in a very direct manner, 
# with all markers brought into a single Bayesian network, and not exploiting 
# the ideas in section 3.4.2. 

# 8 markers simple paternity analysis, for data in Table 1.

markers<-c('D13','D3','D5','D7','FGA','TH01','TPOX','VWA') 

# Evidence: trace=sgt=..
evidence<-list(
	D13='9-14',
	D3='11-17',
	D5='9-11',
	D7='10-10',
	FGA='21-22',
	TH01='7-7',
	TPOX='10-11',
	VWA='18-18'
	)

alleles <-list(
	D13 = c("9", "14", "x"), 
	D3 = c("11", "17", "x" ), 
	D5 = c("9", "11", "x"), 
	D7 = c("10", "x"), 
	FGA = c("21", "22", "x"), 
	TH01 = c("7", "x"), 
	TPOX = c("10", "11", "x"), 
	VWA = c("18", "x")
)

freqsC <-list(
	D13 = c(0.075, 0.048, 0.877), 
	D3 = c(0.002, 0.215, 0.783), 
	D5 = c(0.05, 0.361, 0.589), 
	D7 = c(0.243, 0.757), 
	FGA = c(0.185, 0.219, 0.596), 
	TH01 = c(0.19, 0.81), 
	TPOX = c(0.056, 0.243, 0.701 ), 
	VWA = c(0.2, 0.8)
)

IBD<-function(patt,m,spg,smg,aspg,asmg) 
{ 
# pattern|relation distribution 
zibd<-matrix(0,7,5) 
zibd[1,1]<-1
zibd[c(2,5),2]<-zibd[c(2,4),3]<-zibd[c(1,3),4]<-zibd[c(1,2),5]<-0.5
#
qr<-c('none','1pa2','2pa1','hsma','hspa')%in%vs.relation
zibd<-zibd[,qr,drop=FALSE] 
qp<-apply(zibd,1,sum)>0
zibd<-zibd[qp,,drop=FALSE] 
# 
pool1<-cs('pool1',m);
pool2<-cs('pool2',m); pool3<-cs('pool3',m) 
# IBD patterns
tab(c(patt,'relation'),c(sum(qp),sum(qr)),zibd,
        c('none','1p2p','1m2m','1p2m','1m2p','1p1m','2p2m')[qp])
# correlated founding genes - suspect and alt suspect 
founder(spg)
founder(pool1) 
founder(pool2) 
founder(pool3)
which(aspg,patt,lw=c(2,1,2,2,2,2,2),spg,pool1)
which(smg,patt,lw=c(3,3,3,3,2,1,3),spg,aspg,pool2)
which(asmg,patt,lw=c(4,4,3,1,4,4,2),spg,aspg,smg,pool3) 
}

UAF<-function(M,m,g,Nc='c',Nd='d',Npool='pool',Ntemp='temp')
{
# M is Dirichlet parameter, m is marker suffix, g is vector
# of founder gene names (excluding suffix)
# NG correlated founding genes, Polya urn
# c's and d's are binary selection variables
# pool's and temp's are allele-valued
# for N = 2,3,...,NG
# cN is choice between a new draw from gene pool (poolN)
# and random choice among preceding g's (tempN1)
# d's create tempN1 from successive biased binary choices,
# with tempNM holding intermediate stages
# tempNj is random choice among first j genes out of N

NG<-length(g)

founder(cs(g[1],m))
for(i in 2:NG) founder(cs(Npool,i,m))

for(i in 2:NG) query(cs(Nc,i,m),c(i-1,M)/(M+i-1))

which(cs(g[2],m),cs(Nc,2,m),cs(g[1],m),cs(Npool,2,m))

for(i in 3:NG) 
{
for(j in 2:(i-1)) query(cs(Nd,i,j,m),c(j-1,1)/j)
which(cs(Ntemp,i,2,m),cs(Nd,i,2,m),cs(g[1],m),cs(g[2],m))
if(i>3) for(j in 3:(i-1)) 
	which(cs(Ntemp,i,j,m),cs(Nd,i,j,m),cs(Ntemp,i,j-1,m),cs(g[j],m))
which(cs(g[i],m),cs(Nc,i,m),cs(Ntemp,i,i-1,m),cs(Npool,i,m))
}

}

run<-function(model,use,pars) 
{ 
clean(T)

if(model=='IBD')
{
	alpha<-pars[1]; gamma<-pars[2]
	tab('relation',5,c(1-alpha-gamma,alpha/2,alpha/2,gamma/2,gamma/2),
        c('none','1pa2','2pa1','hsma','hspa'))
} else if(model=='UAF')
{
	M<-pars
} else stop('not a valid model')

query('guilty')

for(m in markers[use]) 
{ 
vs('alleles',alleles[[m]])
gene.freq<<-freqsC[[m]]

trace<-cs('trace',m)
smg<-cs('smg',m); spg<-cs('spg',m); sgt<-cs('sgt',m)
asmg<-cs('asmg',m); aspg<-cs('aspg',m); asgt<-cs('asgt',m)
patt<-cs('patt',m)

if(model=='IBD')
{
	IBD(patt,m,spg,smg,aspg,asmg)
} else if(model=='UAF')
{
	UAF(M,m,c('spg','smg','aspg','asmg'))
} 

genotype(sgt,spg,smg)
genotype(asgt,aspg,asmg)

which(trace,'guilty',sgt,asgt)
}

for(m in markers[use]) 
{
trace<-cs('trace',m); sgt<-cs('sgt',m)
prop.evid(trace,evidence[[m]])
prop.evid(sgt,evidence[[m]])
}

x<-marg(tcq[[acliq('guilty')]],'guilty')
cat(x[1]/x[2],':',markers[use],'\n')
invisible(x)

}

########################

# Computations for 'Baseline' column of Table 5 ###################

# Overall analysis using all 8 markers
run('IBD',1:8,c(0,0)) 

# Marker-by-marker analysis
prod<-1 
for(m in 1:8) 
	{x<-run('IBD',m,c(0,0)); prod<-prod*x[1]/x[2]} 

# Print 'product rule' log-likelihood ratio
log10(prod)

# Computations for 'UAF' column of Table 5 #######################

# Overall analysis using all 8 markers 
# (NB - not necessary, since product rule applies in UAF scenario)
run('UAF',1:8,100) 

# Marker-by-marker analysis
prod<-1 
for(m in 1:8) 
	{x<-run('UAF',m,100); prod<-prod*x[1]/x[2]} 

# Print 'product rule' log-likelihood ratio
log10(prod)

# Computations for 'IBD' column of Table 5 ########################

# Overall analysis using all 8 markers
run('IBD',1:8,c(.05,.05)) 

# Marker-by-marker analysis
prod<-1 
for(m in 1:8) 
	{x<-run('IBD',m,c(.05,.05)); prod<-prod*x[1]/x[2]} 

# Print 'product rule' log-likelihood ratio
log10(prod)