forked from xavierdidelot/TransPhyloMatlab
-
Notifications
You must be signed in to change notification settings - Fork 0
/
tutorial.m
31 lines (22 loc) · 822 Bytes
/
tutorial.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
%Initialisation of random number generators
rand('state',0);randn('state',0);
%Simulation of an outbreak
outbreak=simu(100/365,1000,0.001,1,10);
%Plot the outbreak as a colored phylogeny
plotBothTree(outbreak);
%Extract and plot just the transmission tree
plotTTree(ttreeFromFullTree(outbreak));
%Extract and plot the phylogeny
plotPTree(ptreeFromFullTree(outbreak));
%Inference of transmission tree given a phylogeny
ptree=ptreeFromFullTree(outbreak);
timeLastSeq=max(outbreak(:,1));
result=infer(ptree,10000,timeLastSeq);
%Trace of log-posterior
plot([result(:).pTTree]+[result(:).pPTree]);
%Trace of parameter beta
plot([result(:).beta]);
%Histogram of probabilities to be the source case
hist([result(end/2:end).source],1:10);
%Plot configuration after first half of the MCMC
plotBothTree(result(end/2).tree);