performance - Speeding up nested loop equation in Mathematica -

so i'm generating genetic drift model investigating change in variables {pab,pab,pab,pab} using recursion equation

psel[{pab_, pab_, pab_, pab_},s_] :=      ({pab, pab, pab, pab} {1, 1 + s, 1 + s, (1 + s)^2})     /(1 + s pab + s pab + (2 s + s^2) pab)  pdrift[{pab_, pab_, pab_, pab_}, npop_] :=     if[npop == infinity, {pab, pab, pab, pab},     randomvariate[multinomialdistribution[npop, {pab, pab, pab, pab}]]    /n[npop]];  driftequation[{pab_, pab_, pab_, pab_}, s_, npop_] :=     pdrift[psel[{pab, pab, pab, pab}, s], npop];

this run in for loop tmaxtimes. loop repeated reps times using do loop. imputed altogether function

simulatereps[s_, npop_, tmax_, pint_, reps_] :=        module[{i, pqs, pqsrun, t, driftequs},     pqs = {};     driftequs[{pab_, pab_, pab_, pab_}] =       simplify[driftequation[{pab, pab, pab, pab}, s, npop]];     do[pqsrun = table[{{0, 0, 0, 0}, {0, 0, 0, 0}}, {tmax}];      pqsrun[[1]] = pint;      for[t = 1, t < tmax, t++,        pqsrun[[t + 1]] = driftequs[pqsrun[[t]]]];      pqs = append[pqs, pqsrun];, {reps}];     return[pqs];     ];     pqs = simulatereps[0.001, 10000, 1000, {0.25, 0.25, 0.25, 0.25},10];

with pqs being final result.

this works fine, except takes ~15.982986 seconds run, , need run 9 variations of function. way of speeding loops?

i've read compiling code might work, i'm not sure how use compile function. i've tried reap , sow instead of append, couldn't work.

i'm still new coding basic tips appreciated.

ps: there error message after runningpqs = simulatereps[0.001, 10000, 1000,{0.25, 0.25, 0.25, 0.25},10];, doesn't seem effect results.