| 1 | /* |
|---|
| 2 | * ml3.c |
|---|
| 3 | * |
|---|
| 4 | * |
|---|
| 5 | * Part of TREE-PUZZLE 5.0 (June 2000) |
|---|
| 6 | * |
|---|
| 7 | * (c) 1999-2000 by Heiko A. Schmidt, Korbinian Strimmer, |
|---|
| 8 | * M. Vingron, and Arndt von Haeseler |
|---|
| 9 | * (c) 1995-1999 by Korbinian Strimmer and Arndt von Haeseler |
|---|
| 10 | * |
|---|
| 11 | * All parts of the source except where indicated are distributed under |
|---|
| 12 | * the GNU public licence. See http://www.opensource.org for details. |
|---|
| 13 | */ |
|---|
| 14 | |
|---|
| 15 | |
|---|
| 16 | #define EXTERN extern |
|---|
| 17 | |
|---|
| 18 | |
|---|
| 19 | /* prototypes */ |
|---|
| 20 | #include <stdio.h> |
|---|
| 21 | #include <stdlib.h> |
|---|
| 22 | #include <math.h> |
|---|
| 23 | #include "util.h" |
|---|
| 24 | #include "ml.h" |
|---|
| 25 | #include "gamma.h" |
|---|
| 26 | |
|---|
| 27 | |
|---|
| 28 | |
|---|
| 29 | /******************************************************************************/ |
|---|
| 30 | /* discrete Gamma-distribution and related stuff */ |
|---|
| 31 | /******************************************************************************/ |
|---|
| 32 | |
|---|
| 33 | /* compare general base frequencies with frequencies of taxon i with chi square */ |
|---|
| 34 | double homogentest(int taxon) |
|---|
| 35 | { |
|---|
| 36 | return chi2test(Freqtpm, Basecomp[taxon], gettpmradix(), &chi2fail); |
|---|
| 37 | } |
|---|
| 38 | |
|---|
| 39 | |
|---|
| 40 | /* discrete Gamma according to Yang 1994 (JME 39:306-314) */ |
|---|
| 41 | void YangDiscreteGamma (double shape, int c, dvector x) |
|---|
| 42 | { |
|---|
| 43 | double twoc, mu; |
|---|
| 44 | int i; |
|---|
| 45 | |
|---|
| 46 | twoc = 2.0*c; |
|---|
| 47 | mu = 0.0; |
|---|
| 48 | for (i = 0; i < c; i++) |
|---|
| 49 | { |
|---|
| 50 | /* corresponding rates */ |
|---|
| 51 | x[i] = icdfGamma ( (2.0*i+1.0)/twoc, shape); |
|---|
| 52 | mu += x[i]; |
|---|
| 53 | } |
|---|
| 54 | mu = mu/c; |
|---|
| 55 | |
|---|
| 56 | /* rescale for avarage rate of 1.0 */ |
|---|
| 57 | for (i = 0; i < c; i++) |
|---|
| 58 | { |
|---|
| 59 | x[i] /= mu; |
|---|
| 60 | } |
|---|
| 61 | } |
|---|
| 62 | |
|---|
| 63 | /* compute rates of each category when rates are Gamma-distributed */ |
|---|
| 64 | void updaterates() |
|---|
| 65 | { |
|---|
| 66 | int i; |
|---|
| 67 | double alpha; |
|---|
| 68 | |
|---|
| 69 | if (numcats == 1) |
|---|
| 70 | { |
|---|
| 71 | Rates[0] = 1.0; |
|---|
| 72 | return; |
|---|
| 73 | } |
|---|
| 74 | if (Geta == 0.0) |
|---|
| 75 | { |
|---|
| 76 | for (i = 0; i < numcats; i++) |
|---|
| 77 | Rates[i] = 1.0; |
|---|
| 78 | return; |
|---|
| 79 | } |
|---|
| 80 | alpha = (1.0 - Geta)/Geta; |
|---|
| 81 | |
|---|
| 82 | YangDiscreteGamma (alpha, numcats, Rates); |
|---|
| 83 | |
|---|
| 84 | /* if invariable sites are present */ |
|---|
| 85 | for (i = 0; i < numcats; i++) |
|---|
| 86 | Rates[i] = Rates[i]/(1.0-fracinv); |
|---|
| 87 | |
|---|
| 88 | /* check for very small rates */ |
|---|
| 89 | for (i = 0; i < numcats; i++) |
|---|
| 90 | if (Rates[i] < 0.000001) Rates[i] = 0.000001; |
|---|
| 91 | } |
|---|
| 92 | |
|---|
| 93 | |
|---|
| 94 | |
|---|
| 95 | /******************************************************************************/ |
|---|
| 96 | /* parameter estimation */ |
|---|
| 97 | /******************************************************************************/ |
|---|
| 98 | |
|---|
| 99 | /* compute sample mean and standard deviation of sample mean */ |
|---|
| 100 | void computestat(double *data, int n, double *mean, double *err) |
|---|
| 101 | { |
|---|
| 102 | int i; |
|---|
| 103 | double sum; |
|---|
| 104 | |
|---|
| 105 | sum = 0; |
|---|
| 106 | for (i = 0; i < n; i++) sum += data[i]; |
|---|
| 107 | (*mean) = sum/(double) n; |
|---|
| 108 | |
|---|
| 109 | sum = 0; |
|---|
| 110 | for (i = 0; i < n; i++) sum += (data[i] - (*mean))*(data[i] - (*mean)); |
|---|
| 111 | if (n != 1) |
|---|
| 112 | (*err) = sqrt(sum)/sqrt((double)(n-1)*n); /* unbiased estimator */ |
|---|
| 113 | else |
|---|
| 114 | (*err) = 0.0; /* if n == 1 */ |
|---|
| 115 | } |
|---|
| 116 | |
|---|
| 117 | /* compute ML value of quartet (a,b,c,d) */ |
|---|
| 118 | double quartetml(int a, int b, int c, int d) |
|---|
| 119 | { |
|---|
| 120 | double d1, d2, d3; |
|---|
| 121 | |
|---|
| 122 | /* compute ML for all topologies */ |
|---|
| 123 | if (approxp_optn) { /* approximate parameter mode */ |
|---|
| 124 | d1 = quartet_alklhd(a,b,c,d); /* (a,b)-(c,d) */ |
|---|
| 125 | d2 = quartet_alklhd(a,c,b,d); /* (a,c)-(b,d) */ |
|---|
| 126 | d3 = quartet_alklhd(a,d,b,c); /* (a,d)-(b,c) */ |
|---|
| 127 | } else { |
|---|
| 128 | d1 = quartet_lklhd(a,b,c,d); /* (a,b)-(c,d) */ |
|---|
| 129 | d2 = quartet_lklhd(a,c,b,d); /* (a,c)-(b,d) */ |
|---|
| 130 | d3 = quartet_lklhd(a,d,b,c); /* (a,d)-(b,c) */ |
|---|
| 131 | } |
|---|
| 132 | |
|---|
| 133 | /* looking for max(d1, d2, d3) */ |
|---|
| 134 | if (d1 < d2) { /* d2 > d1 */ |
|---|
| 135 | if (d2 < d3) { /* d3 > d2 > d1 */ |
|---|
| 136 | /* d3 maximum */ |
|---|
| 137 | return d3; |
|---|
| 138 | } else { /* d2 >= d3 > d1 */ |
|---|
| 139 | /* d2 maximum */ |
|---|
| 140 | return d2; |
|---|
| 141 | } |
|---|
| 142 | } else { /* d1 >= d2 */ |
|---|
| 143 | if (d1 < d3) { /* d3 > d1 >= d2 */ |
|---|
| 144 | /* d3 maximum */ |
|---|
| 145 | return d3; |
|---|
| 146 | } else { /* d1 >= d2 && d1 >= d3 */ |
|---|
| 147 | /* d1 maximum */ |
|---|
| 148 | return d1; |
|---|
| 149 | } |
|---|
| 150 | } |
|---|
| 151 | } |
|---|
| 152 | |
|---|
| 153 | /* optimization function TSparam - quartets */ |
|---|
| 154 | double opttsq(double x) |
|---|
| 155 | { |
|---|
| 156 | if (x < MINTS) TSparam = MINTS; |
|---|
| 157 | else if (x > MAXTS) TSparam = MAXTS; |
|---|
| 158 | else TSparam = x; |
|---|
| 159 | tranprobmat(); |
|---|
| 160 | distupdate(qca, qcb, qcc, qcd); |
|---|
| 161 | return (-quartetml(qca, qcb, qcc, qcd)); |
|---|
| 162 | } |
|---|
| 163 | |
|---|
| 164 | /* optimization function YRparam - quartets */ |
|---|
| 165 | double optyrq(double x) |
|---|
| 166 | { |
|---|
| 167 | if (x < MINYR) YRparam = MINYR; |
|---|
| 168 | else if (x > MAXYR) YRparam = MAXYR; |
|---|
| 169 | else YRparam = x; |
|---|
| 170 | tranprobmat(); |
|---|
| 171 | distupdate(qca, qcb, qcc, qcd); |
|---|
| 172 | return (-quartetml(qca, qcb, qcc, qcd)); |
|---|
| 173 | } |
|---|
| 174 | |
|---|
| 175 | /* estimate substitution process parameters - random quartets */ |
|---|
| 176 | void optimseqevolparamsq() |
|---|
| 177 | { |
|---|
| 178 | double tsmeanold, yrmeanold; |
|---|
| 179 | dvector tslist, yrlist; |
|---|
| 180 | int fin; |
|---|
| 181 | ivector taxon; |
|---|
| 182 | uli minqts, maxqts, n; |
|---|
| 183 | |
|---|
| 184 | |
|---|
| 185 | taxon = new_ivector(4); |
|---|
| 186 | |
|---|
| 187 | /* number of quartets to be investigated */ |
|---|
| 188 | minqts = (uli) floor(0.25 * MINPERTAXUM * Maxspc) + 1; |
|---|
| 189 | maxqts = (uli) floor(0.25 * MAXPERTAXUM * Maxspc) + 1; |
|---|
| 190 | if (Maxspc == 4) { |
|---|
| 191 | minqts = (uli) 1; |
|---|
| 192 | maxqts = (uli) 1; |
|---|
| 193 | } |
|---|
| 194 | |
|---|
| 195 | tslist = new_dvector(maxqts); |
|---|
| 196 | yrlist = new_dvector(maxqts); |
|---|
| 197 | |
|---|
| 198 | /* initialize averages */ |
|---|
| 199 | tsmean = TSparam; |
|---|
| 200 | yrmean = YRparam; |
|---|
| 201 | |
|---|
| 202 | fin = FALSE; |
|---|
| 203 | |
|---|
| 204 | /* investigate maxqts random quartets */ |
|---|
| 205 | for (n = 0; n < maxqts; n++) { |
|---|
| 206 | |
|---|
| 207 | /* choose random quartet */ |
|---|
| 208 | chooser(Maxspc, 4, taxon); |
|---|
| 209 | |
|---|
| 210 | /* |
|---|
| 211 | * optimize parameters on this quartet |
|---|
| 212 | */ |
|---|
| 213 | |
|---|
| 214 | qca = taxon[0]; |
|---|
| 215 | qcb = taxon[1]; |
|---|
| 216 | qcc = taxon[2]; |
|---|
| 217 | qcd = taxon[3]; |
|---|
| 218 | |
|---|
| 219 | /* initialize start values with average value */ |
|---|
| 220 | if ((SH_optn || nuc_optn) && optim_optn && (data_optn == 0)) TSparam = tsmean; |
|---|
| 221 | if ((nuc_optn && TN_optn) && optim_optn && (data_optn == 0)) YRparam = yrmean; |
|---|
| 222 | |
|---|
| 223 | /* estimation */ |
|---|
| 224 | twodimenmin(PEPS1, |
|---|
| 225 | (SH_optn || nuc_optn) && optim_optn && (data_optn == 0), |
|---|
| 226 | MINTS, &TSparam, MAXTS, opttsq, &tserr, |
|---|
| 227 | (nuc_optn && TN_optn) && optim_optn && (data_optn == 0), |
|---|
| 228 | MINYR, &YRparam, MAXYR, optyrq, &yrerr); |
|---|
| 229 | |
|---|
| 230 | |
|---|
| 231 | tsmeanold = tsmean; |
|---|
| 232 | yrmeanold = yrmean; |
|---|
| 233 | tslist[n] = TSparam; |
|---|
| 234 | yrlist[n] = YRparam; |
|---|
| 235 | computestat(tslist, n+1 , &tsmean, &tserr); |
|---|
| 236 | computestat(yrlist, n+1 , &yrmean, &yrerr); |
|---|
| 237 | |
|---|
| 238 | /* check whether the means are converging */ |
|---|
| 239 | if (n > minqts-2) { |
|---|
| 240 | if ((fabs(tsmean-tsmeanold) < TSDIFF) && |
|---|
| 241 | (fabs(yrmean-yrmeanold) < YRDIFF)) |
|---|
| 242 | fin = TRUE; |
|---|
| 243 | } |
|---|
| 244 | |
|---|
| 245 | /* investigate at least minqts quartets */ |
|---|
| 246 | if (n > minqts-2 && (fin || n > maxqts-2)) break; |
|---|
| 247 | } |
|---|
| 248 | |
|---|
| 249 | /* round estimated numbers to 2 digits after the decimal point */ |
|---|
| 250 | if (tserr != 0.0) tsmean = floor(100.0*tsmean+0.5)/100.0; |
|---|
| 251 | if (yrerr != 0.0) yrmean = floor(100.0*yrmean+0.5)/100.0; |
|---|
| 252 | |
|---|
| 253 | /* update ML engine */ |
|---|
| 254 | TSparam = tsmean; |
|---|
| 255 | YRparam = yrmean; |
|---|
| 256 | tranprobmat(); |
|---|
| 257 | |
|---|
| 258 | free_ivector(taxon); |
|---|
| 259 | } |
|---|
| 260 | |
|---|
| 261 | /* optimization function TSparam - tree */ |
|---|
| 262 | double opttst(double x) |
|---|
| 263 | { |
|---|
| 264 | double result; |
|---|
| 265 | |
|---|
| 266 | if (x < MINTS) TSparam = MINTS; |
|---|
| 267 | else if (x > MAXTS) TSparam = MAXTS; |
|---|
| 268 | else TSparam = x; |
|---|
| 269 | tranprobmat(); |
|---|
| 270 | computedistan(); |
|---|
| 271 | if (approxp_optn) result = usertree_alklhd(); |
|---|
| 272 | else result = usertree_lklhd(); |
|---|
| 273 | |
|---|
| 274 | return (-result); |
|---|
| 275 | } |
|---|
| 276 | |
|---|
| 277 | /* optimization function YRparam - tree */ |
|---|
| 278 | double optyrt(double x) |
|---|
| 279 | { |
|---|
| 280 | double result; |
|---|
| 281 | |
|---|
| 282 | if (x < MINYR) YRparam = MINYR; |
|---|
| 283 | else if (x > MAXYR) YRparam = MAXYR; |
|---|
| 284 | else YRparam = x; |
|---|
| 285 | tranprobmat(); |
|---|
| 286 | computedistan(); |
|---|
| 287 | if (approxp_optn) result = usertree_alklhd(); |
|---|
| 288 | else result = usertree_lklhd(); |
|---|
| 289 | |
|---|
| 290 | return (-result); |
|---|
| 291 | } |
|---|
| 292 | |
|---|
| 293 | |
|---|
| 294 | /* optimize substitution process parameters - tree */ |
|---|
| 295 | void optimseqevolparamst() |
|---|
| 296 | { |
|---|
| 297 | twodimenmin(PEPS1, |
|---|
| 298 | (SH_optn || nuc_optn) && optim_optn && (data_optn == 0), |
|---|
| 299 | MINTS, &TSparam, MAXTS, opttst, &tserr, |
|---|
| 300 | (nuc_optn && TN_optn) && optim_optn && (data_optn == 0), |
|---|
| 301 | MINYR, &YRparam, MAXYR, optyrt, &yrerr); |
|---|
| 302 | } |
|---|
| 303 | |
|---|
| 304 | |
|---|
| 305 | /* optimization function fracinv */ |
|---|
| 306 | double optfi(double x) |
|---|
| 307 | { |
|---|
| 308 | double result; |
|---|
| 309 | |
|---|
| 310 | if (x < MINFI) fracinv = MINFI; |
|---|
| 311 | else if (x > MAXFI) fracinv = MAXFI; |
|---|
| 312 | else fracinv = x; |
|---|
| 313 | |
|---|
| 314 | computedistan(); |
|---|
| 315 | if (approxp_optn) result = usertree_alklhd(); |
|---|
| 316 | else result = usertree_lklhd(); |
|---|
| 317 | |
|---|
| 318 | return (-result); |
|---|
| 319 | } |
|---|
| 320 | |
|---|
| 321 | |
|---|
| 322 | /* optimization function Geta */ |
|---|
| 323 | double optge(double x) |
|---|
| 324 | { |
|---|
| 325 | double result; |
|---|
| 326 | |
|---|
| 327 | if (x < MINGE) Geta = MINGE; |
|---|
| 328 | else if (x > MAXGE) Geta = MAXGE; |
|---|
| 329 | else Geta = x; |
|---|
| 330 | |
|---|
| 331 | updaterates(); |
|---|
| 332 | |
|---|
| 333 | computedistan(); |
|---|
| 334 | if (approxp_optn) result = usertree_alklhd(); |
|---|
| 335 | else result = usertree_lklhd(); |
|---|
| 336 | |
|---|
| 337 | return (-result); |
|---|
| 338 | } |
|---|
| 339 | |
|---|
| 340 | |
|---|
| 341 | /* optimize rate heterogeneity parameters */ |
|---|
| 342 | void optimrateparams() |
|---|
| 343 | { |
|---|
| 344 | twodimenmin(PEPS2, |
|---|
| 345 | fracinv_optim, |
|---|
| 346 | MINFI, &fracinv, fracconst, optfi, &fierr, |
|---|
| 347 | grate_optim, |
|---|
| 348 | MINGE, &Geta, MAXGE, optge, &geerr); |
|---|
| 349 | |
|---|
| 350 | } |
|---|
