PageRenderTime 19ms CodeModel.GetById 11ms app.highlight 6ms RepoModel.GetById 1ms app.codeStats 0ms

/working/auteur.extended/R/splitormerge.R

http://github.com/eastman/auteur
R | 72 lines | 61 code | 9 blank | 2 comment | 10 complexity | bbe5a97172a9d6b607b6ffc5133bf185 MD5 | raw file
 1#rjmcmc proposal mechanism: proposal mechanism for traversing model complexity (by one parameter at a time)
 2#author: JM EASTMAN 2010
 3
 4splitormerge <-
 5function(cur.delta, cur.values, phy, node.des, lambda=lambda, theta=FALSE, internal.only=FALSE) { 
 6	bb=cur.delta
 7	vv=cur.values
 8	names(vv)<-names(bb)<-phy$edge[,2]
 9	new.bb=choose.one(bb, phy=phy, internal.only)
10	new.vv=vv
11	
12	s=names(new.bb[bb!=new.bb])
13	all.shifts=as.numeric(names(bb[bb>0]))
14	all.D=node.des[[which(names(node.des)==s)]]
15	if(length(all.D)!=0) {
16		untouchable=unlist(lapply(all.shifts[all.shifts>s], function(x)node.des[[which(names(node.des)==x)]]))
17		remain.unchanged=union(all.shifts, untouchable)
18	} else {
19		untouchable=NULL
20		remain.unchanged=list()
21	}
22	
23	marker=match(s, names(new.vv))
24	nn=length(vv)
25	K=sum(bb)
26	N=Ntip(phy)
27	
28	ncat=sum(bb)
29	cur.vv=as.numeric(vv[marker])
30	ca.vv=length(which(vv==cur.vv))
31	
32	if(sum(new.bb)>sum(bb)) {			# add transition: SPLIT
33		decision="split"
34		n.desc=sum(!all.D%in%remain.unchanged)+1
35		n.split=sum(vv==cur.vv)-n.desc
36		if(theta) u=splitvalue(cur.vv=cur.vv, n.desc=n.desc, n.split=n.split, factor=lambda) else u=splitrate(cur.vv=cur.vv, n.desc=n.desc, n.split=n.split, factor=lambda)
37		nr.split=u$nr.split
38		nr.desc=u$nr.desc
39		new.vv[vv==cur.vv]=nr.split
40		if(length(remain.unchanged)==0) {	# assign new value to all descendants
41			new.vv[match(all.D,names(new.vv))] = nr.desc
42		} else {							# assign new value to all 'open' descendants 
43			new.vv[match(all.D[!(all.D%in%remain.unchanged)],names(new.vv))] = nr.desc
44		}
45		new.vv[match(s, names(new.vv))]=nr.desc
46		lnHastingsRatio = log((K+1)/(2*N-2-K)) ### from Drummond and Suchard 2010: where N is tips, K is number of local parms in tree
47		lnPriorRatio = log(ptpois(K+1,lambda,nn)/ptpois(K,lambda,nn))
48		
49	} else {							# drop transition: MERGE
50		decision="merge"
51		anc = get.ancestor.of.node(s, phy)
52		if(!is.root(anc, phy)) {			# base new rate on ancestral rate of selected branch
53			anc.vv=as.numeric(vv[match(anc,names(vv))])
54			na.vv=length(which(vv==anc.vv))
55			nr=(anc.vv*na.vv+cur.vv*ca.vv)/(ca.vv+na.vv)
56			new.vv[vv==cur.vv | vv==anc.vv]=nr
57		} else {							# if ancestor of selected node is root, base new rate on sister node
58			sister.tmp=get.desc.of.node(anc,phy)
59			sister=sister.tmp[sister.tmp!=s]
60			sis.vv=as.numeric(vv[match(sister,names(vv))])
61			ns.vv=length(which(vv==sis.vv))
62			nr=(sis.vv*ns.vv+cur.vv*ca.vv)/(ca.vv+ns.vv)
63			new.vv[vv==cur.vv | vv==sis.vv]=nr			
64		}
65		lnHastingsRatio = log((2*N-2-K+1)/K) ### from Drummond and Suchard 2010: where N is tips, K is number of local parms in tree
66		lnPriorRatio = log(ptpois(K-1,lambda,nn)/ptpois(K,lambda,nn))
67		
68	}
69	
70	return(list(new.delta=new.bb, new.values=new.vv, lnHastingsRatio=lnHastingsRatio, lnPriorRatio=lnPriorRatio, decision=decision))
71}
72