Ulrich,
Martijn gave you some great advice on speeding up your loops. I also suggest that you use SeedRandom[]. If you use SeedRandom you can get the same result every time. This helps with timing the routine and verifying that you do not inadvertently add a bug when improving your code.
For example I can do this and get consistent timing and always the same results:
SeedRandom[3]; AbsoluteTiming[
distributionF1C[3, 5, 3, 10000, 10]]
{0.082987, {{1., 0.261}, {2., 0.329}, {3., 0.3443}, {4.,
0.3518}, {5., 0.3516}, {6., 0.3493}, {7., 0.349}, {8.,
0.3476}, {9., 0.3475}}}
You also asked
I could use the code
r = RandomInteger[{0, d}, n];
s = n - Count[r, 0];
only for the first generation because the lineages that do not die out must be preserved. I don't see how I could use it for generations V>1
The best way to avoid looping is to think in terms of lists. For example, you currently represent a generation with a list:
SeedRandom[3]; first = RandomInteger[3, 4]
Which is {1, 3, 2, 2}
You can then generate a new list with the next generation
Map[RandomInteger[3, #] &, first]
to get {{2}, {0, 0, 2}, {2, 0}, {0, 0}} but you only care about the total in each lineage so use this instead
Map[Total[RandomInteger[3, #]] &, first]
to get your next generation, {2, 2, 2, 0}
Now you can do this for V generations and get a list of every generation on the way using NestList
SeedRandom[3]; NestList[
Function[x, Map[Total[RandomInteger[3, #]] &, x]],
RandomInteger[3, 4], 8 - 2]
to get every generation (7 of them) total for each lineage:
{{1, 3, 2, 2}, {2, 2, 2, 0}, {3, 5, 4, 0}, {7, 7, 2, 0}, {13, 7, 5,
0}, {23, 14, 6, 0}, {43, 23, 14, 0}}
I used V-2 because your code only did 7 generations for entering V=8 and NestList includes the first one automatically so V=1 gives 2 generations.
Now I just need to count up the totals that equal S and run a bunch of samples and divide by the number of samples so here is the final code:
mydist[d_, n_, S_, z_, V_] :=
Total[Table[
Map[Boole[(n - Count[#, 0]) == S] &,
NestList[Function[x, Map[Total[RandomInteger[d, #]] &, x]],
RandomInteger[d, n], V - 2]], z]]/N[z]
To test it I get the same answer as the compiled loop above in 10,000 samples but its 17% faster than the compiled loop version.
SeedRandom[3]; AbsoluteTiming[mydist[3, 5, 3, 10000, 10]]
{0.070848, {0.261, 0.329, 0.3443, 0.3518, 0.3516, 0.3493,
0.349, 0.3476, 0.3475}}
I tried to compile this loop and it did not speed up at all -- in fact its slightly slower! (I am not sure why but maybe MMA is already fast). Maybe Martijn or someone else might see if I did something wrong or if this code can be improved for compiling. Here is that code
mydistc :=
Compile[{{d, _Integer}, {n, _Integer}, {S, _Integer}, {z, _Integer}, \
{V, _Integer}},
Total[Table[
Map[Boole[(n - Count[#, 0]) == S] &,
NestList[Function[x, Map[Total[RandomInteger[d, #]] &, x]],
RandomInteger[d, n], V - 2]], z]]/N[z]]
I hope this helps,
Regards,
Neil
