0
$\begingroup$

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 
  • 0
    Perhaps you should write the algorithm in pseudo code, in case the problem is there and other people can help you too.2012-05-15

1 Answers 1

0

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);

enter image description here

  • 0
    Forget what I said about the laws, it works with everything I tried. Thanks again2012-05-16