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
    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.2012-05-16
  • 0
    Forget what I said about the laws, it works with everything I tried. Thanks again2012-05-16