-/* rTrait.c 2010-07-26 */
+/* rTrait.c 2011-06-25 */
-/* Copyright 2010 Emmanuel Paradis */
+/* Copyright 2010-2011 Emmanuel Paradis */
/* This file is part of the R-package `ape'. */
/* See the file ../COPYING for licensing issues. */
}
break;
case 2 : for (i = *Nedge - 1; i >= 0; i--) {
+ if (alpha[i]) {
+ alphaT = alpha[i] * el[i];
+ M = exp(-alphaT);
+ S = sigma[i] * sqrt((1 - exp(-2 * alphaT))/(2 * alpha[i]));
+ } else {
+ M = 1;
+ S = sqrt(el[i]) * sigma[i];
+ }
GetRNGstate();
- x[edge2[i]] = x[edge1[i]] - alpha[i]*el[i]*(x[edge1[i]] - theta[i]) + sqrt(el[i])*sigma[i]*norm_rand();
- PutRNGstate();
- }
- break;
- case 3 : for (i = *Nedge - 1; i >= 0; i--) {
- GetRNGstate();
- x[edge2[i]] = x[edge1[i]] - (1 - exp(alpha[i]*el[i]))*(x[edge1[i]] - theta[i]) + sqrt(el[i])*sigma[i]*norm_rand();
+ x[edge2[i]] = x[edge1[i]] * M + theta[i] * (1 - M) + S * norm_rand();
PutRNGstate();
}
break;