# Random-Tree simulation

I wanna simulate a Galton-Watson Tree to a maximum of n generation given a reproduction law P. I use Maple but I am unable to create the edges of the tree whenever there are more than two vertices in the tree. I think that the problem is that I am unable to store the number of children that each node has. Here is the code:

I use the law below to test, but I believe it should work with any discrete law on posint. P:=RandomVariable(EmpiricalDistribution(Array([0,1,2]))): Then the procedure to generate the tree is

GWT:=proc(n,P) #Renvoie un Arbre de Galton-Watson avec au plus n génération sous loi P;
local G,g,p,s,N,k;
if n=0 then
G:=Graph(1);
else
g:=GWT(n-1,P);
p:=Sample(P,1);
G:=DisjointUnion(Graph(1),$(g,p[1])); if p[1]=1 then s=0; else s:=((p[1])^n-1)/(p[1]-1); fi; N:=nops(Vertices(G)); G:=RelabelVertices(G,[seq(k,k=1..N)]); if p[1]>0 then G:=AddEdge(G,{seq({1, 2+k*s},k = 0 .. p[1]-1)}); fi; G; fi; end proc;  Just to compare, I was able to code the construction of a m-ary tree up to n generation, given n and m with a pretty similar code:  Arbre:=proc(n,m) # m-ary arbre à n génération local k,T,t,s,N; if n=0 then T:=Graph(1); else t:=Arbre(n-1,m); T:=DisjointUnion(Graph(1),$(t,m));
if m=1 then s:=0; else s:=(m^n-1)/(m-1);
fi;
N:=nops(Vertices(T));
T:=RelabelVertices(T,[seq(k,k=1..N)]);
T:=AddEdge(T,{seq({1, 2+k*s},k = 0 .. m-1)});
T;
fi;
end proc;


I was wondering if anyone could help with the coding? Maybe another way of storing the children. It may also be that it cannot be done with recursion.

Thanks

Edit#1 : Pseudo code would look like:

Input: n:nonneg integer, P: probability law
if n=0 then return 1 vertice
else
create 1 vertice; i=0; while i<=n
for each vertices in generation i generate a random number p from Law P and connect it
with p other vertices (in generation i+1); i=i+1;
Return the graph

-
Perhaps you should write the algorithm in pseudo code, in case the problem is there and other people can help you too. –  Patrick Da Silva May 15 '12 at 21:25

Here's my version, which seems to work well in Maple 16.

GWT:=proc(n,P)
uses Statistics, GraphTheory;
local v,v1,v2,gen,S,T1,T2,E,i;
v:= 1; E:= {};
v1:=1; v2:= 1;
for gen from 1 to n do
if v1 > v2 then break end if;
S:= map(round, Sample(P,v2-v1+1));
T2[v1-1]:= v2;
for i from v1 to v2 do
T1[i]:= T2[i-1]+1;
T2[i]:= T1[i]+S[i-v1+1]-1;
E:= E union {seq({i,j},j=T1[i] .. T2[i])};
end do;
v1:= T1[v1];
v2:= T2[v2];
end do;
Graph([\$1..v2],E);
end proc;


To try it out (with your P):

GraphTheory:-DrawGraph(GWT(5,P), style=tree, root=1);

-
Thanks Robert! I'll look it up tomorrow and try to break down what you've done. I tested it as well and got it to work with Probability laws that take only a finite numbers of values, but I couldnt get it to process with, say a Poisson or a Geometric distribution. –  Jean-Sébastien May 16 '12 at 3:45
Forget what I said about the laws, it works with everything I tried. Thanks again –  Jean-Sébastien May 16 '12 at 3:49