# Asia  
# set up DAG
clean(T)

query('asia',c(0.01,0.99))
query('smoke')
tab(c('tb','asia'),,c(.05,.95,.01,.99),c('yes','no'))
tab(c('cancer','smoke'),,c(.1,.9,.01,.99),c('yes','no'))
tab(c('bronc','smoke'),,c(.6,.4,.3,.7),c('yes','no'))
or('tbcanc','tb','cancer')
tab(c('xray','tbcanc'),,c(.98,.02,.05,.95),c('yes','no'))
tab(c('dysp','tbcanc','bronc'),,c(.9,.1,.8,.2,.7,.3,.1,.9),c('yes','no'))

pnmarg('cancer')

prop.evid('asia','yes')
prop.evid('dysp','yes')
prop.evid('xray','no')
pnmarg('cancer')

#-----------------------------------------------------------
# Simple illustration of Bayes theorem - blood tests  
# set up DAG
clean(T)

tab('disease',,c(0.95,0.05),c('harmless','serious'))
tab('number',4,rep(0.25,4))
p<-c(.04,.98)
P<-array(0,c(5,2,4))
for(i in 1:2) for(n in 1:4) P[,i,n]<-dbinom(0:4,n,p[i])
tab(c('test','disease','number'),c(5,2,4),P,0:4)

prop.evid('number',1)
prop.evid('test',1,usevs=T)
pnmarg('disease')

equil()
prop.evid('number',3)
prop.evid('test',1,usevs=T)
pnmarg('disease')

#-----------------------------------------------------------
# An ordinary Markov chain
clean(T)

tab(0,5,rep(0.2,5))
for(i in 1:20) tab(c(i,i-1),c(5,5),c(
 .7,.3,0,0,0,
 .3,.4,.3,0,0,
 0,.3,.4,.3,0,
 0,0,.3,.4,.3,
 0,0,0,.3,.7))
prop.evid(10,2)
fq(0:20,tr=T)


#-----------------------------------------------------------
# Mix  
# set up DAG
clean(T)

vs('alleles',c('8','10','11','x'))
gene.freq<<-c(.184884,.134884,.233721,.446511)
founder('vmg')
founder('vpg')
genotype('vgt','vmg','vpg')

founder('smg')
founder('spg')
genotype('sgt','smg','spg')

query('T2eqv',,c('V','U'))
query('T1eqs',,c('S','U'))
by('target','T2eqv','T1eqs')

select('T2mg','vmg','T2eqv')
select('T2pg','vpg','T2eqv')

select('T1mg','smg','T1eqs')
select('T1pg','spg','T1eqs')

genotype('T2gt','T2mg','T2pg')
genotype('T1gt','T1mg','T1pg')

mix('mix','T2gt','T1gt')

prop.evid('vgt','8-10')
prop.evid('sgt','8-11')
prop.evid('mix','8-10-11')
pnmarg('target')


#-----------------------------------------------------------
# Simple mutation model  
# set up DAG
clean(T)

vs('alleles',c('10','12','x'))
gene.freq<<-c(0.284,0.259,0.457)

founder('pfmg')
founder('pfpg')
genotype('pfgt','pfmg','pfpg')
query('tfeqpf')
select('tfmg','pfmg','tfeqpf')
select('tfpg','pfpg','tfeqpf')
mendel('copg','tfmg','tfpg')
founder('mmg')
founder('mpg')
genotype('mgt','mmg','mpg')
mendel('comg','mmg','mpg')
select('capg','copg','mutp')
select('camg','comg','mutm')
genotype('cgt','camg','capg')

# loop over values of mutation probability,
# defining rest of model, recompiling, initialising
# and equilibrating, and re-inserting and propagting evidence

pr.mut<-(0:10)/200
for(j in 1:length(pr.mut))
	{
	pr<-pr.mut[j]
	query('mutm',c(1-pr,pr),c('no','yes'))
	query('mutp',c(1-pr,pr),c('no','yes'))

	prop.evid('cgt','10-12',q=T)
	prop.evid('mgt','10-10',q=T)
	prop.evid('pfgt','12-12',q=T)

	cat(pr,nm('tfeqpf'),'\n')
	}


#-----------------------------------------------------------
# Mixtures with an unknown number of contributors
clean(T)

vs('alleles',c('A','B','C'))
gene.freq<<-c(.566,.429,.005)
nmax<-2

for(id in c('v','s',paste('u',1:nmax,sep='')))
      {founder(cs(id,'mg'));  founder(cs(id,'pg'))}
query('v.in.mix',,c('V','U'))
query('s.in.mix',,c('S','U'))
tab('n.unknown',nmax+1,rep(1/(nmax+1),nmax+1),nmax:0)
by('target','v.in.mix','s.in.mix','n.unknown')
for(al in vs.alleles)
      {
      for(id in c('v','s',paste('u',1:nmax,sep='')))
            or(cs(id,'has',al),cs(id,'mg'),cs(id,'pg'),
                  al==vs.alleles,al==vs.alleles)
      for(id in c('s','v'))
            and(cs(al,id),cs(id,'has',al),cs(id,'.in.mix'))
      for(n in 1:nmax)
            {
            id<-cs('u',n)
            and(cs(al,id),cs(id,'has',al),'n.unknown',,(nmax:0)>=n)
            }
      or(cs(al,'m',0),cs(al,'v'),cs(al,'s'))
      for(n in 1:nmax)
            or(cs(al,'m',n),cs(al,'m',n-1),cs(al,'u',n))
      }

prop.evid(cs('Am',nmax),'yes')
prop.evid(cs('Bm',nmax),'yes')
prop.evid(cs('Cm',nmax),'no')

prop.evid('vhasA','yes')
prop.evid('vhasB','yes')

prop.evid('shasA','yes')
prop.evid('shasB','no')
prop.evid('shasC','no')

pnmarg('target')
print(joint(c('v.in.mix','s.in.mix')))
fq('n.unknown',tr=T,va=T)