1 /* mh.c
2
3 Written by Frederic Bois
4
5 Copyright (c) 1996-2018 Free Software Foundation, Inc.
6
7 This file is part of GNU MCSim.
8
9 GNU MCSim is free software; you can redistribute it and/or
10 modify it under the terms of the GNU General Public License
11 as published by the Free Software Foundation; either version 3
12 of the License, or (at your option) any later version.
13
14 GNU MCSim is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18
19 You should have received a copy of the GNU General Public License
20 along with GNU MCSim; if not, see <http://www.gnu.org/licenses/>
21
22 This file calls at least one GNU Scientific Library function, if the
23 symbol HAVE_LIBGSL is defined in config.h (or in the makefile).
24 Otherwise, the corresponding features are disabled (and the program
25 exits with an error message).
26 */
27
28 #include <stdio.h>
29 #include <string.h>
30 #include <inttypes.h>
31 #include <time.h>
32 #include <sys/time.h>
33
34 #ifdef USEMPI
35 #include <mpi.h>
36 #endif
37
38 #include "mh.h"
39 #include "lsodes.h"
40 #include "lexerr.h"
41 #include "simmonte.h"
42
43 #ifdef HAVE_LIBGSL
44 #include <gsl/gsl_cdf.h>
45 #endif
46
47
48 /* Function -------------------------------------------------------------------
49 CalculateMeanAndVariance
50
51 Calculate mean and variance of each component of array x, incrementally.
52 Index parameter n must be received as 0 at the start of a sequence.
53 */
54 #ifdef NDEF
55 /* Frederic's version */
CalculateMeanAndVariance(long * n,double x,double * mean,double * var)56 void CalculateMeanAndVariance (long *n, double x, double *mean, double *var)
57 {
58 static double K = 0.0;
59 static double Ex = 0.0;
60 static double Ex2 = 0.0;
61
62 double dTmp;
63
64 /* shifted data algorithm */
65 if (*n == 0) {
66 K = x;
67 Ex = 0.0;
68 Ex2 = 0.0;
69 }
70 else {
71 dTmp = x - K;
72 Ex += dTmp;
73 Ex2 += dTmp * dTmp;
74 }
75
76 (*n)++;
77
78 dTmp = Ex / *n;
79 *mean = K + dTmp;
80
81 *var = (Ex2 - Ex * dTmp) / (*n - 1);
82
83 #ifdef NDEF
84 // or Welford algorithm:
85
86 // for a new value newValue, compute the new count, new mean, the new M2.
87 // mean accumulates the mean of the entire dataset
88 // M2 aggregates the squared distance from the mean
89 // count aggregates the number of samples seen so far
90 def update(existingAggregate, newValue):
91 (count, mean, M2) = existingAggregate
92 count = count + 1
93 delta = newValue - mean
94 mean = mean + delta / count
95 delta2 = newValue - mean
96 M2 = M2 + delta * delta2
97
98 return (count, mean, M2)
99
100 # retrieve the mean, variance and sample variance from an aggregate
101 def finalize(existingAggregate):
102 (count, mean, M2) = existingAggregate
103 (mean, variance, sampleVariance) = (mean, M2/count, M2/(count - 1))
104 if count < 2:
105 return float('nan')
106 else:
107 return (mean, variance, sampleVariance)
108 #endif
109
110 } /* CalculateMeanAndVariance */
111 #endif
112
113 /* Cathie's version
114 * Variance is calculated using a method described by Knuth
115 * m_k = m_k-1 + (x_k - m_k-1)/k
116 * v_k = v_k-1 + (x_k - m_k-1)(x_k - m_k)
117 * sigma^2 = v_k / (k-1)
118 * This portion of the code only calculates the itermediate
119 * values - before the variance is sent along the final
120 * calculation should be executed (CollectConvInfo)
121 */
122 void CalculateMeanAndVariance (long n, double x, double *xi_bari,
123 double *si_2i) {
124
125 if (n == 1) {
126 *xi_bari = x; // mean
127 *si_2i = 0;
128 return;
129 }
130
131 double mTmp = *xi_bari;
132 *xi_bari = *xi_bari + (x - *xi_bari)/n;
133 *si_2i = *si_2i + (x - mTmp) * (x - *xi_bari);
134
135 } /* End CalculateMeanAndVariance */
136
137
138 int checkConvergence (int nOut, int variableCount, int p_count,
139 double **meansForAll, double **varsForAll, double *Rhat)
140 {
141 int vi, pi;
142 int converged = 0;
143 double varsofvars, meansofmeans, varsofmeans, meansofvars;
144
145 // first calculate the means x_bar = Mean(nChains, xi_bar);
146 // calculate the mean of the variances W = Mean (nChains, si_2);
147 for (vi = 0; vi < variableCount; vi++) {
148 meansofmeans = 0.0;
149 for (pi = 0; pi < p_count; pi++) {
150 // use the function we already wrote for consistency
151 CalculateMeanAndVariance((pi+1),meansForAll[pi][vi],&meansofmeans,
152 &varsofmeans);
153 CalculateMeanAndVariance((pi+1),varsForAll[pi][vi],&meansofvars,
154 &varsofvars);
155 }
156
157 if ((meansofvars == 0) && (varsofmeans == 0)) {
158 *Rhat = 1;
159 converged++;
160 } else {
161 // calculate Vhat and Rhat together (to save memory)
162 // s_2 = ((nOut - 1) * W + B) / nOut; // intra-chain variance part 1
163 double s2 = ((nOut-1) * meansofvars + varsofmeans) / nOut;
164 // Vhat = s_2 + B / (nOut * dChains); // part 2
165 double Vhat = s2 + varsofmeans /(nOut * p_count);
166 // Rhat = Vhat / W; // convergence criterion
167 *Rhat = Vhat / meansofvars;
168 if (*Rhat < 1.05) // we should parameterize this 1.05 criterion
169 converged++;
170 }
171 }
172
173 return converged;
174
175 } /* end checkConvergence */
176
177
178 void CollectConvInfo (PLEVEL plevel, char **args)
179 {
180 double **mean_dest = (double **) args[0];
181 double **var_dest = (double **) args[1];
182 long *n = (long *) args[2];
183 long i;
184 PMCVAR pMCVar;
185
186 /* For all MC vars at this level */
187 for (i = 0; i < plevel->nMCVars; i++) {
188 pMCVar = plevel->rgpMCVars[i];
189 **mean_dest = pMCVar->dVal_mean;
190 /* Keep in mind that the running variance calculation is an
191 intermediate value and we want to communicate the variance here
192 That means that we need to do the final step before sending
193 sigma^2 = v_k / (k-1) */
194 **var_dest = pMCVar->dVal_var/(*n-1);
195 (*mean_dest)++;
196 (*var_dest)++;
197 }
198 } /* end CollectConvInfo */
199
200
201 /* Function -------------------------------------------------------------------
202 AnnounceMarkov
203
204 Print out the type of simulations to be performed.
205 */
206 void AnnounceMarkov (int size, int nSimTypeFlag, long nIter)
207 {
208
209 switch (nSimTypeFlag) {
210
211 case 0:
212 printf("\nDoing %ld Metropolis within Gibbs simulation", nIter);
213 printf((nIter != 1 ? "s" : ""));
214 if (size > 1)
215 printf(" on each of %d processors\n", size);
216 else
217 printf("\n");
218 break;
219
220 case 1:
221 printf("\nPrinting data and predictions for the last line of the "
222 "restart file\n");
223 break;
224
225 case 2:
226 printf("\nDoing %ld Metropolis-Hastings simulation", nIter);
227 printf((nIter != 1 ? "s" : ""));
228 if (size > 1)
229 printf(" on each of %d processors\n", size);
230 else
231 printf("\n");
232 break;
233
234 case 3:
235 printf("\nDoing %ld Metropolis within Gibbs posterior "
236 "tempered simulation", nIter);
237 printf((nIter != 1 ? "s\n" : "\n"));
238 break;
239
240 case 4:
241 printf("\nDoing %ld Metropolis within Gibbs likelihood "
242 "tempered simulation", nIter);
243 printf((nIter != 1 ? "s\n" : "\n"));
244 break;
245
246 case 5:
247 printf("\nDoing Stochastic optimization\n");
248 break;
249 }
250
251 } /* AnnounceMarkov */
252
253
254 /* Function -------------------------------------------------------------------
255 CalculateTotals
256
257 Find total prior, likelihood for all MC vars
258 Called from TraverseLevels
259 */
260 void CalculateTotals (PLEVEL plevel, char **args)
261 {
262 PANALYSIS panal = (PANALYSIS)args[0];
263 double *pdLnPrior = (double*)args[1];
264
265 long n;
266
267 for (n = 0; n < plevel->nMCVars; n++) {
268 *pdLnPrior += LnDensity(plevel->rgpMCVars[n], panal);
269 }
270
271 } /* CalculateTotals */
272
273
274 /* Function -------------------------------------------------------------------
275 CheckForFixed
276
277 It is possible for an MC var to be fixed by a subsequent `=' statement in
278 the level just below the level at which it is declared in the input file;
279 This routine marks such vars as "Fixed" and set their dVal to the
280 prescribed number.
281 Called from TraverseLevels
282 */
283
284 void CheckForFixed (PLEVEL plevel, char **args)
285 {
286 long n, m;
287 PMCVAR pMCVar;
288 PVARMOD pFVar;
289
290 for (n = 0; n < plevel->nMCVars; n++) {
291 pMCVar = plevel->rgpMCVars[n];
292 for (m = 0; m < plevel->nFixedVars; m++) {
293 pFVar = plevel->rgpFixedVars[m];
294 if (pMCVar->hvar == pFVar->hvar) {
295 pMCVar->bIsFixed = TRUE;
296 if (IsInput (pFVar->hvar)) {
297 printf("Error: a sampled parameter cannot be assigned an input\n");
298 exit(0);
299 }
300 else
301 pMCVar->dVal = pFVar->uvar.dVal;
302 }
303 }
304 }
305
306 } /* CheckForFixed */
307
308
309 /* Function -------------------------------------------------------------------
310 CheckPrintStatements
311
312 Tries to check that 'Print' and 'Data' statements are consistent.
313 Note that experiments must have outputs; this is checked in
314 PrepareOutSpec in siminit.c.
315 Called from TraverseLevels
316 */
317
318 void CheckPrintStatements (PLEVEL plevel, char **args)
319 {
320 PANALYSIS panal = (PANALYSIS) args[0];
321 POUTSPEC pos;
322 /* PMCVAR pMCVar;
323 BOOL bFound, bOK; */
324 long i, j /* , k, l, dCount, dIndex */;
325
326 if (plevel->pexpt == NULL)
327 return;
328
329 pos = &(plevel->pexpt->os);
330
331 /* check that the same var does not appear in 2 or more print statements */
332 for (j = 0; j < pos->nOutputs; j++)
333 for (i = j+1; i < pos->nOutputs; i++)
334 if (pos->phvar_out[j] == pos->phvar_out[i])
335 ReportRunTimeError (panal, RE_DUPVARINEXPRT | RE_FATAL,
336 pos->pszOutputNames[i], "Print");
337
338 /* check that the same var does not appear in 2 or more data statements */
339 for (j = 0; j < pos->nData; j++)
340 for (i = j+1; i < pos->nData; i++)
341 if (pos->phvar_dat[j] == pos->phvar_dat[i])
342 ReportRunTimeError (panal, RE_DUPVARINEXPRT | RE_FATAL,
343 pos->pszDataNames[i], "Data");
344
345 /* check that print and matching data lists are of equal length */
346 for (j = 0; j < pos->nOutputs; j++)
347 for (i = 0; i < pos->nData; i++)
348 if ((pos->phvar_out[j] == pos->phvar_dat[i]) &&
349 (pos->pcOutputTimes[j] != pos->pcData[i])) {
350 printf("\nError: unequal times in Print and Data statements for %s\n"
351 "Exiting.\n\n", pos->pszOutputNames[j]);
352 exit(0);
353 }
354
355 } /* CheckPrintStatements */
356
357
358 /* Function -------------------------------------------------------------------
359 CheckAllTransitions
360
361 Check whether all perk transitions rates are high enough.
362 */
363 BOOL CheckAllTransitions (PGIBBSDATA pgd)
364 {
365 BOOL bOK;
366 double AcceptRate;
367 int i;
368
369 i = pgd->startT;
370 bOK = TRUE;
371 while ((i <= pgd->endT - 1) && bOK) {
372 if (pgd->rglTransAttempts[i] < 10) {
373 bOK = FALSE;
374 break;
375 }
376 else
377 AcceptRate = pgd->rglTransAccepts[i] /
378 ((double) pgd->rglTransAttempts[i]);
379
380 bOK = (AcceptRate > 0.15);
381 i++;
382 }
383
384 return bOK;
385
386 } /* end CheckAllTransitions */
387
388
389 /* Function -------------------------------------------------------------------
390 CheckTransitions
391
392 Check whether the first perk transition rate is between reasonable bounds
393 and that enough attemps have been made
394 Return -1 if too low or too few attempts, 0 if OK, +1 if too high
395 */
396 int CheckTransitions (PGIBBSDATA pgd)
397 {
398 double AcceptRate;
399 int i;
400
401 i = pgd->startT;
402 if (pgd->rglTransAttempts[i] < 10)
403 return (-1); /* too low */
404 else
405 AcceptRate = pgd->rglTransAccepts[i] /
406 ((double) pgd->rglTransAttempts[i]);
407
408 if (AcceptRate < 0.30) {
409 return (-1); /* too low */
410 }
411 else {
412 if (AcceptRate < 1) {
413 return (0); /* OK */
414 }
415 else {
416 return (+1); /* too high */
417 }
418 }
419
420 } /* end CheckTransitions */
421
422
423 /* Function -------------------------------------------------------------------
424 CloneLikes
425
426 Called from TraverseLevels
427 First copy the likelihoods in the current plistLikes in the rgpLikes
428 down to the next lower levels. Then eventually override them by definitions
429 in rgpLikes at this level.
430 */
431 void CloneLikes (PLEVEL plevel, char **args)
432 {
433 long nLikes;
434 long i, j, k;
435 PLEVEL pLower;
436 PMCVAR pClone;
437 PMCVAR pMCVar;
438 BOOL bFound;
439
440 for (i = 0; i < plevel->iInstances; i++) {
441 pLower = plevel->pLevels[i];
442 /* The number of likelihoods is overestimated by this, but that will do
443 for now: */
444 pLower->nLikes = nLikes = plevel->nLikes + ListLength(plevel->plistLikes);
445
446 if (pLower->nLikes != 0) {
447 if (!(pLower->rgpLikes =
448 (PMCVAR*) malloc (pLower->nLikes * sizeof(PMCVAR))))
449 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "CloneLikes", NULL);
450 }
451 }
452
453 /* Copy the likelihoods in the plistLikes at this level to the
454 rgpLikes of the lower levels */
455 nLikes = 0;
456 ForAllList3 (plevel->plistLikes, &CloneLikesL, plevel, &nLikes, NULL);
457
458 /* Note: nLikes now equals the number of likelihoods in the current
459 plistLikes */
460
461 /* Copy the likelihoods in rgpLikes at this level down to the lower
462 levels, if they are not already defined */
463 for (i = 0; i < plevel->iInstances; i++) {
464 pLower = plevel->pLevels[i];
465
466 for (j = 0; j < plevel->nLikes; j++) {
467 pMCVar = plevel->rgpLikes[j];
468
469 /* if that likelihood is already in pLower->rgpLikes forget it */
470 bFound = FALSE;
471 k = 0;
472 while ((k < nLikes) && (!bFound)) {
473 bFound = (pMCVar->hvar == pLower->rgpLikes[k]->hvar);
474 if (!bFound) k++;
475 }
476 if (!bFound) {
477 if (!(pClone = (PMCVAR) malloc (sizeof (MCVAR))))
478 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "CloneLikes", NULL);
479
480 memcpy (pClone, pMCVar, sizeof (MCVAR));
481 pLower->rgpLikes[nLikes+j] = pClone;
482 }
483 }
484 }
485
486 } /* CloneLikes */
487
488
489 /* Function -------------------------------------------------------------------
490 CloneLikesL
491
492 Clone and copy the members of a plistLikes from a level to the rgpLikes
493 arrays of the levels below
494 Called from ForAllList3
495 */
496 void CloneLikesL (PVOID pData, PVOID pUser1, PVOID pUser2, PVOID pUser3)
497 {
498 PMCVAR pMCVar = (PMCVAR) pData;
499 PLEVEL plevel = (PLEVEL) pUser1;
500 long *pnLikes = (long *) pUser2;
501 long i;
502 PLEVEL pLower;
503 PMCVAR pClone;
504
505 ++pMCVar->iDepth;
506 for (i = 0; i < plevel->iInstances; i++) {
507 pLower = plevel->pLevels[i];
508 if (!(pClone = (PMCVAR) malloc (sizeof (MCVAR))))
509 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "CloneLikeL", NULL);
510 memcpy (pClone, pMCVar, sizeof (MCVAR));
511 pLower->rgpLikes[*pnLikes] = pClone;
512 }
513 ++(*pnLikes);
514
515 } /* CloneLikesL */
516
517
518 /* Function -------------------------------------------------------------------
519 CloneMCVars
520
521 Called from TraverseLevels
522 For all MC vars in list at given level, add to arrays of all instances of
523 next (lower) level; the instances are created by the cloning routine
524 CloneMCVarsL
525 */
526 void CloneMCVars (PLEVEL plevel, char **args)
527 {
528 long nMCVars = ListLength(plevel->plistMCVars);
529 long n;
530 PLEVEL pLower;
531
532 for (n = 0; n < plevel->iInstances; n++) {
533 pLower = plevel->pLevels[n];
534 pLower->nMCVars = nMCVars;
535 if (nMCVars != 0) {
536 if (!(pLower->rgpMCVars = (PMCVAR*) malloc (nMCVars * sizeof(PMCVAR))))
537 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "CloneMCVars", NULL);
538 }
539 }
540
541 nMCVars = 0;
542 ForAllList3 (plevel->plistMCVars, &CloneMCVarsL, plevel, &nMCVars, NULL);
543
544 } /* CloneMCVars */
545
546
547 /* Function -------------------------------------------------------------------
548 CloneMCVarsL
549
550 Clone and copy the members of a plistMCVars
551 Called from ForAllList3
552 */
553 void CloneMCVarsL (PVOID pData, PVOID pUser1, PVOID pUser2, PVOID pUser3)
554 {
555 PMCVAR pMCVar = (PMCVAR) pData;
556 PLEVEL plevel = (PLEVEL) pUser1;
557 long *pnMCVars = (long *) pUser2;
558 long i;
559 PLEVEL pLower;
560 PMCVAR pClone;
561
562 ++pMCVar->iDepth;
563 for (i = 0; i < plevel->iInstances; i++) {
564 pLower = plevel->pLevels[i];
565 if (!(pClone = (PMCVAR) malloc(sizeof (MCVAR))))
566 ReportError(NULL, RE_OUTOFMEM | RE_FATAL, "CloneMCVarsL", NULL);
567 memcpy (pClone, pMCVar, sizeof (MCVAR));
568 pClone->plistDependents = InitList();
569 pLower->rgpMCVars[*pnMCVars] = pClone;
570 }
571 ++(*pnMCVars);
572
573 } /* CloneMCVarsL */
574
575
576 /* Function -------------------------------------------------------------------
577 CloseMarkovFiles
578
579 Closes output files associated with the Markov sampler.
580 The restart file has already been closed by the ReadRestart
581
582 */
583 void CloseMarkovFiles (PGIBBSDATA pgd)
584 {
585 /* If tempered MCMC, close perk scale recording file */
586 if ((pgd->nSimTypeFlag == 3) || (pgd->nSimTypeFlag == 4)) {
587 char szFileName[MAX_FILENAMESIZE+6];
588 sprintf(szFileName, "%s%s", pgd->szGout, ".perks");
589 fclose(pgd->pfilePerks);
590 printf("\nWrote perks to \"%s\"\n", szFileName);
591 }
592
593 if (pgd->pfileOut) {
594 fclose(pgd->pfileOut);
595 printf("Wrote MCMC sample to \"%s\"\n", pgd->szGout);
596 }
597
598 printf("\n");
599
600 } /* CloseMarkovFiles */
601
602
603 /* Function -------------------------------------------------------------------
604 ConvertLists
605
606 Converts lists to arrays
607 Called from TraverseLevels
608 */
609 void ConvertLists(PLEVEL plevel, char **args)
610 {
611 PANALYSIS panal = (PANALYSIS)args[0];
612 long m, n;
613 PMCVAR pMCVar;
614
615 if (plevel->pexpt == NULL)
616 ListToPVArray (panal, plevel->plistVars, &plevel->nFixedVars,
617 &plevel->rgpFixedVars);
618 else
619 ListToPVArray (panal, plevel->pexpt->plistParmMods, &plevel->nFixedVars,
620 &plevel->rgpFixedVars);
621
622 for (n = 0; n < plevel->nMCVars; n++) {
623 pMCVar = plevel->rgpMCVars[n];
624 ListToPMCArray (panal, pMCVar->plistDependents,
625 &pMCVar->nDependents, &pMCVar->rgpDependents);
626
627 /* if there are no dependent, experiments are most likely to be dependent
628 on that parameter (complete checking would require verifying the
629 model file - FB 2 May 1998 */
630 if (pMCVar->nDependents == 0)
631 pMCVar->bExptIsDep = TRUE;
632 else {
633 /* if there are dependent, experiments may be dependent on that
634 parameter if they are not shielded by an intermediate instance of
635 the same parameter */
636 m = 0;
637 while ((m < pMCVar->nDependents) &&
638 !(pMCVar->bExptIsDep = (strcmp(pMCVar->rgpDependents[m]->pszName,
639 pMCVar->pszName) ? TRUE : FALSE)))
640 m++;
641 }
642 }
643
644 } /* ConvertLists */
645
646
647 /* Function -------------------------------------------------------------------
648 DoMarkov
649
650 Core routine of the MCMC sampler
651 */
652 void DoMarkov (PANALYSIS panal)
653 {
654 #ifdef USEMPI
655 //double startMarkov = time(NULL);
656 #endif
657 PGIBBSDATA pgd = &panal->gd;
658 PLEVEL pLevel0 = panal->pLevels[0];
659 long nThetas, nUpdateAt, nTotal;
660 long i, iter = 0;
661 long nIter = pgd->nMaxIter; /* scheduled iterations of the sampler */
662 double *pdMCVarVals = NULL; /* read in values of thetas */
663 double *pdSum = NULL; /* column sums of read thetas */
664 double **prgdSumProd = NULL; /* column sums of thetas cross-products */
665 double dTmp, dLnPrior = 0, dLnData = 0;
666
667 if (panal->rank == 0)
668 AnnounceMarkov(panal->size, pgd->nSimTypeFlag, nIter);
669
670 OpenMarkovFiles(panal);
671
672 ReadDataFile(panal); /* (if needed) */
673
674 /* MC variables must be placed in arrays at the next lower level */
675 TraverseLevels(pLevel0, CloneMCVars, NULL);
676
677 /* likelihoods must percolate down to the experiment level */
678 TraverseLevels(pLevel0, CloneLikes, NULL);
679
680 /* find the parents and dependents of the MC vars */
681 TraverseLevels(pLevel0, FindMCParents, panal, NULL);
682 TraverseLevels(pLevel0, FindMCDependents, panal, NULL);
683
684 /* find the parents of the parameters in likelihood statements */
685 TraverseLevels(pLevel0, FindLikeParents, panal, NULL);
686
687 /* convert the rest of the lists to arrays */
688 TraverseLevels(pLevel0, ConvertLists, panal, NULL);
689
690 /* check for MC vars that have been fixed */
691 TraverseLevels(pLevel0, CheckForFixed, NULL);
692
693 /* check variables in statements of type
694 'Distrib(<var1>, <distrib>, Prediction, <var2>)'
695 for identical 'Print' statements */
696 TraverseLevels(pLevel0, CheckPrintStatements, panal, NULL);
697
698 /* if required, print out the structure for debugging */
699 if ((panal->rank == 0) && (panal->bDependents)) {
700 printf("Hierarchical structure:\n\n");
701 TraverseLevels(pLevel0, PrintDeps, NULL);
702 printf("\nDone.\n\n");
703 return;
704 }
705
706 /* change the MC vars hvar pointers from pointing to model parameters to
707 pointing to the parents' dVal */
708 TraverseLevels(pLevel0, SetPointers, panal, NULL);
709
710 /* get the initial values of the MC vars by reading the restart file if
711 one is given */
712 if (pgd->szGrestart) {
713
714 /* count the variables to read */
715 nThetas = 0;
716 TraverseLevels(pLevel0, GetNumberOfMCVars, &nThetas);
717
718 /* read the starting values in the order they are printed and
719 close the file when finished */
720 if (((pgd->nSimTypeFlag == 3) || (pgd->nSimTypeFlag == 4)) &&
721 (pgd->nPerks != 0)) {
722 /* if the user has required tempering and specified a perk scale,
723 then read also the pseudo-priors */
724 ReadRestartTemper(pgd->pfileRestart, nThetas, pgd->nPerks,
725 &pdMCVarVals, &pdSum, &prgdSumProd, &iter,
726 &pgd->indexT, pgd->rgdlnPi);
727 }
728 else
729 ReadRestart(pgd->pfileRestart, nThetas, &pdMCVarVals, &pdSum,
730 &prgdSumProd, &iter);
731
732 /* set the dVals of the variables to the values read in */
733 nThetas = 0; /* variables will be recounted */
734 if (!TraverseLevels1 (pLevel0, SetMCVars, pdMCVarVals, &nThetas, NULL)) {
735 printf("\nError: some read-in parameters are out of bounds - "
736 "Exiting\n\n");
737 exit(0);
738 }
739
740 /* If nSimTypeFlag is 1 print just the predictions and data and exit */
741 if (pgd->nSimTypeFlag == 1) {
742 if (panal->rank == 0) { /* do it on only one processor */
743 PrintAllExpts(pLevel0, panal, pgd->pfileOut);
744 CloseMarkovFiles (pgd);
745 }
746 return;
747 }
748
749 /* If component jumps, tempered or stochastic optimization: read eventually
750 the kernel file */
751 if ((pgd->nSimTypeFlag == 0) || (pgd->nSimTypeFlag >= 3)) {
752
753 char szKernelFile[MAX_FILENAMESIZE+12];
754 /* prefix the filename with the rank of the process if more than one
755 process is used, postfix it with .kernel in any case */
756 if (panal->size > 1) {
757 sprintf(szKernelFile, "%04d_%s%s", panal->rank, panal->gd.szGrestart,
758 ".kernel");
759 }
760 else {
761 sprintf(szKernelFile, "%s%s", panal->gd.szGrestart, ".kernel");
762 }
763
764 FILE *pfile = fopen(szKernelFile, "r");
765
766 if (pfile) {
767 /* Read kernel sizes from file */
768 printf("Reading kernel file %s\n", szKernelFile);
769 TraverseLevels(pLevel0, ReadKernel, pfile, NULL);
770 }
771 else {
772 /* Set the jumping kernel's SD automatically */
773 TraverseLevels(pLevel0, SetKernel, 1, pdMCVarVals, NULL);
774
775 /* Free the now useless array */
776 free(pdMCVarVals);
777 }
778 }
779 else { /* vector jumping, initialize prior */
780 TraverseLevels(pLevel0, CalculateTotals, panal, &dLnPrior, NULL);
781 }
782
783 /* Initialize the predictions arrays by running all experiments */
784 if (!RunAllExpts (panal, &dLnData)) {
785 /* error, cannot compute */
786 printf("\nError: cannot compute at the starting point - Exiting\n\n");
787 exit(0);
788 }
789
790 /* Save the data likelihoods */
791 TraverseLevels1 (pLevel0, SaveLikelihoods, NULL);
792
793 /* If tempered MCMC and no user-specified perks, set the perk
794 scale and the pseudo-priors automatically */
795 InitPerks(panal);
796
797 /* Now that we have the MC vars etc. right, write the output file header */
798 WriteHeader(panal);
799
800 } /* end if restart file */
801
802 else { /* no restart file, init by sampling */
803
804 /* If nSimTypeFlag is 1 print an error message and exit */
805 if (pgd->nSimTypeFlag == 1) {
806 printf("\nError: a restart file must be given to print data and"
807 " predictions - Exiting.\n\n");
808 exit(0);
809 }
810
811 /* Set the jumping kernel's SD */
812 TraverseLevels(pLevel0, SetKernel, 2, pdMCVarVals, NULL);
813
814 /* initialize the thetas by sampling from the prior */
815 TraverseLevels(pLevel0, InitMCVars, NULL);
816
817 #ifdef USEMPI
818 /* count the variables if we monitor their convergence */
819 nThetas = 0;
820 TraverseLevels(pLevel0, GetNumberOfMCVars, &nThetas);
821 #endif
822
823 /* Calculate the total prior */
824 TraverseLevels(pLevel0, CalculateTotals, panal, &dLnPrior, NULL);
825
826 /* Initialize the predictions arrays by running all experiments */
827 if (!RunAllExpts (panal, &dLnData)) {
828 /* Error, cannot compute */
829 printf("\nError: cannot compute at the starting point - Exiting\n\n");
830 exit(0);
831 }
832
833 /* Save the data likelihoods */
834 TraverseLevels1(pLevel0, SaveLikelihoods, NULL);
835
836 /* If tempered MCMC and no user-specified perks, set the perk
837 scale and the pseudo-priors automatically */
838 InitPerks(panal);
839
840 /* Now that we have the MC vars etc. right, write the output file header */
841 WriteHeader(panal);
842
843 /* Write the current thetas to the output file */
844 fprintf(pgd->pfileOut, "0\t");
845 TraverseLevels(pLevel0, WriteMCVars, pgd->pfileOut, NULL);
846
847 /* Output indexT and pseudo-priors if needed */
848 if ((pgd->nSimTypeFlag == 3) || (pgd->nSimTypeFlag == 4)) {
849 fprintf(pgd->pfileOut, "%d\t", pgd->indexT);
850 for (i = 0; i < pgd->nPerks; i++)
851 fprintf(pgd->pfileOut, "%e\t", pgd->rgdlnPi[i]);
852 }
853
854 /* Output prior and likelihood */
855 fprintf(pgd->pfileOut, "%e\t%e\t%e\n", dLnPrior, dLnData,
856 dLnPrior + dLnData);
857 fflush(pgd->pfileOut);
858
859 } /* end no restart file */
860
861 /* Initializations are finished, let's do the iterations */
862 nTotal = UPDATE_BASE;
863 nUpdateAt = iter + nTotal; /* kernel will be updated at that iteration */
864
865 #ifdef USEMPI
866 double **meansForAll;
867 double **varsForAll;
868 double *means;
869 double *vars;
870 double Rhat;
871 if (panal->bPrintConvergence) {
872 /* when running in parallel we need to pack the running mean and
873 convergence values for each variable into an array that will be
874 sent to the root rank; this code allocates the space for those values;
875 we also have to prep the 0 rank processor to calculate the variance
876 of variances and mean of means for convergence */
877 if (panal->rank == 0) {
878 meansForAll = (double**) malloc(sizeof(double*) * panal->size);
879 varsForAll = (double**) malloc(sizeof(double*) * panal->size);
880 int sender;
881 for (sender = 0; sender < panal->size; sender++){
882 meansForAll[sender] = (double*) malloc(sizeof(double) * nThetas);
883 varsForAll[sender] = (double*) malloc(sizeof(double) * nThetas);
884 }
885 means = meansForAll[0];
886 vars = varsForAll[0];
887 } else {
888 means = (double*) malloc(sizeof(double) * nThetas);
889 vars = (double*) malloc(sizeof(double) * nThetas);
890 }
891 }
892 //fprintf(stderr, "starting %ld iterations\n", nIter);
893 //double startIter = time(NULL);
894 //fprintf(stderr, "Time to warm up %lf\n", startIter - startMarkov);
895
896 MPI_Barrier(MPI_COMM_WORLD);
897 #endif
898
899 while (iter < nIter) {
900
901 #ifdef USEMPI
902 if (panal->bPrintConvergence) { /* convergence check */
903 if ((iter+1) % pgd->nPrintFreq == 0) {
904 double *tempMeans, *tempVars;
905 tempMeans = means;
906 tempVars = vars;
907 /* all processors must collect their info into an array */
908 TraverseLevels(pLevel0, CollectConvInfo, &tempMeans, &tempVars, &iter);
909 /* the root process needs to collect everyone's data and process it */
910 if (panal->rank == 0 ) {
911 int sender = 0;
912 MPI_Status status;
913 for (sender = 1; sender < panal->size; sender++) {
914 MPI_Recv(meansForAll[sender], nThetas, MPI_DOUBLE, sender,
915 0, MPI_COMM_WORLD, &status);
916 }
917 int converged = checkConvergence(iter+1, nThetas, panal->size,
918 meansForAll, varsForAll, &Rhat);
919 if (iter > 10) /* avoid pathologies at start */
920 fprintf(stderr, "Iteration %ld,\tRhat %g,\tconverged %d\n",
921 iter + 1, Rhat, converged);
922 if (iter == nIter - 1)
923 printf("\n");
924 /* everyone else just sends their data */
925 } else {
926 MPI_Send(means, nThetas, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
927 }
928 }
929 }
930 #endif
931
932 /* Output to screen, if required as a command-line option */
933 if (panal->bOutputIter && ((iter+1) % panal->nOutputFreq == 0)) {
934 if (panal->size > 1)
935 printf("Processor %d, Iteration %ld\n", panal->rank, iter + 1);
936 else
937 printf("Iteration %ld\n", iter + 1);
938 if (iter == nIter - 1)
939 printf("\n");
940 }
941
942 /* Start output to file, eventually */
943 if (((iter + 1) % pgd->nPrintFreq == 0) &&
944 (iter >= pgd->nMaxIter - pgd->nPrintIter))
945 fprintf(pgd->pfileOut, "%ld\t", iter + 1);
946
947 /* take this path for 0 and 5 */
948 if ((pgd->nSimTypeFlag == 0) || (pgd->nSimTypeFlag == 5)) {
949 /* component by component jumps or stochastic optimization */
950
951 TraverseLevels(pLevel0, SampleThetas, panal, pgd, &iter, &nUpdateAt,
952 &nTotal, NULL);
953
954 /* Output log-densities, eventually */
955 if (((iter + 1) % pgd->nPrintFreq == 0) &&
956 (iter >= pgd->nMaxIter - pgd->nPrintIter)) {
957 dLnPrior = 0.0;
958 TraverseLevels(pLevel0, CalculateTotals, panal, &dLnPrior, NULL);
959 dLnData = 0.0;
960 TraverseLevels1(pLevel0, SumAllExpts, &dLnData, NULL);
961 fprintf(pgd->pfileOut, "%e\t%e\t%e\n", dLnPrior, dLnData,
962 dLnPrior + dLnData);
963 fflush (pgd->pfileOut);
964 }
965 }
966 else {
967 /* this is the path for 3 and 4 */
968 if ((pgd->nSimTypeFlag == 3) || (pgd->nSimTypeFlag == 4)) {
969 /* tempered */
970 TraverseLevels(pLevel0, SampleThetasTempered, panal, pgd, &iter,
971 &nUpdateAt, &nTotal, &pgd->indexT, NULL);
972
973 dLnPrior = 0.0;
974 TraverseLevels(pLevel0, CalculateTotals, panal, &dLnPrior, NULL);
975 dLnData = 0.0;
976 TraverseLevels1(pLevel0, SumAllExpts, &dLnData, NULL);
977
978 /* Robbins-Monro updating of the pseudo prior */
979 dTmp = pgd->dCZero / (iter + pgd->dNZero);
980 for (i = 0; i < pgd->nPerks; i++) {
981 if (i == pgd->indexT)
982 pgd->rgdlnPi[i] -= dTmp;
983 else
984 pgd->rgdlnPi[i] += dTmp / pgd->nPerks;
985 }
986
987 /* Output indexT and log-densities, eventually */
988 if (((iter + 1) % pgd->nPrintFreq == 0) &&
989 (iter >= pgd->nMaxIter - pgd->nPrintIter)) {
990 fprintf(pgd->pfileOut, "%d\t", pgd->indexT);
991 for (i = 0; i < pgd->nPerks; i++)
992 fprintf(pgd->pfileOut, "%e\t", pgd->rgdlnPi[i]);
993 fprintf(pgd->pfileOut, "%e\t%e\t%e\n", dLnPrior, dLnData,
994 dLnPrior + dLnData);
995 fflush(pgd->pfileOut);
996 }
997
998 /* update population count of current temperature */
999 pgd->rglCount[pgd->indexT] = pgd->rglCount[pgd->indexT]+1;
1000
1001 /* test the temperature and change indexT if necessary */
1002 pgd->indexT = SampleTemperature2 (pgd, dLnPrior, dLnData);
1003
1004 } /* end tempered */
1005 else { /* vector jumps (option 2) */
1006
1007 SampleThetaVector(pLevel0, panal, nThetas, pdMCVarVals, pdSum,
1008 prgdSumProd, iter, nUpdateAt, nTotal, &dLnPrior,
1009 &dLnData);
1010
1011 /* Output, eventually */
1012 if (((iter + 1) % pgd->nPrintFreq == 0) &&
1013 (iter >= pgd->nMaxIter - pgd->nPrintIter)) {
1014
1015 for (i = 0; i < nThetas; i++)
1016 fprintf(pgd->pfileOut, "%5g\t", pdMCVarVals[i]);
1017
1018 fprintf(pgd->pfileOut, "%e\t%e\t%e\n", dLnPrior, dLnData,
1019 dLnPrior + dLnData);
1020 fflush(pgd->pfileOut);
1021 }
1022 }
1023 }
1024
1025 /* Adjust the update time, eventually */
1026 if (iter == nUpdateAt) {
1027 nTotal = (nTotal * 3) / 2;
1028 nUpdateAt = iter + nTotal;
1029 }
1030
1031 /* Increment the iteration counter */
1032 iter = iter + 1;
1033
1034 } /* while iter */
1035
1036 #ifdef USEMPI
1037 //double stopIter = time(NULL);
1038 //fprintf(stderr,"Time Iterating: %lf\n", (stopIter-startIter));
1039 #endif
1040
1041 /* If tempered MCMC: print perk scale */
1042 if ((pgd->nSimTypeFlag == 3) || (pgd->nSimTypeFlag == 4)) {
1043 /* PrintTemperatureDiagnostics (stdout, pgd); */
1044 PrintTemperatureDiagnostics(pgd->pfilePerks, pgd);
1045 }
1046
1047 /* If component jumps, tempered or stochastic optimization: save kernel */
1048 if ((pgd->nSimTypeFlag == 0) || (pgd->nSimTypeFlag >= 3)) {
1049
1050 char szKernelFile[MAX_FILENAMESIZE+7];
1051 sprintf(szKernelFile, "%s%s", panal->gd.szGout, ".kernel");
1052
1053 FILE *pfile = fopen (szKernelFile, "w");
1054 if (!pfile) {
1055 printf("Cannot create kernel saving file '%s'\n", panal->gd.szGdata);
1056 exit(0);
1057 }
1058
1059 /* Write out the jumping kernel SDs */
1060 TraverseLevels (pLevel0, WriteKernel, pfile, NULL);
1061
1062 fprintf(pfile, "\n");
1063 fclose(pfile);
1064 printf("Wrote kernel SDs to \"%s\"\n", szKernelFile);
1065 }
1066
1067 CloseMarkovFiles(pgd);
1068
1069 #ifdef USEMPI
1070 //double stopMarkov = time(NULL);
1071 //fprintf(stderr,"Time Markov: %lf\n", (stopMarkov - startMarkov));
1072 #endif
1073
1074 } /* DoMarkov */
1075
1076
1077 /* Function -------------------------------------------------------------------
1078 EqualSlopes
1079
1080 Check whether three points in a row are aligned (within a tolerance range).
1081 Return TRUE if they do, FALSE otherwise.
1082 The relative difference is tested in fact.
1083 */
1084 BOOL EqualSlopes (PDOUBLE x, PDOUBLE y, int i)
1085 {
1086 #define SL_EPSILON 0.01
1087 double s1, s2;
1088
1089 s1 = (y[i+1] - y[i]) / (x[i+1] - x[i]);
1090
1091 s2 = (y[i+2] - y[i]) / (x[i+2] - x[i]);
1092
1093 return (fabs(s2/s1 - 1) < SL_EPSILON);
1094
1095 } /* EqualSlopes */
1096
1097
1098 /* Function -------------------------------------------------------------------
1099 Extrapolate
1100
1101 Extrapolate linearly to dTargetX from a pair of
1102 (perk, pseudo-prior) points.
1103 */
1104 double Extrapolate (PGIBBSDATA pgd, double dTargetX, int i1, int i2)
1105 {
1106 return (pgd->rgdlnPi[i1] -
1107 (pgd->rgdPerks[i1] - dTargetX) *
1108 (pgd->rgdlnPi [i2] - pgd->rgdlnPi [i1]) /
1109 (pgd->rgdPerks[i2] - pgd->rgdPerks[i1]));
1110
1111 } /* Extrapolate */
1112
1113
1114 /* Function -------------------------------------------------------------------
1115 FindLikeParents
1116
1117 Called from TraverseLevels
1118 Find the parents of the MC vars corresponding to likelihood at this level
1119 by looking at sampled parameters at this and all previous levels. Note that
1120 dependencies are not checked among likelihoods.
1121 */
1122 void FindLikeParents (PLEVEL plevel, char **args)
1123 {
1124 PANALYSIS panal = (PANALYSIS)args[0];
1125 long k, l, m, n;
1126 PLEVEL pPrevLev;
1127 PMCVAR pMCVar1, pMCVar2;
1128 BOOL bFound;
1129
1130 /* Set the current level array as we pass through */
1131 panal->pCurrentLevel[plevel->iDepth] = plevel;
1132
1133 for (k = 0; k < plevel->nLikes; k++) {
1134 pMCVar1 = plevel->rgpLikes[k]; /* current likelihood */
1135
1136 for (l = 0; l < 4; l++) {
1137 if (pMCVar1->iParmType[l] == MCVP_PARM) {
1138 /* if there is a parent, find it, but parents must appear before
1139 current pMCVar */
1140 bFound = FALSE;
1141
1142 /* First, this level */
1143 for (m = 0; m < plevel->nMCVars; m++) { /* for all previous distrib */
1144 pMCVar2 = plevel->rgpMCVars[m];
1145 if (pMCVar1->hParm[l] == pMCVar2->hvar) {
1146 pMCVar1->pMCVParent[l] = pMCVar2;
1147 bFound = TRUE;
1148 }
1149 } /* for m */
1150
1151 /* Now, all previous levels */
1152 if (!bFound) {
1153 for (n = plevel->iDepth - 1; n >= 0; n--) {
1154 pPrevLev = panal->pCurrentLevel[n];
1155
1156 for (m = 0; m < pPrevLev->nMCVars; m++) {
1157 pMCVar2 = pPrevLev->rgpMCVars[m];
1158 if (pMCVar1->hParm[l] == pMCVar2->hvar) {
1159 pMCVar1->pMCVParent[l] = pMCVar2;
1160 bFound = TRUE;
1161 }
1162 } /* for m */
1163 } /* for n */
1164 }
1165
1166 if (!bFound) { /* oops, parent not found, error */
1167 printf("\n"
1168 "Error: parent in position %ld of %s must be\n"
1169 " declared before it when creating\n"
1170 " sampling dependencies - Exiting.\n\n",
1171 l, pMCVar1->pszName);
1172 exit(0);
1173 }
1174
1175 }
1176 } /* for l */
1177 } /* for k */
1178
1179 } /* FindLikeParents */
1180
1181
1182 /* Function -------------------------------------------------------------------
1183 FindMCDependents
1184
1185 Called from TraverseLevels
1186 Find the direct dependents of the MC vars at this level by looking at this
1187 and all lower levels
1188 */
1189 void FindMCDependents (PLEVEL plevel, char **args)
1190 {
1191 long i, j;
1192 PMCVAR pMCVar;
1193
1194 for (i = 0; i < plevel->nMCVars; i++) {
1195 pMCVar = plevel->rgpMCVars[i];
1196 for (j = 0; j < 4; j++)
1197 if ((pMCVar->pMCVParent[j] != NULL) &&
1198 (pMCVar->pMCVParent[j]->hvar == pMCVar->hParm[j]))
1199 QueueListItem(pMCVar->pMCVParent[j]->plistDependents, pMCVar);
1200 }
1201
1202 } /*FindMCDependents */
1203
1204
1205 /* Function -------------------------------------------------------------------
1206 FindMCParents
1207
1208 Called from TraverseLevels
1209 Find the parents of the MC vars at this level by looking at this and
1210 all previous levels.
1211 */
1212 void FindMCParents (PLEVEL plevel, char **args)
1213 {
1214 PANALYSIS panal = (PANALYSIS)args[0];
1215 long k, l, m, n;
1216 PLEVEL pPrevLev;
1217 PMCVAR pMCVar1, pMCVar2;
1218 BOOL bFound;
1219
1220 /* Set the current level array as we pass through */
1221 panal->pCurrentLevel[plevel->iDepth] = plevel;
1222
1223 for (k = 0; k < plevel->nMCVars; k++) {
1224 pMCVar1 = plevel->rgpMCVars[k]; /* current distrib */
1225
1226 for (l = 0; l < 4; l++) {
1227 if (pMCVar1->iParmType[l] == MCVP_PARM) {
1228 /* if there is a parent, find it, but parents must appear before
1229 current pMCVar */
1230 bFound = FALSE;
1231
1232 /* First, this level */
1233 for (m = 0; m < k; m++) { /* for all previous distrib */
1234 pMCVar2 = plevel->rgpMCVars[m];
1235 if (pMCVar1->hParm[l] == pMCVar2->hvar) {
1236 pMCVar1->pMCVParent[l] = pMCVar2;
1237 bFound = TRUE;
1238 }
1239 } /* for m */
1240
1241 /* Now, all previous levels */
1242 if (!bFound) {
1243 for (n = plevel->iDepth - 1; n >= 0; n--) {
1244 pPrevLev = panal->pCurrentLevel[n];
1245
1246 for (m = 0; m < pPrevLev->nMCVars; m++) {
1247 pMCVar2 = pPrevLev->rgpMCVars[m];
1248 if (pMCVar1->hParm[l] == pMCVar2->hvar) {
1249 pMCVar1->pMCVParent[l] = pMCVar2;
1250 bFound = TRUE;
1251 }
1252 } /* for m */
1253 } /* for n */
1254 }
1255
1256 if (!bFound) { /* oops, parent not found, error */
1257 printf("\n"
1258 "Error: parent in position %ld of %s must be\n"
1259 " declared before it when creating\n"
1260 " sampling dependencies - Exiting.\n\n",
1261 l, pMCVar1->pszName);
1262 exit(0);
1263 }
1264
1265 }
1266 } /* for l */
1267 } /* for k */
1268
1269 } /* FindMCParents */
1270
1271
1272 /* Function -------------------------------------------------------------------
1273 GetNumberOfMCVars
1274
1275 Find the total number of MC vars
1276 Called from TraverseLevels
1277 */
1278 void GetNumberOfMCVars (PLEVEL plevel, char **args)
1279 {
1280 long *pnThetas = (long *) args[0];
1281
1282 *pnThetas += plevel->nMCVars;
1283
1284 } /* GetNumberOfMCVars */
1285
1286
1287 /* Function -------------------------------------------------------------------
1288 InitMCVars
1289
1290 Sample initial values of thetas if not fixed
1291 Called from TraverseLevels
1292 */
1293 void InitMCVars(PLEVEL plevel, char **args)
1294 {
1295 long n;
1296
1297 for (n = 0; n < plevel->nMCVars; n++)
1298 if ( !(plevel->rgpMCVars[n]->bIsFixed))
1299 CalculateOneMCParm (plevel->rgpMCVars[n]);
1300
1301 } /* InitMCVars */
1302
1303
1304 /* Function -------------------------------------------------------------------
1305 ListToPMCArray
1306
1307 Convert list to array
1308 */
1309 void ListToPMCArray (PANALYSIS panal, PLIST plist,
1310 long *pnMCVars, PMCVAR **rgpMCVars)
1311 {
1312 if ((*pnMCVars = ListLength(plist)) == 0)
1313 return;
1314
1315 if (!(*rgpMCVars = (PMCVAR*) malloc (*pnMCVars * sizeof(PMCVAR))))
1316 ReportRunTimeError (panal, RE_OUTOFMEM | RE_FATAL, "ListToPMCArray");
1317
1318 *pnMCVars = 0;
1319 ForAllList3 (plist, &ListToPMCArrayL, pnMCVars, *rgpMCVars, NULL);
1320
1321 } /*ListToPMCArray */
1322
1323
1324 /* Function -------------------------------------------------------------------
1325 ListToPMCArrayL
1326 */
1327 void ListToPMCArrayL (PVOID pData, PVOID pUser1, PVOID pUser2, PVOID pUser3)
1328 {
1329 PMCVAR pMCVar = (PMCVAR)pData;
1330 long *pnMCVars = (long*)pUser1;
1331 PMCVAR *rgpMCVars = (PMCVAR*)pUser2;
1332
1333 rgpMCVars[(*pnMCVars)++] = pMCVar;
1334
1335 } /* ListToPMCArrayL */
1336
1337
1338 /* Function -------------------------------------------------------------------
1339 ListToPVArray
1340 */
1341 void ListToPVArray (PANALYSIS panal, PLIST plist,
1342 long *pnFixedVars, PVARMOD **rgpFixedVars)
1343 {
1344 if ((*pnFixedVars = ListLength (plist)) == 0)
1345 return;
1346
1347 if (!(*rgpFixedVars = (PVARMOD*) malloc (*pnFixedVars * sizeof(PVARMOD))))
1348 ReportRunTimeError (panal, RE_OUTOFMEM | RE_FATAL, "ListToPVArray");
1349
1350 *pnFixedVars = 0;
1351 ForAllList3 (plist, &ListToPVArrayL, pnFixedVars, *rgpFixedVars, NULL);
1352
1353 } /*ListToPVArray */
1354
1355
1356 /* Function -------------------------------------------------------------------
1357 ListToPVArrayL
1358 */
1359 void ListToPVArrayL (PVOID pData, PVOID pUser1, PVOID pUser2, PVOID pUser3)
1360 {
1361 PVARMOD pVar = (PVARMOD)pData;
1362 long *pnFixedVars = (long*)pUser1;
1363 PVARMOD *rgpFixedVars = (PVARMOD*)pUser2;
1364
1365 rgpFixedVars[(*pnFixedVars)++] = pVar;
1366
1367 } /* ListToPVArrayL */
1368
1369
1370 /* Function -------------------------------------------------------------------
1371 LnDensity
1372
1373 Returns the log of the (exact) density of variate under its distribution.
1374 */
1375 #define LN2PI 1.837877066409345339082
1376 #define LNSQRT2PI 0.918938533204672669541
1377 #define LNINVERPI -1.144729885849400163877
1378 #define LN2OVERPI -0.4515827052894548221396
1379
1380 double LnDensity (PMCVAR pMCVar, PANALYSIS panal)
1381 {
1382 double dTmp = 0, density;
1383 double dTmp2, dTmp3, dTmp4;
1384 double dParm1 = *(pMCVar->pdParm[0]);
1385 double dParm2 = *(pMCVar->pdParm[1]);
1386 double dMin = *(pMCVar->pdParm[2]);
1387 double dMax = *(pMCVar->pdParm[3]);
1388 double dTheta = pMCVar->dVal;
1389 char str[10];
1390
1391 /* This should take care of all dTheta checking */
1392 if (pMCVar->iType == MCV_BINOMIALBETA) {
1393 if (dTheta < 0) {
1394 printf("Error: variate out of bounds in LnDensity\n");
1395 exit(0);
1396 }
1397 }
1398 else if (pMCVar->iType == MCV_GENLOGNORMAL ||
1399 pMCVar->iType == MCV_STUDENTT) {
1400 if (dParm1 < 0) {
1401 printf("Error: parameter %g out of bounds in LnDensity\n",dParm1);
1402 exit(0);
1403 }
1404 }
1405 else {
1406 if (dTheta > dMax || dTheta < dMin) {
1407 return (NULL_SUPPORT);
1408 }
1409 }
1410
1411 switch (pMCVar->iType) {
1412
1413 case MCV_UNIFORM:
1414 if (dTheta > dParm2 || dTheta < dParm1)
1415 return (NULL_SUPPORT);
1416 if (dParm2 <= dParm1)
1417 ReportRunTimeError (panal, RE_BADUNIFORMDIST | RE_FATAL,
1418 pMCVar->pszName, "LnDensity");
1419 return -log(dParm2 - dParm1);
1420
1421 case MCV_LOGUNIFORM:
1422 if (dTheta > dParm2 || dTheta < dParm1)
1423 return (NULL_SUPPORT);
1424 if (dParm2 <= dParm1)
1425 ReportRunTimeError (panal, RE_BADUNIFORMDIST | RE_FATAL,
1426 pMCVar->pszName, "LnDensity");
1427 return -log(dTheta * (dParm2 - dParm1));
1428
1429 case MCV_NORMALV: dParm2 = sqrt (dParm2);
1430 return lnDFNormal (dTheta, dParm1, dParm2);
1431
1432 case MCV_NORMALCV: dParm2 = fabs(dParm1 * dParm2);
1433 return lnDFNormal (dTheta, dParm1, dParm2);
1434
1435 case MCV_NORMAL:
1436 case MCV_HALFNORMAL:
1437 return lnDFNormal (dTheta, dParm1, dParm2);
1438
1439 case MCV_LOGNORMALV: dParm2 = exp(sqrt(dParm2)); /* fall thru */
1440 case MCV_LOGNORMAL:
1441 if (dParm1 <= 0.0) {
1442 sprintf(str, "%5.2e", dParm1);
1443 ReportRunTimeError(panal, RE_BADLOGNORMALMEAN | RE_FATAL,
1444 pMCVar->pszName, str, "LnDensity");
1445 }
1446 return (lnDFNormal(log(dTheta),log(dParm1),log(dParm2)) - log(dTheta));
1447
1448 case MCV_TRUNCNORMALV: dParm2 = sqrt (dParm2);
1449 return lnDFNormal (dTheta, dParm1, dParm2) -
1450 log(CDFNormal ((dMax - dParm1) / dParm2) -
1451 CDFNormal ((dMin - dParm1) / dParm2));
1452
1453 case MCV_TRUNCNORMALCV: dParm2 = fabs(dParm1 * dParm2);
1454 return lnDFNormal (dTheta, dParm1, dParm2) -
1455 log(CDFNormal ((dMax - dParm1) / dParm2) -
1456 CDFNormal ((dMin - dParm1) / dParm2));
1457
1458 case MCV_TRUNCNORMAL:
1459 if (dParm2 <= 0.0) {
1460 sprintf(str, "%5.2e", dParm2);
1461 ReportRunTimeError(panal, RE_BADNORMALSD | RE_FATAL,
1462 pMCVar->pszName, str, "LnDensity");
1463 }
1464 return lnDFNormal (dTheta, dParm1, dParm2) -
1465 log(CDFNormal ((dMax - dParm1) / dParm2) -
1466 CDFNormal ((dMin - dParm1) / dParm2));
1467
1468 case MCV_TRUNCLOGNORMALV: dParm2 = exp(sqrt(dParm2)); /* fall thru */
1469 case MCV_TRUNCLOGNORMAL:
1470 if (dParm1 <= 0.0 ) {
1471 sprintf(str, "%5.2e", dParm1);
1472 ReportRunTimeError(panal, RE_BADLOGNORMALMEAN | RE_FATAL,
1473 pMCVar->pszName, str, "LnDensity");
1474 }
1475 if (dParm2 <= 1.0 ) {
1476 sprintf(str, "%5.2e", dParm2);
1477 ReportRunTimeError(panal, RE_BADLOGNORMALSD | RE_FATAL,
1478 pMCVar->pszName, str, "LnDensity");
1479 }
1480 dTmp = log(dParm2);
1481 return lnDFNormal (log(dTheta), log(dParm1), dTmp) - log(dTheta) -
1482 log(CDFNormal (log(dMax / dParm1) / dTmp) -
1483 CDFNormal (log(dMin / dParm1) / dTmp));
1484
1485 case MCV_BETA:
1486 return lnDFBeta (dTheta, dParm1, dParm2, dMin, dMax);
1487
1488 case MCV_CHI2:
1489 dTmp = 0.5 * dParm1;
1490 return (dTmp - 1) * log(dTheta) - 0.5 * dTheta +
1491 dTmp * (-6.9314718056E-01) - lnGamma (dTmp);
1492
1493 case MCV_BINOMIAL:
1494 if ((dParm1 < 0) || (dParm1 > 1)) {
1495 printf("Error: bad p for binomial variate in LnDensity\n");
1496 exit (0);
1497 }
1498 if (dTheta > dParm2) {
1499 return (NULL_SUPPORT);
1500 }
1501 /* log binomial coefficient n! / (x!(n-x)!) */
1502 dTmp = lnGamma (dParm2 + 1) - lnGamma (dTheta + 1) -
1503 lnGamma (dParm2 - dTheta + 1);
1504
1505 if (dParm1 == 0) { /* revised FB 18/07/97 */
1506 if (dTheta != 0)
1507 return (NULL_SUPPORT); /* should be -INF */
1508 }
1509 else
1510 dTmp = dTmp + dTheta * log(dParm1);
1511
1512 if (dParm1 == 1) { /* revised FB 18/07/97 */
1513 if ((dParm2 - dTheta) == 0)
1514 return dTmp; /* because log(0^0) is 0 */
1515 else
1516 return (NULL_SUPPORT); /* should be -INF */
1517 }
1518 else
1519 return dTmp + (dParm2 - dTheta) * log(1 - dParm1);
1520
1521 case MCV_NEGATIVEBINOM:
1522 /* dParm1: r, number of runs before success,
1523 dParm2: p, probability of success */
1524 if ((dParm2 < 0) || (dParm2 > 1)) {
1525 printf("Error: bad p for negative binomial variate in LnDensity\n");
1526 exit (0);
1527 }
1528 /* log binomial coefficient ((x + r - 1)! / (x!(r - 1)!) */
1529 dTmp = lnGamma (dTheta + dParm1) - lnGamma (dTheta + 1) -
1530 lnGamma (dParm1);
1531
1532 if ((dParm2 == 0) && (dTheta != 0))
1533 return (NULL_SUPPORT); /* should be -INF */
1534 else
1535 dTmp = dTmp + dTheta * log(dParm2);
1536
1537 if (dParm2 == 1) {
1538 if (dParm1 == 0)
1539 return dTmp; /* because log(0^0) is 0 */
1540 else
1541 return (NULL_SUPPORT); /* should be -INF */
1542 }
1543 else
1544 return (dTmp + dParm1 * log(1 - dParm2));
1545
1546 case MCV_PIECEWISE:
1547 density = 2 / (dMax + dParm2 - dParm1 - dMin);
1548
1549 if (dTheta <= dParm1)
1550 return log(density * (dTheta - dMin) / (dParm1 - dMin));
1551
1552 else
1553 if (dTheta <= dParm2)
1554 return log(density);
1555 else
1556 return log(density * (dMax - dTheta) / (dMax - dParm2));
1557
1558 case MCV_EXPONENTIAL:
1559 if (dParm1 <= 0) {
1560 printf ("Error: negative or null inverse scale (%g) for exponential "
1561 "variate in LnDensity\n", dParm1);
1562 exit (0);
1563 }
1564 return -dTheta * dParm1 + log(dParm1);
1565
1566 case MCV_GGAMMA:
1567 if (dParm2 <= 0) {
1568 printf ("Error: bad inv. scale for gamma variate in LnDensity\n");
1569 exit (0);
1570 }
1571 return (dParm1 - 1) * log(dTheta) - dParm2 * dTheta +
1572 dParm1 * log(dParm2) - lnGamma (dParm1);
1573
1574 case MCV_TRUNCINVGGAMMA:
1575 #ifdef HAVE_LIBGSL
1576 dTmp = gsl_cdf_gamma_P (1/dMax, dParm1, dParm2) -
1577 gsl_cdf_gamma_P (1/dMin, dParm1, dParm2); /* fall thru */
1578 #else
1579 printf("Error: Truncated inverse gamma density cannot be evaluated\n");
1580 printf(" if the GNU Scientific Library is not installed\n");
1581 exit(0);
1582 #endif
1583 case MCV_INVGGAMMA:
1584 if (dParm2 <= 0) {
1585 printf("Error: bad scale for inv. gamma variate in LnDensity\n");
1586 exit(0);
1587 }
1588 return (-dParm1 - 1) * log(dTheta) - dParm2 / dTheta +
1589 dParm1 * log(dParm2) - lnGamma (dParm1) - dTmp;
1590
1591 case MCV_POISSON:
1592 if (dParm1 <= 0) {
1593 printf("Error: bad rate for Poisson variate in LnDensity\n");
1594 exit(0);
1595 }
1596 return dTheta * log(dParm1) - dParm1 - lnGamma (dTheta + 1);
1597
1598 case MCV_BINOMIALBETA:
1599 if (dParm1 < 0) {
1600 printf("Error: bad expectation for BinomialBeta variate "
1601 "in LnDensity\n");
1602 exit(0);
1603 }
1604 if (dParm2 <= 0) {
1605 printf("Error: bad alpha for BinomialBeta variate in LnDensity\n");
1606 exit(0);
1607 }
1608 if (dMin <= 0) {
1609 printf("Error: bad beta for BinomialBeta variate in LnDensity\n");
1610 exit(0);
1611 }
1612 dTmp = floor (0.5 + dParm1 + dParm1 * dMin / dParm2); /* this is N */
1613 if (dTheta > dTmp)
1614 return (NULL_SUPPORT);
1615 else {
1616 if ((dParm2 == 0.5) && (dMin == 0.5))
1617 dTmp = lnGamma (0.5 + dTheta) +
1618 lnGamma (0.5 + dTmp - dTheta) -
1619 lnGamma (dTheta + 1) - lnGamma (dTmp - dTheta + 1);
1620 else
1621 dTmp = lnGamma (dParm2 + dMin) + lnGamma (dTmp + 1) +
1622 lnGamma (dParm2 + dTheta) +
1623 lnGamma (dMin + dTmp - dTheta) -
1624 lnGamma (dTheta + 1) - lnGamma (dTmp - dTheta + 1) -
1625 lnGamma (dParm2) - lnGamma(dMin) -
1626 lnGamma (dParm2 + dMin + dTmp);
1627 return dTmp;
1628 }
1629
1630 case MCV_GENLOGNORMAL: /* Generalized LogNormal */
1631 if (dParm1 < 0) {
1632 printf("Error: bad expectation for GenLogNormal variate "
1633 "in LnDensity\n");
1634 exit(0);
1635 }
1636 /* This is relative stdev of lognormal part -- dMin is sigma for the
1637 lognormal part */
1638 dTmp = sqrt(exp(pow(dMin,2)) * (exp(pow(dMin,2)) - 1));
1639 /* This is the transformation parameter lambda */
1640 dTmp2 = pow(dParm2/dTmp,2);
1641 /* This is the transformed mean */
1642 dTmp3 = log(dParm1 + sqrt(pow(dParm1,2) + dTmp2));
1643 /* This is the transformed theta -- use Taylor series for
1644 negative dTheta when lambda/dTheta^2 < 0.01 */
1645 if ((dTheta < 0) && (dTmp2 < (0.01*dTheta*dTheta)))
1646 dTmp4 = log(dTmp2/(-2*dTheta)*(1+0.25*dTmp2/pow(dTheta,2)));
1647 else
1648 dTmp4 = log(dTheta + sqrt(pow(dTheta,2) + dTmp2));
1649 /* Normal log-density minus log-jacobian */
1650 return lnDFNormal( dTmp4, dTmp3, dTmp ) - 0.5*log(pow(dTmp4,2) + dTmp2);
1651
1652 case MCV_STUDENTT: /* Student t,
1653 dParm1 is dof, dParm2 is m, dMin is sigma */
1654 if (dParm1 <= 0) {
1655 printf("Error: bad dof for Student-T variate"
1656 "in LnDensity\n");
1657 exit(0);
1658 }
1659 dTmp = (dParm1 + 1)/ 2;
1660 return (lnGamma(dTmp) - lnGamma(dParm1 / 2)
1661 -0.5 * log(dParm1 * PI *dMin *dMin)
1662 -dTmp * log(1 + pow((dTheta - dParm2) / dMin, 2) / dParm1));
1663
1664 case MCV_CAUCHY:
1665 return (LNINVERPI - log(dParm1 + dTheta * dTheta / dParm1));
1666
1667 case MCV_HALFCAUCHY: /* dParm1 is scale */
1668 return (LN2OVERPI - log(dParm1 + dTheta * dTheta / dParm1));
1669
1670 case MCV_USERLL: /* dParm1 is the model computed log-likelihood */
1671 return (dParm1);
1672
1673 default:
1674 ReportRunTimeError(panal, RE_UNKNOWNDIST | RE_FATAL, "LnDensity");
1675
1676 } /* switch */
1677
1678 /* Not reached */
1679 return 0.0 ;
1680
1681 } /* LnDensity */
1682
1683
1684 /* Function -------------------------------------------------------------------
1685 LnLike
1686
1687 returns the log-likelihood of the dependents given the parameter specified
1688 by pMCVar
1689 */
1690 double LnLike (PMCVAR pMCVar, PANALYSIS panal)
1691 {
1692 long n;
1693 double dDensity, dLnLike = 0.0;
1694
1695 for (n = 0; n < pMCVar->nDependents; n++) {
1696 dDensity = LnDensity(pMCVar->rgpDependents[n], panal);
1697 if (dDensity == NULL_SUPPORT)
1698 return NULL_SUPPORT;
1699 else
1700 dLnLike += dDensity;
1701 }
1702 return dLnLike;
1703
1704 } /* LnLike */
1705
1706
1707 /* Function -------------------------------------------------------------------
1708 LnLikeData
1709
1710 Likelihood of the data for one experiment
1711 */
1712 double LnLikeData (PLEVEL plevel, PANALYSIS panal)
1713 {
1714 PMCVAR pMCVar;
1715 long i, j, k;
1716 double dLnLike = 0.0;
1717 double dTmp;
1718 BOOL bMissData, bMissOutp;
1719 static PDOUBLE pdBase[4];
1720
1721 /* For all the likelihoods seen at the experiment level */
1722 for (i = 0; i < plevel->nLikes; i++) {
1723 pMCVar = plevel->rgpLikes[i];
1724
1725 for (k = 0; k < 4; k++)
1726 pdBase[k] = pMCVar->pdParm[k]; /* store the base pointers of pdParm */
1727
1728 /* Set dVal of pMCVar to the current data value */
1729 for (j = 0; j < pMCVar->lCount; j++) { /* for each data value */
1730
1731 pMCVar->dVal = pMCVar->pdVal[j]; /* point to jth data value */
1732
1733 /* If the main data is coded as missing do nothing, and otherwise: */
1734 if (pMCVar->dVal != INPUT_MISSING_VALUE) {
1735 /* Set the pdParms of pMCVar to point to data or output values */
1736 bMissData = FALSE; /* initialize */
1737 bMissOutp = FALSE; /* initialize */
1738 for (k = 0; k < 4; k++) {
1739 if (pMCVar->iParmType[k] == MCVP_PRED) {
1740 /* Advance in the prediction array */
1741 pMCVar->pdParm[k] = pdBase[k] + j;
1742 bMissOutp = bMissOutp + (*(pMCVar->pdParm[k]) == MISSING_VALUE);
1743 }
1744 else if (pMCVar->iParmType[k] == MCVP_DATA) {
1745 /* Advance in the data array */
1746 pMCVar->pdParm[k] = pdBase[k] + j;
1747 bMissData = bMissData +
1748 (*(pMCVar->pdParm[k]) == INPUT_MISSING_VALUE);
1749 }
1750 } /* end for k */
1751
1752 /* If no missing data among the k parameters */
1753 if (bMissData == FALSE) {
1754 /* If no missing model output among the k parameters */
1755 if (bMissOutp == FALSE) {
1756 dTmp = LnDensity (pMCVar, panal);
1757 if (dTmp == NULL_SUPPORT) {
1758 /* Reset the pdParms to the beginning of arrays FB 9/6/99 */
1759 for (k = 0; k < 4; k++)
1760 pMCVar->pdParm[k] = pdBase[k];
1761 return (NULL_SUPPORT);
1762 }
1763 else
1764 dLnLike = dLnLike + dTmp;
1765 }
1766 else /* missing output, cannot compute */
1767 ReportRunTimeError (panal, RE_BADMODEL | RE_FATAL, "LnLikeData");
1768 } /* end if bMissData */
1769
1770 } /* end if ! INPUT_MISSING_VALUE */
1771 } /* end for j */
1772
1773 /* Reset the pdParms to the beginning of the output or data arrays */
1774 for (k = 0; k < 4; k++)
1775 pMCVar->pdParm[k] = pdBase[k];
1776
1777 } /* end for i */
1778
1779 return (dLnLike);
1780
1781 } /* LnLikeData */
1782
1783
1784 /* Function -------------------------------------------------------------------
1785 MaxMCVar
1786
1787 Get the upper bound of an MCVar (random variable).
1788 */
1789 double MaxMCVar (PMCVAR pMCVar)
1790 {
1791 /* if the parameter has a discrete distribution, round it - FB 12/06/97 */
1792 if (pMCVar->iType == MCV_BINOMIAL || pMCVar->iType == MCV_POISSON ) {
1793 return( *(pMCVar->pdParm[3]));
1794 }
1795 else { /* FB fixed the uniform case - FB 30/06/97 */
1796 if (pMCVar->iType == MCV_UNIFORM || pMCVar->iType == MCV_LOGUNIFORM ) {
1797 return( *(pMCVar->pdParm[1]));
1798 }
1799 else {
1800 return( *(pMCVar->pdParm[3]));
1801 }
1802 }
1803
1804 } /* MaxMCVar*/
1805
1806
1807 /* Function -------------------------------------------------------------------
1808 MinMCVar
1809
1810 Get the lower bound of an MCVar (a random variable)
1811 */
1812 double MinMCVar (PMCVAR pMCVar)
1813 {
1814 /* if the parameter has a discrete distribution, round it - FB 12/06/97 */
1815 if (pMCVar->iType == MCV_BINOMIAL || pMCVar->iType == MCV_POISSON) {
1816 return(*(pMCVar->pdParm[2]));
1817 }
1818 else { /* FB fixed the uniform case - FB 30/06/97 */
1819 if (pMCVar->iType == MCV_UNIFORM || pMCVar->iType == MCV_LOGUNIFORM) {
1820 return( *(pMCVar->pdParm[0]));
1821 }
1822 else {
1823 return( *(pMCVar->pdParm[2]));
1824 }
1825 }
1826
1827 } /* MinMCVar */
1828
1829
1830 /* Function -------------------------------------------------------------------
1831 RunTemperingBlock
1832
1833 Run a block of MCMC simulations to adjust the perk scale and pseudo-priors
1834 */
1835 void RunTemperingBlock (PANALYSIS panal, long lRunLength, PLONG iter)
1836 {
1837 PGIBBSDATA pgd = &panal->gd;
1838 PLEVEL pLevel0 = panal->pLevels[0];
1839 double dTmp, dLnPrior = 0, dLnData = 0;
1840 long i, j;
1841 long nUpdateAt, nTotal;
1842
1843 for (i = 0; i < lRunLength; i++) { /* run block */
1844
1845 nTotal = UPDATE_BASE;
1846 nUpdateAt = *iter + nTotal;
1847
1848 TraverseLevels (pLevel0, SampleThetasTempered, panal, pgd, &i,
1849 &nUpdateAt, &nTotal, &pgd->indexT, NULL);
1850
1851 dLnPrior = 0.0;
1852 TraverseLevels (pLevel0, CalculateTotals, panal, &dLnPrior, NULL);
1853 dLnData = 0.0;
1854 TraverseLevels1 (pLevel0, SumAllExpts, &dLnData, NULL);
1855
1856 /* special? Robbins-Munro updating of the pseudo-priors */
1857 dTmp = pgd->dCZero / (double) (i + pgd->dNZero);
1858 /*if (pgd->indexT == pgd->startT) {
1859 pgd->rgdlnPi[pgd->startT] -= dTmp;
1860 for (j = pgd->startT + 1; j <= pgd->endT; j++)
1861 pgd->rgdlnPi[j] += dTmp / (double) pgd->nPerks;
1862 }
1863 else {
1864 pgd->rgdlnPi[pgd->startT] += dTmp;
1865 for (j = pgd->startT + 1; j <= pgd->endT; j++)
1866 pgd->rgdlnPi[j] -= dTmp / (double) pgd->nPerks;
1867 }
1868 */
1869
1870 for (j = pgd->startT; j <= pgd->endT; j++) {
1871 if (j == pgd->indexT)
1872 pgd->rgdlnPi[j] -= dTmp;
1873 else
1874 pgd->rgdlnPi[j] += dTmp / (double) pgd->nPerks;
1875 }
1876
1877 /* update population count of current temperature */
1878 pgd->rglCount[pgd->indexT] = pgd->rglCount[pgd->indexT]+1;
1879
1880 /* test the temperature and change indexT if necessary */
1881 pgd->indexT = SampleTemperature2 (pgd, dLnPrior, dLnData);
1882
1883 /* Adjust the update time eventually */
1884 if (i == nUpdateAt) {
1885 nTotal = (nTotal * 3) / 2;
1886 nUpdateAt = i + nTotal;
1887 }
1888
1889 /* Increment the total iteration counter */
1890 (*iter)++;
1891
1892 } /* end for */
1893
1894 } /* RunTemperingBlock */
1895
1896
1897 /* Function -------------------------------------------------------------------
1898 NextDown
1899
1900 Return from a table the element just inferior to the argument
1901 */
1902 double NextDown (double Perk)
1903 {
1904 int i;
1905 static double PTable[21] = {0, 1E-6, 1E-5, 1E-4, 1E-3, 1E-2,
1906 0.1, 0.2, 0.3, 0.5, 0.6, 0.7, 0.8, 0.9,
1907 0.95, 0.97, 0.99, 0.999, 0.9999, 0.99999, 1};
1908
1909 i = 0;
1910 while (Perk > PTable[i]) {
1911 i++;
1912 }
1913
1914 return (i == 0 ? PTable[i] : PTable[i-1]);
1915
1916 } /* NextDown */
1917
1918
1919 /* Function -------------------------------------------------------------------
1920 InitPerks
1921
1922 Adjusts automatically the perk schedule and pseudo-priors
1923 for tempered MCMC, if the user has not specified values.
1924 */
1925 void InitPerks (PANALYSIS panal)
1926 {
1927 PGIBBSDATA pgd = &panal->gd;
1928 long i, j, k, iter = 0;
1929 double dTmp;
1930 int bTrans;
1931 BOOL bHappy, bTooManyTrials;
1932
1933 /* if we are doing tempered MCMC */
1934 if ((pgd->nSimTypeFlag == 3) || (pgd->nSimTypeFlag == 4)) {
1935
1936 /* allocate transition counts array */
1937 if (!(pgd->rglTransAttempts = InitlVector (NTEMP)) ||
1938 !(pgd->rglTransAccepts = InitlVector (NTEMP)))
1939 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "InitPerks", NULL);
1940
1941 /* initialize them at zero */
1942 for (i = 0; i < NTEMP; i++)
1943 pgd->rglTransAttempts[i] = pgd->rglTransAccepts[i] = 0;
1944
1945 /* if perks not set by the user */
1946 if (pgd->nPerks == 0) {
1947
1948 /* Screen message */
1949 printf ("Setting perks (inverse temperatures).\n");
1950
1951 pgd->nPerks = NTEMP;
1952
1953 /* allocate perk, pseudo-prior and population count arrays */
1954 if (!(pgd->rgdPerks = InitdVector (NTEMP)) ||
1955 !(pgd->rgdlnPi = InitdVector (NTEMP)) ||
1956 !(pgd->rglCount = InitlVector (NTEMP)))
1957 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "InitPerks", NULL);
1958
1959 for (i = 0; i < NTEMP; i++) {
1960 pgd->rgdPerks[i] = pgd->rgdlnPi[i] = pgd->rglCount[i] = 0;
1961 }
1962
1963 /* start at perk 1 and just below it */
1964 pgd->endT = NTEMP - 1;
1965 pgd->startT = NTEMP - 2;
1966 pgd->indexT = pgd->startT;
1967
1968 double dEPSILON = 0.99;
1969 double dUP = 2.0;
1970 pgd->rgdPerks[pgd->startT] = dEPSILON;
1971 pgd->rgdPerks[pgd->endT] = 1.00;
1972
1973 long nOldPrintIter = pgd->nPrintIter;
1974 pgd->nPrintIter = -pgd->nMaxPerkSetIter;
1975 int lRunLength = 100;
1976 double dBoundary = 0.0;
1977
1978 do { /* loop over running a block, checking results, and adjusting */
1979
1980 pgd->indexT = pgd->startT;
1981
1982 /* run a batch of tempered MCMC simulations */
1983 RunTemperingBlock (panal, lRunLength, &iter);
1984
1985 /* post-run diagnostic printing */
1986 PrintTemperatureDiagnostics (stdout, pgd);
1987 PrintTemperatureDiagnostics (pgd->pfilePerks, pgd);
1988
1989 /* check block */
1990 bTrans = CheckTransitions (pgd);
1991 if (pgd->rgdPerks[pgd->startT] == dBoundary) {
1992 /* boundary perk reached */
1993 bHappy = (bTrans > -1); /* transition OK or too high */
1994 }
1995 else { /* boundary perk not reached */
1996 bHappy = FALSE;
1997 }
1998 bTooManyTrials = (iter > pgd->nMaxPerkSetIter);
1999
2000 /* if unhappy: adjust */
2001 if (!bHappy) {
2002
2003 if (bTrans == -1) { /* acceptance rate too low */
2004
2005 printf ("acceptance rate 1<->2 too low, stepping back up\n");
2006
2007 /* go back up half way to the point above */
2008 dTmp = (pgd->rgdPerks[pgd->startT] +
2009 pgd->rgdPerks[pgd->startT+1]) / dUP;
2010
2011 /* adjust pseudo-prior with an average of old and extrapolated */
2012 pgd->rgdlnPi[pgd->startT] = (pgd->rgdlnPi[pgd->startT] +
2013 Extrapolate (pgd, dTmp, pgd->startT, pgd->startT+1)) / 2;
2014
2015 pgd->rgdPerks[pgd->startT] = dTmp;
2016
2017 } /* end too low */
2018
2019 if (bTrans == +1) { /* acceptance rate too high */
2020
2021 printf ("acceptance rate 1<->2 too high, moving down\n");
2022
2023 dTmp = NextDown(pgd->rgdPerks[pgd->startT]);
2024
2025 /* adjust pseudo-prior with an average of old and extrapolated */
2026 pgd->rgdlnPi[pgd->startT] = (pgd->rgdlnPi[pgd->startT] +
2027 Extrapolate (pgd, dTmp, pgd->startT, pgd->startT+1)) / 2;
2028
2029 pgd->rgdPerks[pgd->startT] = dTmp;
2030
2031 } /* end too high */
2032
2033 if (bTrans == 0) {
2034 /* acceptance rate ok, but target perk not reached
2035 check if we have space to add a point below */
2036 if (pgd->startT > 0) { /* we do have space */
2037
2038 printf ("acceptance rate 1<->2 ok, adding a new point\n");
2039
2040 /* if all transitions rates are OK the pseudo-priors should be
2041 about right and we might be able to take some intermediate
2042 points off the scale */
2043 if (CheckAllTransitions(pgd)) {
2044 /* remove recursively center point if 3 points are aligned */
2045 int i = pgd->endT;
2046 int j = i - 2;
2047 while (j >= pgd->startT) {
2048 if (EqualSlopes(pgd->rgdPerks, pgd->rgdlnPi, j)) {
2049 /* remove point i - 1 from the scale, repacking
2050 the scale that implies copying perks,
2051 pseudo-prior, and counts from positions start
2052 to (i - 2) to positions (start + 1) to (i - 1)
2053 and setting start to (start + 1). start has
2054 moved up, so j does not need to be updated */
2055 for (k = j; k >= pgd->startT; k--) {
2056 pgd->rgdPerks[k+1] = pgd->rgdPerks[k];
2057 pgd->rgdlnPi [k+1] = pgd->rgdlnPi [k];
2058 pgd->rglCount[k+1] = pgd->rglCount[k];
2059 }
2060 pgd->startT++;
2061 if (pgd->indexT <= j)
2062 pgd->indexT++;
2063 lRunLength = lRunLength - 100;
2064 printf("Scale has been reduced.\n");
2065 }
2066 else { /* slopes not equal, move i and j down */
2067 i--;
2068 j--;
2069 }
2070 }
2071 }
2072
2073 pgd->startT = pgd->startT - 1;
2074 pgd->indexT = pgd->startT;
2075
2076 /* give perk to the new point */
2077 pgd->rgdPerks[pgd->startT] =
2078 NextDown(pgd->rgdPerks[pgd->startT+1]);
2079
2080 /* pseudo-prior of new point is an average of above and
2081 extrapolated */
2082 pgd->rgdlnPi[pgd->startT] = (pgd->rgdlnPi[pgd->startT+1] +
2083 Extrapolate (pgd, pgd->rgdPerks[pgd->startT],
2084 pgd->startT+1, pgd->startT+2)) / 2;
2085
2086 lRunLength = lRunLength + 100;
2087 }
2088 else /* no more space in scale, stop */
2089 bTooManyTrials = TRUE;
2090
2091 } /* end OK */
2092
2093 /* re-initialize counts at zero */
2094 for (i = pgd->startT; i <= pgd->endT; i++) {
2095 pgd->rglCount[i] = 0;
2096 pgd->rglTransAttempts[i] = pgd->rglTransAccepts[i] = 0;
2097 }
2098
2099 } /* end if !bHappy */
2100
2101 } while ((!bHappy) && (!bTooManyTrials));
2102
2103 if (pgd->rgdPerks[pgd->startT] == dBoundary)
2104 printf ("Perk %lg reached in %ld iterations.\n", dBoundary, iter);
2105 else
2106 printf ("Perk %lg not reached in %ld iterations...\n", dBoundary, iter);
2107
2108 /* restore */
2109 pgd->nPrintIter = nOldPrintIter;
2110
2111 /* prepare the actual runs */
2112 int iCount = pgd->endT - pgd->startT + 1;
2113 if (iCount != NTEMP) { /* shift array contents to start at 0 */
2114 pgd->nPerks = iCount;
2115 pgd->indexT = pgd->indexT - pgd->startT;
2116 for (i = 0; i < iCount; i++) {
2117 j = pgd->startT + i;
2118 pgd->rgdPerks[i] = pgd->rgdPerks[j];
2119 pgd->rgdlnPi[i] = pgd->rgdlnPi[j];
2120 pgd->rglCount[i] = 0;
2121 }
2122 pgd->startT = 0;
2123 pgd->endT = iCount - 1;
2124 }
2125
2126 printf ("Done with InitPerks - Continuing.\n\n");
2127
2128 } /* if perks not set by user */
2129 } /* if tempered */
2130
2131 } /* InitPerks */
2132
2133
2134 /* Function -------------------------------------------------------------------
2135 OpenMarkovFiles
2136
2137 Opens the output file and the restart file.
2138 */
2139 void OpenMarkovFiles (PANALYSIS panal)
2140 {
2141 PGIBBSDATA pgd = &panal->gd;
2142 char* with_rank;
2143
2144 /* if we are debugging the structure, do nothing here */
2145 if (panal->bDependents) return;
2146
2147 /* Take care of the output file first */
2148
2149 /* Use command line spec if given */
2150 if (panal->bCommandLineSpec) {
2151 free (pgd->szGout);
2152 panal->bAllocatedFileName = FALSE;
2153 pgd->szGout = panal->szOutfilename;
2154 }
2155
2156 /* Default if none given */
2157 else if (!(pgd->szGout))
2158 pgd->szGout = "MCMC.default.out";
2159
2160 if (panal->size > 1) {
2161 with_rank = malloc(sizeof(char)*(6+strlen(pgd->szGout)));
2162 sprintf(with_rank,"%04d_%s",panal->rank,pgd->szGout);
2163 pgd->szGout = with_rank;
2164 }
2165
2166 /* Eventually open the restart file before crushing the output file */
2167 if (pgd->szGrestart)
2168 if (!(pgd->pfileRestart)
2169 && !(pgd->pfileRestart = fopen (pgd->szGrestart, "r")))
2170 ReportRunTimeError(panal, RE_FATAL | RE_CANNOTOPEN,
2171 pgd->szGrestart, "OpenMarkovFiles");
2172
2173 if (!(pgd->pfileOut)
2174 && !(pgd->pfileOut = fopen (pgd->szGout, "w")))
2175 ReportRunTimeError(panal, RE_FATAL | RE_CANNOTOPEN,
2176 pgd->szGout, "OpenMarkovFiles");
2177
2178 /* If tempered MCMC, open perk scale recording file */
2179 if ((pgd->nSimTypeFlag == 3) || (pgd->nSimTypeFlag == 4)) {
2180 char szFileName[MAX_FILENAMESIZE+6];
2181 sprintf(szFileName, "%s%s", pgd->szGout, ".perks");
2182 if (!(pgd->pfilePerks)
2183 && !(pgd->pfilePerks = fopen (szFileName, "w")))
2184 ReportRunTimeError(panal, RE_FATAL | RE_CANNOTOPEN,
2185 szFileName, "OpenMarkovFiles");
2186 }
2187
2188 } /* OpenMarkovFiles */
2189
2190
2191 /* Function -------------------------------------------------------------------
2192 PrintAllExpts
2193
2194 Run and print all experiments at level below `plevel'
2195 */
2196 void PrintAllExpts (PLEVEL plevel, PANALYSIS panal, PFILE pOutFile)
2197 {
2198 long n;
2199
2200 for (n = 0; n < plevel->iInstances; n++)
2201 TraverseLevels1 (plevel->pLevels[n], PrintExpt, panal, pOutFile, NULL);
2202
2203 } /* PrintAllExpts */
2204
2205
2206 /* Function -------------------------------------------------------------------
2207 PrintDeps
2208
2209 Called from TraverseLevels
2210
2211 For debugging, print the variables, parents, and dependencies
2212 */
2213 void PrintDeps (PLEVEL plevel, char **args)
2214 {
2215 long n, m;
2216 PMCVAR pMCVar;
2217
2218 printf("Depth %d; Instance %d\n", plevel->iDepth, plevel->iSequence);
2219
2220 for (n = 0; n < plevel->nMCVars; n++) {
2221 pMCVar = plevel->rgpMCVars[n];
2222
2223 printf ("Variable %s (%d) [%" PRIxPTR "]\n",
2224 pMCVar->pszName, pMCVar->iDepth, (intptr_t) pMCVar);
2225
2226 for (m = 0; m < 4; m++)
2227 if (pMCVar->pMCVParent[m] != NULL)
2228 printf (" Parent %ld: %s (%d) [%" PRIxPTR "]\n", m,
2229 pMCVar->pMCVParent[m]->pszName, pMCVar->pMCVParent[m]->iDepth,
2230 (intptr_t) pMCVar->pMCVParent[m]);
2231
2232 for (m = 0; m < pMCVar->nDependents; m++)
2233 printf (" Dependent: %s (%d) [%" PRIxPTR "]\n",
2234 pMCVar->rgpDependents[m]->pszName,
2235 pMCVar->rgpDependents[m]->iDepth,
2236 (intptr_t) pMCVar->rgpDependents[m]);
2237
2238 if (pMCVar->bExptIsDep)
2239 printf(" This variable influences experiments directly\n");
2240 }
2241
2242 } /* PrintDeps */
2243
2244
2245 /* Function -------------------------------------------------------------------
2246 PrintExpt
2247
2248 Run all experiments and print the level code, experiment number, time, data,
2249 and predictions to the output file.
2250
2251 Called from TraverseLevels1
2252 */
2253 int PrintExpt (PLEVEL plevel, char **args)
2254 {
2255 PANALYSIS panal = (PANALYSIS)args[0];
2256 PFILE pOutFile = (PFILE)args[1];
2257 long k, l, m, n;
2258 PEXPERIMENT pExpt = plevel->pexpt;
2259 POUTSPEC pos;
2260 static long printed_head = 0;
2261
2262 if (!printed_head) {
2263 fprintf (pOutFile,
2264 "Level\tSimulation\tOutput_Var\tTime\tData\tPrediction\n");
2265 printed_head = 1;
2266 }
2267
2268 /* Set level sequence */
2269 panal->pCurrentLevel[plevel->iDepth] = plevel;
2270
2271 panal->iInstance[plevel->iDepth] = plevel->iSequence;
2272
2273 if (pExpt != NULL) {
2274 InitModel ();
2275
2276 /* Set the model vars that have been sampled in this experiment and
2277 above levels */
2278 for (n = 0; n <= plevel->iDepth; n++) {
2279 SetModelVars (panal->pCurrentLevel[n]);
2280 SetFixedVars (panal->pCurrentLevel[n]);
2281 }
2282
2283 if (!DoOneExperiment (pExpt)) {
2284 /* Error */
2285 printf ("Warning: DoOneExperiment failed\n");
2286 return 0;
2287 }
2288 else {
2289 pos = &pExpt->os;
2290 for (m = 0; m < pos->nOutputs; m++) {
2291 /* find the corresponding data index in case Print and Data are not
2292 in the same order */
2293 for (k = 0; k < pos->nData; k++)
2294 if (!strcmp(pos->pszDataNames[k], pos->pszOutputNames[m]))
2295 break;
2296
2297 for (l = 0; l < pos->pcOutputTimes[m]; l++) {
2298 for (n = 1; n < plevel->iDepth; n++)
2299 fprintf (pOutFile, "%d_", panal->iInstance[n]);
2300 fprintf (pOutFile, "%d\t", panal->iInstance[plevel->iDepth]);
2301
2302 if (k != pos->nData) /* Data found */
2303 fprintf (pOutFile, "%d\t%s\t%g\t%g\t%g\n", pExpt->iExp,
2304 pos->pszOutputNames[m], pos->prgdOutputTimes[m][l],
2305 pos->prgdDataVals[k][l], pos->prgdOutputVals[m][l]);
2306 else /* data not found, empty field */
2307 fprintf (pOutFile, "%d\t%s\t%g\t\t%g\n", pExpt->iExp,
2308 pos->pszOutputNames[m], pos->prgdOutputTimes[m][l],
2309 pos->prgdOutputVals[m][l]);
2310 } /* end for l */
2311 fprintf (pOutFile, "\n");
2312
2313 } /* end for m */
2314 fprintf (pOutFile, "\n");
2315
2316 } /* end else */
2317
2318 } /* end if pExpt */
2319
2320 return (1);
2321
2322 } /* PrintExpt*/
2323
2324
2325 /* Function -------------------------------------------------------------------
2326 PrintTemperatureDiagnostics
2327
2328 For debugging and information, print the state of the tempering algorithm
2329 */
2330 void PrintTemperatureDiagnostics (PFILE fOut, PGIBBSDATA pgd)
2331 {
2332 register int i;
2333
2334 fprintf (fOut, "\nPerks:");
2335 for (i = pgd->startT; i <= pgd->endT; i++) {
2336 fprintf (fOut, "\t%g", pgd->rgdPerks[i]);
2337 }
2338 fprintf (fOut, "\nCounts:");
2339 for (i = pgd->startT; i <= pgd->endT; i++) {
2340 fprintf (fOut, "\t%ld", pgd->rglCount[i]);
2341 }
2342 fprintf (fOut, "\nLnPi(i):");
2343 for (i = pgd->startT; i <= pgd->endT; i++) {
2344 fprintf (fOut, "\t%g", pgd->rgdlnPi[i]);
2345 }
2346 fprintf (fOut, "\nTried Jumps:\t");
2347 for (i = pgd->startT; i <= pgd->endT - 1; i++) {
2348 fprintf (fOut, "\t%ld", pgd->rglTransAttempts[i]);
2349 }
2350 fprintf (fOut, "\nAccepted Jumps:\t");
2351 for (i = pgd->startT; i <= pgd->endT - 1; i++) {
2352 fprintf (fOut, "\t%ld", pgd->rglTransAccepts[i]);
2353 }
2354 fprintf(fOut, "\n\n");
2355 fflush(fOut);
2356
2357 #ifdef ndef
2358 /* This can be computed by the user if s/he wishes */
2359 if (eGeyer == Geyer) {
2360 /* Geyer's proposed adjustment of pseudo-priors */
2361 BOOL bZeroes;
2362 for (i = pgd->startT; i <= pgd->endT; i++) { /* check for zero counts */
2363 bZeroes = (pgd->rglCount[i] == 0);
2364 if (bZeroes) break; /* a zero count found, stop */
2365 }
2366 if (bZeroes) {
2367 fprintf (fOut, "Adjusted LnPi(i) not computable, zero-counts found.\n");
2368 }
2369 else {
2370 fprintf (fOut, "Adjusted LnPi(i): ");
2371 for (i = pgd->startT; i <= pgd->endT; i++) {
2372 pgd->rgdlnPi[i] = pgd->rgdlnPi[i] - log(pgd->rglCount[i]);
2373 fprintf (fOut, "% 10.6lE", pgd->rgdlnPi[i]);
2374 }
2375 fprintf (fOut, "\n");
2376 }
2377 }
2378 #endif
2379
2380 } /* PrintTemperatureDiagnostics */
2381
2382
2383 /* Function -------------------------------------------------------------------
2384 ReadData
2385
2386 Reads data for an experiment. There need to be one data item
2387 per output specified, and in the order given by the Print or PrintStep
2388 statements.
2389 Upgraded by FB 27/03/99 to deal with improved data handling.
2390 Called from TraverseLevels.
2391 */
2392 void ReadData (PLEVEL plevel, char **args)
2393 {
2394 FILE *pfileData = (FILE *) args[0];
2395 POUTSPEC pos;
2396 int cDat, i, j;
2397
2398 if (plevel->pexpt == NULL)
2399 return;
2400
2401 pos = &(plevel->pexpt->os);
2402
2403 /* here the number of data must equal the number of outputs */
2404 cDat = pos->nOutputs;
2405 pos->prgdDataVals = InitpdVector (cDat);
2406 pos->pcData = InitiVector (cDat); /* FB 5/11/99 */
2407 pos->pszDataNames = (PSTR *) malloc (cDat * sizeof(PSTR));
2408 pos->phvar_dat = (HVAR *) malloc (cDat * sizeof(HVAR));
2409
2410 if (pos->prgdDataVals == NULL || pos->phvar_dat == NULL ||
2411 pos->pszDataNames == NULL || pos->pcData == NULL)
2412 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "ReadData()", NULL);
2413 else {
2414 pos->nData = cDat; /* FB 27/03/99: Set count of data */
2415
2416 /* scan all data values for data file */
2417 for (i = 0; i < cDat; i++) {
2418 /* allocate space for pos->prgdDataVals[i] */
2419 if ( !(pos->prgdDataVals[i] = InitdVector (pos->pcOutputTimes[i])))
2420 ReportError (NULL, RE_OUTOFMEM | RE_FATAL, "ReadData()", NULL);;
2421
2422 /* for each requested output time */
2423 for (j = 0; j < pos->pcOutputTimes[i]; j++)
2424 if (fscanf(pfileData, "%lg", &(pos->prgdDataVals[i][j])) == EOF) {
2425 printf ("Error: incorrect length for data file - Exiting\n");
2426 exit(0);
2427 }
2428 pos->pcData[i] = j; /* FB 5/11/99 */
2429 pos->phvar_dat[i] = pos->phvar_out[i];
2430 pos->pszDataNames[i] = pos->pszOutputNames[i];
2431
2432 } /* for i */
2433
2434 } /* else */
2435
2436 } /* ReadData */
2437
2438
2439 /* Function -------------------------------------------------------------------
2440 ReadDataFile
2441
2442 Reads in data from a specified data file, if specified.
2443 */
2444 void ReadDataFile (PANALYSIS panal)
2445 {
2446 if (panal->gd.szGdata) {
2447
2448 char c;
2449
2450 FILE *pfile = fopen (panal->gd.szGdata, "r");
2451
2452 if (!pfile) {
2453 printf ("Cannot open data file '%s'\n", panal->gd.szGdata);
2454 exit (0);
2455 }
2456
2457 /* skip the first line */
2458 do { c = getc(pfile); } while (c != '\n');
2459
2460 TraverseLevels (panal->pLevels[0], ReadData, pfile, NULL);
2461
2462 fclose (pfile);
2463 }
2464
2465 } /* ReadDataFile */
2466
2467
2468 /* Function -------------------------------------------------------------------
2469 ReadRestart
2470
2471 initialize the parameters by reading them in the restart file.
2472 */
2473 void ReadRestart (FILE *pfileRestart, long nThetas,
2474 PDOUBLE *pdTheta, PDOUBLE *pdSum, PDOUBLE **prgdSumProd,
2475 long *pIter)
2476 {
2477 register char c;
2478 register long i, j;
2479
2480 if (*pdTheta == NULL)
2481 if ( !(*pdTheta = InitdVector(nThetas)) )
2482 ReportRunTimeError (NULL, RE_OUTOFMEM | RE_FATAL, "ReadRestart");
2483
2484 if (*pdSum == NULL)
2485 if ( !(*pdSum = InitdVector(nThetas)) )
2486 ReportRunTimeError (NULL, RE_OUTOFMEM | RE_FATAL, "ReadRestart");
2487
2488 if (*prgdSumProd == NULL)
2489 if ( !(*prgdSumProd = InitdMatrix (nThetas, nThetas)) )
2490 ReportRunTimeError (NULL, RE_OUTOFMEM | RE_FATAL, "ReadRestart");
2491
2492 *pIter = -1;
2493
2494 for (i = 0; i < nThetas; i++) {
2495 (*pdSum)[i] = 0.0;
2496 for (j = 0; j < nThetas; j++)
2497 (*prgdSumProd)[i][j] = 0.0;
2498 }
2499
2500 /* skip the first line. This allows a MC output file to be used
2501 directly as a restart file. */
2502 do { c = getc(pfileRestart); } while (c != '\n');
2503
2504 /* as long as we have not reached the end of the file we keep reading lines
2505 and overwriting the thetas, they keep only their last value.
2506 We throw away first field, and we keep incrementing the global iteration
2507 counter iter:
2508 */
2509 while (!(feof (pfileRestart) ||
2510 (fscanf (pfileRestart, "%*s") == EOF))) {
2511 for (i = 0; i < nThetas; i++) {
2512 if ( fscanf(pfileRestart, "%lg", &((*pdTheta)[i])) == EOF ) {
2513 printf ("Error: incorrect length for restart file - Exiting\n");
2514 exit(0);
2515 }
2516 else { /* reading ok */
2517 /* update pdSum, the column sum of pdTheta */
2518 (*pdSum)[i] += (*pdTheta)[i];
2519 }
2520 }
2521
2522 /* Throw away remainder of line. This allows a MC output file to be used
2523 directly as a restart file. */
2524 do { c = getc(pfileRestart); } while (c != '\n');
2525
2526 /* update prgdSumProd */
2527 for (i = 0; i < nThetas; i++)
2528 for (j = 0; j < nThetas; j++)
2529 (*prgdSumProd)[i][j] += (*pdTheta)[i] * (*pdTheta)[j];
2530
2531
2532 /* increment pIter */
2533 *pIter = *pIter + 1;
2534
2535 } /* end while */
2536
2537 /* note that the theta returned is the last parameter set read */
2538
2539 fclose (pfileRestart);
2540
2541 } /* ReadRestart */
2542
2543
2544 /* Function -------------------------------------------------------------------
2545 ReadRestartTemper
2546
2547 initialize parameters and temperature stuff (pseudoprior, indexT) by reading
2548 them in the restart file.
2549 */
2550 void ReadRestartTemper (FILE *pfileRestart, long nThetas, int nPerks,
2551 PDOUBLE *pdTheta, PDOUBLE *pdSum, PDOUBLE **prgdSumProd,
2552 long *pIter, int *pindexT, double *pdlnPi)
2553 {
2554 register char c;
2555 register long i, j;
2556
2557 if (*pdTheta == NULL)
2558 if ( !(*pdTheta = InitdVector(nThetas)) )
2559 ReportRunTimeError (NULL, RE_OUTOFMEM | RE_FATAL, "ReadRestart");
2560
2561 if (*pdSum == NULL)
2562 if ( !(*pdSum = InitdVector(nThetas)) )
2563 ReportRunTimeError (NULL, RE_OUTOFMEM | RE_FATAL, "ReadRestart");
2564
2565 if (*prgdSumProd == NULL)
2566 if ( !(*prgdSumProd = InitdMatrix (nThetas, nThetas)) )
2567 ReportRunTimeError (NULL, RE_OUTOFMEM | RE_FATAL, "ReadRestart");
2568
2569 *pIter = -1;
2570
2571 for (i = 0; i < nThetas; i++) {
2572 (*pdSum)[i] = 0.0;
2573 for (j = 0; j < nThetas; j++)
2574 (*prgdSumProd)[i][j] = 0.0;
2575 }
2576
2577 /* skip the first line. This allows a MC output file to be used
2578 directly as a restart file. */
2579 do { c = getc(pfileRestart); } while (c != '\n');
2580
2581 /* as long as we have not reached the end of the file we keep reading lines
2582 and overwriting the thetas, they keep only their last value.
2583 We throw away first field, and we keep incrementing the global iteration
2584 counter iter: */
2585 while (!(feof (pfileRestart) ||
2586 (fscanf (pfileRestart, "%*s") == EOF))) {
2587 for (i = 0; i < nThetas; i++) {
2588 if ( fscanf(pfileRestart, "%lg", &((*pdTheta)[i])) == EOF ) {
2589 printf ("Error: incorrect length for restart file - Exiting\n");
2590 exit(0);
2591 }
2592 else { /* reading ok */
2593 /* update pdSum, the column sum of pdTheta */
2594 (*pdSum)[i] += (*pdTheta)[i];
2595 }
2596 }
2597
2598 if (fscanf(pfileRestart,"%d", pindexT) == EOF) {
2599 printf ("Error: incorrect length for restart file - Exiting\n");
2600 exit(0);
2601 }
2602
2603 for (i = 0; i < nPerks; i++) {
2604 if (fscanf(pfileRestart,"%lg", &(pdlnPi[i])) == EOF) {
2605 printf ("Error: incorrect length for restart file - Exiting\n");
2606 exit(0);
2607 }
2608 }
2609
2610 /* Throw away remainder of line. This allows a MC output file to be used
2611 directly as a restart file. */
2612 do { c = getc(pfileRestart); } while (c != '\n');
2613
2614 /* update prgdSumProd */
2615 for (i = 0; i < nThetas; i++)
2616 for (j = 0; j < nThetas; j++)
2617 (*prgdSumProd)[i][j] += (*pdTheta)[i] * (*pdTheta)[j];
2618
2619 /* increment pIter */
2620 *pIter = *pIter + 1;
2621
2622 } /* end while */
2623
2624 /* note that the theta returned is the last parameter set read */
2625
2626 fclose (pfileRestart);
2627
2628 } /* ReadRestartTemper */
2629
2630
2631 /* Function -------------------------------------------------------------------
2632 RestoreLikelihoods
2633
2634 Called from TraverseLevels1
2635 */
2636 int RestoreLikelihoods (PLEVEL plevel, char **args)
2637 {
2638 PEXPERIMENT pExpt = plevel->pexpt;
2639
2640 if (pExpt != NULL) {
2641 pExpt->dLnLike = pExpt->dLnLikeSave;
2642 }
2643
2644 return (1);
2645
2646 } /* RestoreLikelihoods */
2647
2648
2649 /* Function -------------------------------------------------------------------
2650 RunAllExpts
2651
2652 Run all experiments (i.e. experiments at all levels below
2653 panal->pLevels[0]).
2654 */
2655 int RunAllExpts (PANALYSIS panal, PDOUBLE pdLnData)
2656 {
2657 PLEVEL plevel0 = panal->pLevels[0];
2658 long n;
2659
2660 for (n = 0; n < plevel0->iInstances; n++) {
2661 if (!TraverseLevels1 (plevel0->pLevels[n], RunExpt, panal,
2662 pdLnData, NULL)) {
2663 /* error */
2664 return(0);
2665 }
2666 }
2667
2668 return(1);
2669
2670 } /* RunAllExpts */
2671
2672
2673 /* Function -------------------------------------------------------------------
2674 RunExpt
2675
2676 If `plevel' has experiments, modify the variables that have been sampled
2677 at this level and above using its two lists, and run the experiments.
2678 Return 1 on success and 0 if failure.
2679 Computes the data likelihood in case of success.
2680
2681 Called from TraverseLevels1
2682 */
2683 int RunExpt (PLEVEL plevel, char **args)
2684 {
2685 PANALYSIS panal = (PANALYSIS) args[0];
2686 double *pdLnData = (double *) args[1];
2687 long i;
2688 PEXPERIMENT pExpt = plevel->pexpt;
2689
2690 /* Set level sequence */
2691 panal->pCurrentLevel[plevel->iDepth] = plevel;
2692
2693 if (pExpt != NULL) {
2694 InitModel ();
2695
2696 /* Set the model vars that have been sampled in this level and above */
2697 for (i = 0; i <= plevel->iDepth; i++) {
2698 SetModelVars (panal->pCurrentLevel[i]);
2699 SetFixedVars (panal->pCurrentLevel[i]);
2700 }
2701
2702 if (!DoOneExperiment (pExpt)) {
2703 /* Error */
2704 printf ("Warning: DoOneExperiment failed\n");
2705 return 0;
2706 }
2707 else {
2708 pExpt->dLnLike = LnLikeData (plevel, panal);
2709 *pdLnData = (*pdLnData) + pExpt->dLnLike;
2710 }
2711 } /* if */
2712
2713 return (1);
2714
2715 } /* RunExpt */
2716
2717
2718 /* Function -------------------------------------------------------------------
2719 SampleTemperature
2720 */
2721 long SampleTemperature (PGIBBSDATA pgd, double dLnPrior, double dLnData)
2722 {
2723 int indexT = pgd->indexT;
2724 int indexT_new;
2725
2726 /* Propose a new perk */
2727 if (indexT == 0) indexT_new = 1;
2728 else {
2729 if (indexT == pgd->nPerks - 1) indexT_new = indexT - 1;
2730 else {
2731 if (Randoms() > 0.5) indexT_new = indexT + 1;
2732 else indexT_new = indexT - 1;
2733 }
2734 }
2735
2736 /* Test the temperature */
2737 if (TestTemper (pgd, indexT, indexT_new, dLnPrior, dLnData,
2738 pgd->rgdlnPi[indexT], pgd->rgdlnPi[indexT_new])) {
2739 /* jump */
2740 return (indexT_new);
2741 }
2742 else
2743 return (indexT);
2744
2745 } /* SampleTemperature */
2746
2747
2748 /* Function -------------------------------------------------------------------
2749 SampleTemperature2
2750
2751 Records the count of attempted and accepted temperature jumps
2752 */
2753 long SampleTemperature2 (PGIBBSDATA pgd, double dLnPrior, double dLnData)
2754 {
2755 int indexT = pgd->indexT;
2756 int indexT_new;
2757
2758 /* Propose a new perk */
2759 if (indexT == pgd->startT) indexT_new = indexT + 1;
2760 else {
2761 if (indexT == pgd->endT) indexT_new = indexT - 1;
2762 else {
2763 if (Randoms() > 0.5) indexT_new = indexT + 1;
2764 else indexT_new = indexT - 1;
2765 }
2766 }
2767
2768 int minI = (indexT < indexT_new ? indexT : indexT_new);
2769 pgd->rglTransAttempts[minI]++;
2770
2771 /* Test the temperature */
2772 if (TestTemper (pgd, indexT, indexT_new, dLnPrior, dLnData,
2773 pgd->rgdlnPi[indexT], pgd->rgdlnPi[indexT_new])) {
2774 /* jump */
2775 pgd->rglTransAccepts[minI]++;
2776 return (indexT_new);
2777 }
2778 else
2779 return (indexT);
2780
2781 } /* SampleTemperature2 */
2782
2783
2784 /* Function -------------------------------------------------------------------
2785 SampleTheta, samples from a normal kernel
2786 */
2787 double SampleTheta (PMCVAR pMCVar) {
2788
2789 /* if the parameter has a discrete distribution, round it - FB 12/06/97 */
2790 if (pMCVar->iType == MCV_BINOMIAL || pMCVar->iType == MCV_POISSON) {
2791 return floor (0.5 + TruncNormalRandom (pMCVar->dVal, pMCVar->dKernelSD,
2792 MinMCVar(pMCVar),
2793 MaxMCVar(pMCVar)));
2794 }
2795 else { /* FB fixed the uniform case - FB 30/06/97 */
2796 return TruncNormalRandom (pMCVar->dVal, pMCVar->dKernelSD,
2797 MinMCVar(pMCVar), MaxMCVar(pMCVar));
2798 }
2799
2800 } /* SampleTheta */
2801
2802
2803 /* Function -------------------------------------------------------------------
2804 SampleThetaUnif, samples from a uniform kernel
2805 */
2806 double SampleThetaUnif (PMCVAR pMCVar) {
2807
2808 /* if the parameter has a discrete distribution, round it */
2809 if (pMCVar->iType == MCV_BINOMIAL || pMCVar->iType == MCV_POISSON)
2810 return floor (0.5 + UniformRandom (MinMCVar(pMCVar), MaxMCVar(pMCVar)));
2811 else
2812 return UniformRandom (MinMCVar(pMCVar), MaxMCVar(pMCVar));
2813
2814 } /* SampleThetaUnif */
2815
2816
2817 /* Function -------------------------------------------------------------------
2818 SampleThetas
2819
2820 Perform a Metropolis step on all the random (MC) variables at the
2821 level passed as argument 1.
2822 Sample thetas in sequence - test using prior and likelihood -
2823 restore old values if necessary - write new values to output file.
2824
2825 Called from TraverseLevels
2826 */
2827 void SampleThetas (PLEVEL plevel, char **args)
2828 {
2829 PANALYSIS panal = (PANALYSIS) args[0];
2830 PGIBBSDATA pgd = (PGIBBSDATA) args[1];
2831 long *pnIter = (long *) args[2];
2832 long *pnUpdateAt = (long *) args[3];
2833 long *pnTotal = (long *) args[4];
2834
2835 double dLnPrior, dLnLike, dLnData, dLnKern;
2836 double dLnPriorNew, dLnLikeNew, dLnDataNew, dLnKernNew;
2837 double dTheta, dJumps;
2838 PMCVAR pMCVar;
2839 long n;
2840
2841 /* Set level sequence. This is needed to call RunExpt from the middle
2842 of the tree with the right initialization of the nodes above */
2843 panal->pCurrentLevel[plevel->iDepth] = plevel;
2844
2845 /* For all MC vars at this level */
2846 for (n = 0; n < plevel->nMCVars; n++) {
2847
2848 pMCVar = plevel->rgpMCVars[n];
2849
2850 /* If the MC var is fixed, no sampling is made, just write it out */
2851 if (pMCVar->bIsFixed)
2852 goto WriteIt;
2853
2854 /* Compute prior and likelihood */
2855 dLnPrior = LnDensity (pMCVar, panal);
2856 dLnLike = LnLike (pMCVar, panal);
2857
2858 dLnData = 0.0;
2859
2860 /* If data are dependent compute the data likelihood */
2861 if (pMCVar->bExptIsDep) {
2862 /* Form the likelihood of all experiments at this level or beneath.*/
2863 TraverseLevels1 (plevel, SumAllExpts, &dLnData, NULL);
2864 }
2865
2866 /* Save current value */
2867 dTheta = pMCVar->dVal;
2868
2869 /* Adjust the jumping kernel SD, depending on acceptance rates,
2870 make sure it does not exceed DBL_MAX or a third of the range */
2871 if (*pnIter == *pnUpdateAt) {
2872
2873 dJumps = (double) pMCVar->lJumps / (double) (*pnTotal);
2874
2875 if (dJumps > 0.3) { /* we will increase the kernel spread */
2876 if (dJumps == 1) {
2877 /* pathological case: chances are the parameter has such a
2878 small kernel that it does not change values, this leads
2879 to a "jump" each time. Drastic increase of the kernel is
2880 attempted */
2881 if (pMCVar->dKernelSD < sqrt(DBL_MAX)) {
2882 if (pMCVar->dKernelSD > 2)
2883 pMCVar->dKernelSD = pMCVar->dKernelSD * pMCVar->dKernelSD;
2884 else /* FB 28/03/1999 */
2885 pMCVar->dKernelSD = pMCVar->dKernelSD * 20;
2886 }
2887 else
2888 pMCVar->dKernelSD = pMCVar->dMaxKernelSD;
2889 }
2890 else { /* more normal case */
2891 if (pMCVar->dKernelSD < DBL_MAX / 2)
2892 pMCVar->dKernelSD = pMCVar->dKernelSD * 2;
2893 else
2894 pMCVar->dKernelSD = pMCVar->dMaxKernelSD;
2895 }
2896
2897 /* check that kernel SD does not increase wildly */
2898 if (pMCVar->dKernelSD > pMCVar->dMaxKernelSD)
2899 pMCVar->dKernelSD = pMCVar->dMaxKernelSD;
2900 }
2901 else { /* we will decrease the kernel spread */
2902 if (dJumps == 0) {
2903 /* pathological case: chances are the parameter has such a
2904 big kernel that it will never jump. Drastic decrease of
2905 the kernel is attempted, dissymetric from the increase */
2906 if (pMCVar->dKernelSD > pow(DBL_MIN, 0.45)) {
2907 if (pMCVar->dKernelSD > 2) { /* FB 28/03/1999 */
2908 pMCVar->dKernelSD = pow(pMCVar->dKernelSD, 0.45);
2909 }
2910 else {
2911 pMCVar->dKernelSD = pMCVar->dKernelSD * 0.04;
2912 }
2913 }
2914 else
2915 pMCVar->dKernelSD = DBL_MIN;
2916 }
2917 else { /* more normal case */
2918 /* the kernel should not be infinitely decreased; decrease
2919 is otherwise slighly different from the increase to avoid
2920 oscillations */
2921 if (pMCVar->dKernelSD > DBL_MIN / 0.4)
2922 pMCVar->dKernelSD = pMCVar->dKernelSD * 0.4;
2923 else
2924 pMCVar->dKernelSD = DBL_MIN;
2925 }
2926 }
2927
2928 pMCVar->lJumps = 0; /* reset the jumps counter */
2929
2930 } /* end of kernel stuff */
2931
2932 /* compute the integral of the jump kernel */
2933 dLnKern = log(CDFNormal ((MaxMCVar(pMCVar) - pMCVar->dVal) /
2934 pMCVar->dKernelSD) -
2935 CDFNormal ((MinMCVar(pMCVar) - pMCVar->dVal) /
2936 pMCVar->dKernelSD));
2937
2938 /* sample a new value */
2939 pMCVar->dVal = SampleTheta (pMCVar);
2940
2941 /* update the integral of the jump kernel */
2942 dLnKernNew = log(CDFNormal ((MaxMCVar(pMCVar) - pMCVar->dVal) /
2943 pMCVar->dKernelSD) -
2944 CDFNormal ((MinMCVar(pMCVar) - pMCVar->dVal) /
2945 pMCVar->dKernelSD));
2946
2947 /* update prior and likelihood */
2948 dLnPriorNew = LnDensity (pMCVar, panal);
2949 dLnLikeNew = LnLike (pMCVar, panal);
2950
2951 dLnDataNew = 0.0;
2952
2953 /* If data are dependent compute the data likelihood */
2954 if (pMCVar->bExptIsDep) {
2955 /* Run all experiments at this level or beneath.
2956 We should in fact run only the dependent experiments ! */
2957
2958 if (!TraverseLevels1 (plevel, RunExpt, panal, &dLnDataNew, NULL)) {
2959 /* If running experiments fails, do not jump */
2960 pMCVar->dVal = dTheta;
2961 TraverseLevels1 (plevel, RestoreLikelihoods, NULL);
2962 goto WriteIt;
2963 }
2964 }
2965
2966 /* Test the results and act accordingly */
2967 if (!TestImpRatio (pgd, pMCVar->bExptIsDep, dLnKern, dLnKernNew,
2968 dLnPrior, dLnPriorNew,
2969 dLnLike, dLnLikeNew, dLnData, dLnDataNew)) {
2970 /* reject, restore */
2971 pMCVar->dVal = dTheta;
2972 if (pMCVar->bExptIsDep)
2973 TraverseLevels1 (plevel, RestoreLikelihoods, NULL);
2974 }
2975 else {
2976 /* accept, save likelihoods */
2977 pMCVar->lJumps = pMCVar->lJumps + 1;
2978
2979 if(pMCVar->bExptIsDep)
2980 TraverseLevels1 (plevel, SaveLikelihoods, NULL);
2981 }
2982
2983 /* calculate the running mean and variance for this value */
2984 CalculateMeanAndVariance((*pnIter+1),pMCVar->dVal,&pMCVar->dVal_mean,
2985 &pMCVar->dVal_var);
2986
2987 WriteIt: /* Write the MC var value to output file */
2988
2989 if (((*pnIter+1) % pgd->nPrintFreq == 0) &&
2990 (*pnIter >= pgd->nMaxIter - pgd->nPrintIter)) {
2991 fprintf(pgd->pfileOut, "%5g\t", pMCVar->dVal);
2992 }
2993 } /* end for all pMCVar */
2994
2995 } /* SampleThetas */
2996
2997
2998 /* Function -------------------------------------------------------------------
2999 SampleThetasTempered
3000
3001 Use annealing MCMC
3002 Sample thetas in sequence - test using prior and likelihood -
3003 restore old values if necessary - write new values to output file
3004 Called from TraverseLevels
3005
3006 */
3007 void SampleThetasTempered (PLEVEL plevel, char **args)
3008 {
3009 PANALYSIS panal = (PANALYSIS) args[0];
3010 PGIBBSDATA pgd = (PGIBBSDATA) args[1];
3011 long *pnIter = (long *) args[2];
3012 long *pnUpdateAt = (long *) args[3];
3013 long *pnTotal = (long *) args[4];
3014 long *pindexT = (long *) args[5];
3015
3016 double dLnPrior, dLnLike, dLnData, dLnKern;
3017 double dLnPriorNew, dLnLikeNew, dLnDataNew, dLnKernNew;
3018 double dTheta, dJumps, old_dKernelSD;
3019 PMCVAR pMCVar;
3020 long n;
3021
3022 /* Set level sequence. This is needed to call RunExpt from the middle
3023 of the tree with the right initialization of the nodes above */
3024 panal->pCurrentLevel[plevel->iDepth] = plevel;
3025
3026 /* For all MC vars at this level */
3027 for (n = 0; n < plevel->nMCVars; n++) {
3028
3029 pMCVar = plevel->rgpMCVars[n];
3030
3031 /* If the MC var is fixed, no sampling is made, just write it out */
3032 if (pMCVar->bIsFixed)
3033 goto WriteIt;
3034
3035 /* Compute prior and likelihood */
3036 dLnPrior = LnDensity (pMCVar, panal);
3037 dLnLike = LnLike (pMCVar, panal);
3038
3039 dLnData = 0.0;
3040
3041 /* If data are dependent compute the data likelihood */
3042 if (pMCVar->bExptIsDep) {
3043 /* Form the likelihood of all experiments at this level or beneath.*/
3044 TraverseLevels1 (plevel, SumAllExpts, &dLnData, NULL);
3045 }
3046
3047 /* Save current value */
3048 dTheta = pMCVar->dVal;
3049
3050 /* Adjust the jumping kernel SD, depending on acceptance rates,
3051 make sure it does not exceed DBL_MAX or a third of the range */
3052 if (*pnIter == *pnUpdateAt) {
3053
3054 dJumps = (double) pMCVar->lJumps / (double) (*pnTotal);
3055
3056 if (dJumps > 0.3) { /* we will increase the kernel spread */
3057 if (dJumps == 1) {
3058 /* pathological case: chances are the parameter has such a
3059 small kernel that it does not change values, this leads
3060 to a "jump" each time. Drastic increase of the kernel is
3061 attempted */
3062 if (pMCVar->dKernelSD < sqrt(DBL_MAX)) {
3063 if (pMCVar->dKernelSD > 2)
3064 pMCVar->dKernelSD = pMCVar->dKernelSD * pMCVar->dKernelSD;
3065 else /* FB 28/03/1999 */
3066 pMCVar->dKernelSD = pMCVar->dKernelSD * 20;
3067 }
3068 else
3069 pMCVar->dKernelSD = DBL_MAX;
3070 }
3071 else { /* more normal case */
3072 if (pMCVar->dKernelSD < DBL_MAX / 2)
3073 pMCVar->dKernelSD = pMCVar->dKernelSD * 2;
3074 else
3075 pMCVar->dKernelSD = DBL_MAX;
3076 }
3077
3078 /* check that kernel SD does not increase wildly */
3079 if (pMCVar->dKernelSD > pMCVar->dMaxKernelSD)
3080 pMCVar->dKernelSD = pMCVar->dMaxKernelSD;
3081 }
3082 else { /* we will decrease the kernel spread */
3083 if (dJumps == 0) {
3084 /* pathological case: chances are the parameter has such a
3085 big kernel that it will never jump. Drastic decrease of
3086 the kernel is attempted, dissymetric from the increase */
3087 if (pMCVar->dKernelSD > pow(DBL_MIN, 0.45)) {
3088 if (pMCVar->dKernelSD > 2) { /* FB 28/03/1999 */
3089 pMCVar->dKernelSD = pow(pMCVar->dKernelSD, 0.45);
3090 }
3091 else {
3092 pMCVar->dKernelSD = pMCVar->dKernelSD * 0.04;
3093 }
3094 }
3095 else
3096 pMCVar->dKernelSD = DBL_MIN;
3097 }
3098 else { /* more normal case */
3099 /* the kernel should not be infinitely decreased; decrease
3100 is otherwise slighly different from the increase to avoid
3101 oscillations */
3102 if (pMCVar->dKernelSD > DBL_MIN / 0.4)
3103 pMCVar->dKernelSD = pMCVar->dKernelSD * 0.4;
3104 else
3105 pMCVar->dKernelSD = DBL_MIN;
3106 }
3107 }
3108
3109 pMCVar->lJumps = 0; /* reset the jumps counter */
3110
3111 } /* end kernel updating */
3112
3113 /* sample a new value */
3114
3115 if (pgd->rgdPerks[*pindexT] > 0) { /* usual case */
3116
3117 /* first scale temporarily the kernelSD according to the temperature */
3118
3119 /* save the current kernelSD */
3120 old_dKernelSD = pMCVar->dKernelSD;
3121
3122 /* update the kernelSD according to the temperature, the formula
3123 comes from that of the variance of the powered standard normal */
3124 pMCVar->dKernelSD = pMCVar->dKernelSD *
3125 exp((1 - pgd->rgdPerks[*pindexT]) * LN2PI * 0.25 -
3126 0.75 * log(pgd->rgdPerks[*pindexT]));
3127
3128 /* check that kernel SD does not increase wildly */
3129 if (pMCVar->dKernelSD > pMCVar->dMaxKernelSD)
3130 pMCVar->dKernelSD = pMCVar->dMaxKernelSD;
3131
3132 /* compute the integral of the jump kernel (for truncation) */
3133 dLnKern = log(CDFNormal((MaxMCVar(pMCVar) - pMCVar->dVal) /
3134 pMCVar->dKernelSD) -
3135 CDFNormal((MinMCVar(pMCVar) - pMCVar->dVal) /
3136 pMCVar->dKernelSD));
3137
3138 /* sample a new value */
3139 pMCVar->dVal = SampleTheta (pMCVar);
3140
3141 /* update the integral of the jump kernel (for truncation) */
3142 dLnKernNew = log(CDFNormal((MaxMCVar(pMCVar) - pMCVar->dVal) /
3143 pMCVar->dKernelSD) -
3144 CDFNormal((MinMCVar(pMCVar) - pMCVar->dVal) /
3145 pMCVar->dKernelSD));
3146 }
3147 else {
3148 /* inverse temperature is zero, sample uniformly or from the prior,
3149 the jump will always be accepted */
3150 if (pgd->nSimTypeFlag == 3) { /* the posterior is tempered */
3151 /* sample a new value uniformly in its range */
3152 pMCVar->dVal = SampleThetaUnif(pMCVar);
3153 }
3154 else { /* the likelihood is tempered, sample from the prior */
3155 CalculateOneMCParm(pMCVar);
3156 }
3157 }
3158
3159 /* recompute prior and likelihood */
3160 dLnPriorNew = LnDensity (pMCVar, panal);
3161 dLnLikeNew = LnLike (pMCVar, panal);
3162
3163 dLnDataNew = 0.0;
3164
3165 /* If data are dependent compute the data likelihood */
3166 if (pMCVar->bExptIsDep) {
3167 /* Run all experiments beneath this level.
3168 We should in fact run only the dependent experiments ! */
3169 if (!TraverseLevels1 (plevel, RunExpt, panal, &dLnDataNew, NULL)) {
3170 /* If running experiments fails, do not jump */
3171 pMCVar->dVal = dTheta; /* restore */
3172 TraverseLevels1 (plevel, RestoreLikelihoods, NULL);
3173 goto WriteIt;
3174 }
3175 }
3176
3177 /* Test the results and act accordingly */
3178 if (!TestImpRatioTemper (pgd, pMCVar->bExptIsDep, dLnKern, dLnKernNew,
3179 dLnPrior, dLnPriorNew,
3180 dLnLike, dLnLikeNew, dLnData, dLnDataNew,
3181 *pindexT)) {
3182 /* reject, restore */
3183 pMCVar->dVal = dTheta;
3184 if (pMCVar->bExptIsDep)
3185 TraverseLevels1 (plevel, RestoreLikelihoods, NULL);
3186 }
3187 else {
3188 /* accept, save likelihoods */
3189 pMCVar->lJumps = pMCVar->lJumps + 1;
3190 if (pMCVar->bExptIsDep)
3191 TraverseLevels1 (plevel, SaveLikelihoods, NULL);
3192 }
3193
3194 if (pgd->rgdPerks[*pindexT] > 0) /* restore kernelSD */
3195 pMCVar->dKernelSD = old_dKernelSD;
3196
3197 WriteIt: /* Write the current MC var value to output file */
3198
3199 if (((*pnIter+1) % pgd->nPrintFreq == 0) &&
3200 (*pnIter >= pgd->nMaxIter - pgd->nPrintIter)) {
3201 fprintf(pgd->pfileOut, "%5g\t", pMCVar->dVal);
3202 }
3203 } /* Metrolopolis-Hastings */
3204
3205 } /* SampleThetasTempered */
3206
3207
3208 /* Function -------------------------------------------------------------------
3209 SampleThetaVector
3210
3211 Sample thetas in block, do Metropolis test.
3212 */
3213 void SampleThetaVector (PLEVEL pLevel, PANALYSIS panal, long nThetas,
3214 double *pdTheta, double *pdSum, double **prgdSumProd,
3215 long iter, long nUpdateAt, long nTotal,
3216 PDOUBLE pdLnPrior, PDOUBLE pdLnData)
3217 {
3218 register long i, j;
3219 double dTmp, dAccept, dLnPrior_old, dLnData_old;
3220 BOOL bInBounds;
3221
3222 static long lAccepted = 0;
3223 static double dJumpSpread;
3224 static PDOUBLE pdTheta_old = NULL; /* previous model parameters values */
3225 static PDOUBLE *prgdComponent;
3226 static PDOUBLE *prgdVariance;
3227 static PDOUBLE dNormVar; /* storage for nParms normal deviates */
3228
3229 if ((pdTheta_old == NULL) || (iter == nUpdateAt)) {
3230
3231 if (pdTheta_old == NULL) { /* initialize */
3232 if ( !(pdTheta_old = InitdVector (nThetas)) ||
3233 !(dNormVar = InitdVector (nThetas)) ||
3234 !(prgdVariance = InitdMatrix (nThetas, nThetas)) ||
3235 !(prgdComponent = InitdMatrix (nThetas, nThetas)))
3236 ReportRunTimeError (panal, RE_OUTOFMEM | RE_FATAL,
3237 "SampleThetaVector");
3238
3239 /* initialize dJumpSpread */
3240 dJumpSpread = 2.4 / sqrt(nThetas); /* Gelman's Normal theory result */
3241 }
3242 else {
3243 /* done if iter == nUpdateAt, but not at start:
3244 if some vector samplings have been made, check that the current
3245 dJumpSpread leads to an acceptation rate of 15 to 30% over the
3246 last batch of simulations. Adjust eventually. */
3247 dAccept = ((double) lAccepted) / (double) (nTotal);
3248
3249 if ( dAccept > 0.3) dJumpSpread = dJumpSpread * 1.5;
3250 else if (dAccept < 0.15) dJumpSpread = dJumpSpread * 0.7;
3251
3252 printf ("Monitoring: iter\t%ld\t", iter);
3253 printf ("success rate\t%g\tspread\t%g\n", dAccept, dJumpSpread);
3254 lAccepted = 0; /* reset the counter */
3255 }
3256
3257 /* other updates: */
3258
3259 /* generate the covariance matrix */
3260 for (i = 0; i < nThetas; i++)
3261 for (j = 0; j < nThetas; j++)
3262 prgdVariance[i][j] = (prgdSumProd[i][j] -
3263 pdSum[i] * pdSum[j] / (double) (iter+1)) /
3264 (double) iter;
3265
3266 /* do the Cholesky decomposition of the covariance matrix */
3267 if (!Cholesky (prgdVariance, prgdComponent, nThetas)) {
3268 /* try to save the computation by zeroing all non-diagonal elements */
3269 for (i = 0; i < nThetas; i++)
3270 for (j = 0; j < nThetas; j++) {
3271 if (i == j)
3272 prgdVariance[i][j] = prgdSumProd[i][j] / (double) (iter);
3273 else
3274 prgdVariance[i][j] = 0.0;
3275 }
3276
3277 /* if it still does not work, exit */
3278 if (!Cholesky (prgdVariance, prgdComponent, nThetas)) {
3279 printf ("Error: impossible to compute a jumping kernel - Exiting."
3280 "(You should check or change the restart file).\n\n");
3281 exit (0);
3282 }
3283 }
3284 }
3285
3286 /* keep the value of all thetas */
3287 for (i = 0; i < nThetas; i++)
3288 pdTheta_old[i] = pdTheta[i];
3289
3290 /* keep old prior and old likelihood */
3291 dLnPrior_old = *pdLnPrior;
3292 dLnData_old = *pdLnData;
3293
3294 /* generate new pdTheta vector */
3295 for (i = 0; i < nThetas; i++)
3296 dNormVar[i] = NormalRandom(0.0, 1.0);
3297
3298 for (i = 0; i < nThetas; i++) {
3299 dTmp = 0;
3300 for (j = 0; j <= i; j++) /* only the non-zero part of prgdComponent */
3301 dTmp = dTmp + dNormVar[j] * prgdComponent[i][j];
3302
3303 pdTheta[i] = pdTheta_old[i] + dJumpSpread * dTmp;
3304 }
3305
3306 /* Set the dVals of the variables to the values sampled and check that
3307 we are within bounds */
3308 long iTmp = 0; /* use a dummy variable to avoid resetting nThetas */
3309 bInBounds = TraverseLevels1 (pLevel, SetMCVars, pdTheta, &iTmp, NULL);
3310
3311 if (!bInBounds) { /* reject */
3312 for (i = 0; i < nThetas; i++)
3313 pdTheta[i] = pdTheta_old[i]; /* restore */
3314 goto DontJump;
3315 }
3316
3317 /* Calculate the new prior */
3318 *pdLnPrior = 0.0;
3319 TraverseLevels (pLevel, CalculateTotals, panal, pdLnPrior, NULL);
3320
3321 /* compute the model at the newly drawn point and the likelihood */
3322 *pdLnData = 0.0;
3323 if (!RunAllExpts (panal, pdLnData)) {
3324 /* computation failed, don't jump, restore */
3325 for (i = 0; i < nThetas; i++)
3326 pdTheta[i] = pdTheta_old[i];
3327 *pdLnPrior = dLnPrior_old;
3328 *pdLnData = dLnData_old;
3329 }
3330 else {
3331 /* Test */
3332 if (log(Randoms()) >
3333 ((*pdLnPrior) + (*pdLnData) - dLnPrior_old - dLnData_old)) {
3334 /* don't jump, restore */
3335 for (i = 0; i < nThetas; i++)
3336 pdTheta[i] = pdTheta_old[i];
3337 *pdLnPrior = dLnPrior_old;
3338 *pdLnData = dLnData_old;
3339 }
3340 else { /* jump */
3341 lAccepted++; /* this is used above to adjust the acceptation rate */
3342 }
3343 }
3344
3345 DontJump:
3346
3347 /* update arrays */
3348 for (i = 0; i < nThetas; i++) {
3349 pdSum[i] += pdTheta[i];
3350 for (j = 0; j < nThetas; j++)
3351 prgdSumProd[i][j] += pdTheta[i] * pdTheta[j];
3352 }
3353
3354 } /* SampleThetaVector */
3355
3356
3357 /* Function -------------------------------------------------------------------
3358 SaveLikelihoods
3359
3360 Called from TraverseLevels1
3361 */
3362 int SaveLikelihoods (PLEVEL plevel, char **args)
3363 {
3364 PEXPERIMENT pExpt = plevel->pexpt;
3365
3366 if (pExpt != NULL) {
3367 pExpt->dLnLikeSave = pExpt->dLnLike;
3368 }
3369
3370 return (1);
3371
3372 } /* SaveLikelihoods */
3373
3374
3375 /* Function -------------------------------------------------------------------
3376 SetFixedVars
3377
3378 Set the array of fixed variables
3379 */
3380 void SetFixedVars (PLEVEL plevel)
3381 {
3382 long n;
3383 PVARMOD pFVar;
3384
3385 for (n = 0; n < plevel->nFixedVars; n++) {
3386 pFVar = plevel->rgpFixedVars[n];
3387 if (IsInput (pFVar->hvar))
3388 SetInput (pFVar->hvar, pFVar->uvar.pifn);
3389 else
3390 SetVar (pFVar->hvar, pFVar->uvar.dVal);
3391 }
3392
3393 } /* SetFixedVars */
3394
3395
3396 /* Function -------------------------------------------------------------------
3397 SetKernel
3398
3399 Set initial values of the MCMC jumping kernel and eventually
3400 initializes the sampled parameters (contained in plevel->rgpMCVars[]).
3401
3402 The first argument of the list is a flag indicating whether the
3403 values of the sampled parameters should be restored (to the values
3404 passed as a second arument in an array) (case 1), or left at their
3405 last sampled value (case 2).
3406 */
3407 void SetKernel (PLEVEL plevel, char **args)
3408 {
3409 intptr_t useMCVarVals = (intptr_t) args[0]; /* 1 to restore dVals, else 2 */
3410 double *pdMCVarVals = (double *) args[1];
3411 double dMin, dMax, dTmp;
3412 long n, m;
3413 static long nThetas;
3414 PMCVAR pMCVar;
3415
3416 /* set the jumping kernel's SD: sample 4 variates and take the range */
3417 for (n = 0; n < plevel->nMCVars; n++) {
3418
3419 if (!(plevel->rgpMCVars[n]->bIsFixed)) {
3420
3421 pMCVar = plevel->rgpMCVars[n];
3422 CalculateOneMCParm (pMCVar);
3423
3424 /* set the maximum kernel SD value to about half of the SD of a
3425 uniform distribution having the same range */
3426 if (pMCVar->iType == MCV_UNIFORM || pMCVar->iType == MCV_LOGUNIFORM )
3427 pMCVar->dMaxKernelSD = (*(pMCVar->pdParm[1]) / 6.0) -
3428 (*(pMCVar->pdParm[0]) / 6.0);
3429 else
3430 pMCVar->dMaxKernelSD = (*(pMCVar->pdParm[3]) / 6.0) -
3431 (*(pMCVar->pdParm[2]) / 6.0);
3432
3433 /* we want a kernel SD compatible with the SD of the prior; this could
3434 be done analytically, depending on the prior. We approximate it by
3435 sampling 4 variates from the prior */
3436 dMin = dMax = pMCVar->dVal;
3437 for (m = 0; m < 3; m++) {
3438 CalculateOneMCParm(pMCVar);
3439 dTmp = pMCVar->dVal;
3440 if (dMin >= dTmp) dMin = dTmp;
3441 else if (dMax < dTmp) dMax = dTmp;
3442 }
3443
3444 /* set the range safely because max - min could be too large */
3445 if ((*(pMCVar->pdParm[2]) == -DBL_MAX) ||
3446 (*(pMCVar->pdParm[3]) == DBL_MAX))
3447 pMCVar->dKernelSD = (0.5 * dMax) - (0.5 * dMin);
3448 else
3449 pMCVar->dKernelSD = dMax - dMin;
3450
3451 /* take care of the case in which the range is zero
3452 (can happens for discrete variables) */
3453 if (pMCVar->dKernelSD == 0) {
3454 pMCVar->dKernelSD = pMCVar->dMaxKernelSD;
3455 }
3456
3457 /* check that we do not exceed the maximum SD */
3458 if (pMCVar->dKernelSD > pMCVar->dMaxKernelSD) {
3459 pMCVar->dKernelSD = pMCVar->dMaxKernelSD;
3460 }
3461 }
3462
3463 /* restore the value of the variable - FB 02/07/97 */
3464 if (useMCVarVals == 1)
3465 plevel->rgpMCVars[n]->dVal = pdMCVarVals[nThetas++];
3466
3467 }
3468
3469 } /* SetKernel */
3470
3471
3472 /* Function -------------------------------------------------------------------
3473 WriteKernel
3474
3475 Write values of the MCMC jumping kernel
3476 */
3477 void WriteKernel (PLEVEL plevel, char **args)
3478 {
3479 FILE *pfile = (FILE *) args[0];
3480 long n;
3481 PMCVAR pMCVar;
3482
3483 for (n = 0; n < plevel->nMCVars; n++) {
3484 if ( !(plevel->rgpMCVars[n]->bIsFixed)) {
3485 pMCVar = plevel->rgpMCVars[n];
3486 fprintf(pfile, "%lg\t", pMCVar->dKernelSD);
3487 }
3488 }
3489
3490 } /* WriteKernel */
3491
3492
3493 /* Function -------------------------------------------------------------------
3494 ReadKernel
3495
3496 Read values of the MCMC jumping kernel
3497 */
3498 void ReadKernel (PLEVEL plevel, char **args)
3499 {
3500 FILE *pfile = (FILE *) args[0];
3501 long n;
3502 PMCVAR pMCVar;
3503
3504 for (n = 0; n < plevel->nMCVars; n++) {
3505 if ( !(plevel->rgpMCVars[n]->bIsFixed)) {
3506 pMCVar = plevel->rgpMCVars[n];
3507
3508 /* set the maximum kernel SD value */
3509 pMCVar->dMaxKernelSD = (MaxMCVar(pMCVar) - MinMCVar(pMCVar)) / 6.0;
3510
3511 if (!fscanf(pfile, "%lg", &(pMCVar->dKernelSD))) {
3512 ReportError (NULL, RE_READERROR | RE_FATAL, "kernel file", NULL);
3513 }
3514 }
3515 }
3516
3517 } /* ReadKernel */
3518
3519
3520 /* Function -------------------------------------------------------------------
3521 SetModelVars
3522
3523 Sets the array of model variables to the sampled values. Does not set fixed
3524 variables. That has to be done by SetFixedVars.
3525 */
3526 void SetModelVars(PLEVEL plevel)
3527 {
3528 long n;
3529 PMCVAR pMCVar;
3530
3531 for (n = 0; n < plevel->nMCVars; n++) {
3532 pMCVar = plevel->rgpMCVars[n];
3533 if (!(pMCVar->bIsFixed) && (IsParm (pMCVar->hvar)))
3534 SetVar (pMCVar->hvar, pMCVar->dVal);
3535 }
3536
3537 } /* SetModelVars */
3538
3539
3540 /* Function -------------------------------------------------------------------
3541 SetMCVars
3542
3543 Set initial values of thetas after reading input file or sampling them.
3544 Values are assumed to be in proper order.
3545 It also checks that the ranges are respected.
3546 Called from TraverseLevels.
3547 */
3548 int SetMCVars (PLEVEL plevel, char **args)
3549 {
3550 double *pdMCVarVals = (double *) args[0];
3551 long *nThetas = (long *) args[1];
3552 PMCVAR pMCVar;
3553 double dVar;
3554 long n;
3555
3556 for (n = 0; n < plevel->nMCVars; n++) {
3557 dVar = pdMCVarVals[(*nThetas)++];
3558 pMCVar = plevel->rgpMCVars[n];
3559 if ((pMCVar->iType == MCV_UNIFORM) || (pMCVar->iType == MCV_LOGUNIFORM)) {
3560 if ((dVar < *(pMCVar->pdParm[0])) || (dVar > *(pMCVar->pdParm[1]))) {
3561 /* error */
3562 return (0);
3563 }
3564 }
3565 else {
3566 if ((dVar < *(pMCVar->pdParm[2])) || (dVar > *(pMCVar->pdParm[3]))) {
3567 /* error */
3568 return (0);
3569 }
3570 }
3571
3572 pMCVar->dVal = dVar;
3573 }
3574
3575 /* success */
3576 return (1);
3577
3578 } /* SetMCVars */
3579
3580
3581 /* Function -------------------------------------------------------------------
3582 SetPointers
3583
3584 Called from TraverseLevels
3585
3586 FB 26 nov 96: For each Monte Carlo variable, pdParms are set to point to the
3587 parent's dVals rather than to model parameters. If there is no parent,
3588 pdParms point to their own dParms.
3589 The pdVals and pdParms of the data and output arrays are also set.
3590 */
3591 void SetPointers (PLEVEL plevel, char **args)
3592 {
3593 long i, j, k;
3594 PMCVAR pMCVar;
3595 POUTSPEC pos;
3596 BOOL bFound;
3597
3598 for (i = 0; i < plevel->nMCVars; i++) { /* for parameters */
3599 pMCVar = plevel->rgpMCVars[i];
3600
3601 /* For each distribution parameter */
3602 for (j = 0; j < 4; j++) {
3603 if (pMCVar->pMCVParent[j] == NULL) /* Point to its own values */
3604 pMCVar->pdParm[j] = &(pMCVar->dParm[j]);
3605 else /* Point to the parent dVal */
3606 pMCVar->pdParm[j] = &(pMCVar->pMCVParent[j]->dVal);
3607 }
3608 }
3609
3610 if (plevel->pexpt != NULL) /* if the level has experiments */
3611 for (i = 0; i < plevel->nLikes; i++) { /* for likelihoods */
3612 pMCVar = plevel->rgpLikes[i];
3613 pos = &(plevel->pexpt->os);
3614
3615 /* Set pdVal of pMCVar to a data array, first find it */
3616 bFound = FALSE;
3617 j = 0;
3618 while ((j < pos->nData) && (!bFound)) {
3619 bFound = (pMCVar->hvar == pos->phvar_dat[j]);
3620 if (!bFound) j++;
3621 }
3622 if (bFound) {
3623 pMCVar->pdVal = pos->prgdDataVals[j];
3624 pMCVar->lCount = pos->pcData[j]; /* number of data points */
3625 }
3626 else { /* no corresponding Data found: error */
3627 printf ("Error: no Data for %s in Simulation %d - Exiting.\n\n",
3628 pMCVar->pszName, plevel->pexpt->iExp);
3629 exit (0);
3630 }
3631
3632 /* we also have to set the pdParms of pMCVar to either data or output
3633 arrays */
3634 for (j = 0; j < 4; j++) {
3635
3636 if (pMCVar->iParmType[j] == MCVP_PRED) { /* it's an output */
3637 /* set pdParm[j] to an output, first find it */
3638 bFound = FALSE;
3639 k = 0;
3640 while ((k < pos->nOutputs) && (!bFound)) {
3641 bFound = (pMCVar->hParm[j] == pos->phvar_out[k]);
3642 if (!bFound) k++;
3643 }
3644 if (bFound) {
3645 pMCVar->pdParm[j] = &(pos->prgdOutputVals[k][0]);
3646 }
3647 else { /* no corresponding Print statement found: error */
3648 printf ("Error: missing Print statement for parameter number %ld\n"
3649 "of %s distribution - Exiting.\n\n", j, pMCVar->pszName);
3650 exit (0);
3651 }
3652 } /* end if MCVP_PRED */
3653
3654 else if (pMCVar->iParmType[j] == MCVP_DATA) { /* it's data */
3655 /* set pdParm[j] to a data array, first find it */
3656 bFound = FALSE;
3657 k = 0;
3658 while ((k < pos->nData) && (!bFound)) {
3659 bFound = (pMCVar->hParm[j] == pos->phvar_dat[k]);
3660 if (!bFound) k++;
3661 }
3662 if (bFound) {
3663 pMCVar->pdParm[j] = &(pos->prgdDataVals[k][0]);
3664 }
3665 else { /* no Data found: error */
3666 printf ("Error: no Data for %s in Simulation %d - Exiting.\n\n",
3667 pMCVar->pszName, plevel->pexpt->iExp);
3668 exit (0);
3669 }
3670 } /* end if MCVP_DATA */
3671
3672 else { /* it's either a parameter or a numeric */
3673 if (pMCVar->pMCVParent[j] == NULL) /* Point to its own values */
3674 pMCVar->pdParm[j] = &(pMCVar->dParm[j]);
3675 else /* Point to the parent dVal */
3676 pMCVar->pdParm[j] = &(pMCVar->pMCVParent[j]->dVal);
3677 }
3678 } /* end for j */
3679 }
3680
3681 } /* SetPointers */
3682
3683
3684 /* Function -------------------------------------------------------------------
3685 SumAllExpts
3686
3687 If `plevel' has experiments, add the current Likelihood to the total
3688 passed as argument 2.
3689
3690 Called from TraverseLevels1
3691 */
3692 int SumAllExpts (PLEVEL plevel, char **args)
3693 {
3694 double *pdLnData = (double*)args[0];
3695 PEXPERIMENT pExpt = plevel->pexpt;
3696
3697 if (pExpt != NULL) {
3698 *pdLnData += pExpt->dLnLike;
3699 }
3700 return (1);
3701
3702 } /* SumAllExpts */
3703
3704
3705 /* Function -------------------------------------------------------------------
3706 TestImpRatio
3707
3708 Test prior, likelihood against random number between 0 and 1
3709 */
3710 BOOL TestImpRatio (PGIBBSDATA pgd, BOOL bExptIsDep,
3711 double dLnKern, double dLnKernNew,
3712 double dLnPrior, double dLnPriorNew,
3713 double dLnLike, double dLnLikeNew,
3714 double dLnData, double dLnDataNew)
3715 {
3716 double dPjump;
3717
3718 if (dLnKernNew == NULL_SUPPORT ||
3719 dLnPriorNew == NULL_SUPPORT ||
3720 dLnLikeNew == NULL_SUPPORT ||
3721 dLnDataNew == NULL_SUPPORT)
3722 return FALSE;
3723
3724 dPjump = dLnPriorNew - dLnPrior + dLnLikeNew - dLnLike +
3725 dLnKern - dLnKernNew;
3726
3727 if (bExptIsDep)
3728 dPjump += dLnDataNew - dLnData;
3729
3730 if (pgd->nSimTypeFlag == 0)
3731 return ((BOOL) (log(Randoms()) <= dPjump));
3732 else {
3733 if (pgd->nSimTypeFlag == 5)
3734 return ((BOOL) (0 <= dPjump));
3735 else {
3736 printf ("Error: simTypeFlag should be 0 or 5 in TestImpRatio "
3737 "- Exiting.\n\n");
3738 exit (0);
3739 }
3740 }
3741
3742 } /* TestImpRatio */
3743
3744
3745 /* Function -------------------------------------------------------------------
3746 TestImpRatioTemper
3747
3748 Test prior, likelihood against a random number between 0 and 1 according to
3749 the temperature
3750 */
3751 BOOL TestImpRatioTemper (PGIBBSDATA pgd, BOOL bExptIsDep,
3752 double dLnKern, double dLnKernNew,
3753 double dLnPrior, double dLnPriorNew,
3754 double dLnLike, double dLnLikeNew,
3755 double dLnData, double dLnDataNew, long indexT)
3756 {
3757 double dPjump;
3758
3759 if (dLnPriorNew == NULL_SUPPORT ||
3760 dLnLikeNew == NULL_SUPPORT ||
3761 dLnDataNew == NULL_SUPPORT)
3762 return FALSE;
3763
3764 if (pgd->rgdPerks[indexT] == 0) /* always accept jumps at perk zero */
3765 return TRUE;
3766
3767 if (pgd->nSimTypeFlag == 3) { /* posterior is tempered */
3768 dPjump = pgd->rgdPerks[indexT] *
3769 (dLnPriorNew - dLnPrior + dLnLikeNew - dLnLike) +
3770 dLnKern - dLnKernNew;
3771 }
3772 else { /* only the likelihood is tempered */
3773 dPjump = dLnPriorNew - dLnPrior +
3774 pgd->rgdPerks[indexT] * (dLnLikeNew - dLnLike) +
3775 dLnKern - dLnKernNew;
3776 }
3777
3778 if (bExptIsDep)
3779 dPjump += pgd->rgdPerks[indexT] * (dLnDataNew - dLnData);
3780
3781 return ((BOOL) (log(Randoms()) <= dPjump));
3782
3783 } /* TestImpRatioTemper */
3784
3785
3786 /* Function -------------------------------------------------------------------
3787 TestTemper
3788
3789 Test temperature against a random number between 0 and 1
3790 */
3791 BOOL TestTemper (PGIBBSDATA pgd, long indexT, long indexT_new, double dLnPrior,
3792 double dLnData, double pseudo, double pseudonew)
3793 {
3794 double dPjump = 0;
3795 #define MINUSLN2 -0.6931471805599452862268
3796
3797 if (dLnPrior + dLnData == NULL_SUPPORT)
3798 return FALSE;
3799
3800 if (pgd->nSimTypeFlag == 3) { /* the posterior is tempered */
3801 dPjump = (pgd->rgdPerks[indexT_new] -
3802 pgd->rgdPerks[indexT]) * (dLnPrior + dLnData) +
3803 pseudonew - pseudo +
3804 ((indexT_new == 0) || (indexT_new == pgd->nPerks - 1) ?
3805 0 : MINUSLN2) -
3806 ((indexT == 0) || (indexT == pgd->nPerks - 1) ?
3807 0 : MINUSLN2);
3808 }
3809 else { /* only the likelihood is tempered, the prior cancels out */
3810 dPjump = (pgd->rgdPerks[indexT_new] -
3811 pgd->rgdPerks[indexT]) * dLnData +
3812 pseudonew - pseudo +
3813 ((indexT_new == 0) || (indexT_new == pgd->nPerks - 1) ?
3814 0 : MINUSLN2) -
3815 ((indexT == 0) || (indexT == pgd->nPerks - 1) ?
3816 0 : MINUSLN2);
3817 }
3818
3819 return ((BOOL) (log(Randoms()) <= dPjump));
3820
3821 } /* TestTemper */
3822
3823
3824 /* Function -------------------------------------------------------------------
3825 TraverseLevels (recursive)
3826
3827 Called with variable argument list ending in NULL;
3828 arguments should be pointers only; if you call this with a value
3829 that can be zero, you will be very sorry
3830
3831 Find all allocated levels, execute `routinePtr' for each, starting at the
3832 top, passing the argument list as char**
3833
3834 The argument list is saved from the initial call; on recursive calls the
3835 list is NULL
3836 */
3837 void TraverseLevels (PLEVEL plevel,
3838 void (*routinePtr)(PLEVEL plevel, char **args), ...)
3839 {
3840 va_list ap;
3841 static char *arg[MAX_ARGS], **args = arg;
3842 char *arg1;
3843 long n, nargs = 0;
3844
3845 va_start(ap, routinePtr);
3846 if ((arg1 = va_arg (ap, char*)) != NULL) {
3847 arg[0] = arg1;
3848 while ((arg[++nargs] = va_arg(ap, char*)) != NULL) {};
3849 }
3850 va_end (ap);
3851
3852 routinePtr (plevel, args);
3853
3854 for (n = 0; n < plevel->iInstances; n++)
3855 TraverseLevels (plevel->pLevels[n], routinePtr, NULL);
3856
3857 } /* TraverseLevels */
3858
3859
3860 /* Function -------------------------------------------------------------------
3861 TraverseLevels1 (recursive)
3862
3863 Same as TraverseLevels, but checks the return code of the routine passed
3864 and returns the same code (0 if error, 1 if success).
3865 */
3866 int TraverseLevels1 (PLEVEL plevel,
3867 int (*routinePtr)(PLEVEL plevel, char **args), ...)
3868 {
3869 va_list ap;
3870 static char *arg[MAX_ARGS], **args = arg;
3871 char *arg1;
3872 long n, nargs = 0;
3873
3874 va_start (ap, routinePtr);
3875 if ((arg1 = va_arg (ap, char*)) != NULL) {
3876 arg[0] = arg1;
3877 while ((arg[++nargs] = va_arg(ap, char*)) != NULL) {};
3878 }
3879 va_end (ap);
3880
3881 if (routinePtr (plevel, args)) {
3882 for (n = 0; n < plevel->iInstances; n++) {
3883 if (!TraverseLevels1(plevel->pLevels[n], routinePtr, NULL)) {
3884 /* error */
3885 return (0);
3886 }
3887 }
3888 }
3889 else /* error */
3890 return (0);
3891
3892 /* success */
3893 return (1);
3894
3895 } /* TraverseLevels1 */
3896
3897
3898 /* Function -------------------------------------------------------------------
3899 WriteHeader
3900
3901 Write the complete output file header
3902 */
3903 void WriteHeader (PANALYSIS panal)
3904 {
3905 PGIBBSDATA pgd = &panal->gd;
3906 long i;
3907
3908 fprintf(pgd->pfileOut, "iter\t");
3909 TraverseLevels (panal->pLevels[0], WriteParameterNames, panal,
3910 pgd->pfileOut, NULL);
3911 if ((pgd->nSimTypeFlag == 3) || (pgd->nSimTypeFlag == 4)) {
3912 /* if the user has required tempering and specified a perk scale,
3913 then print header for temperature index and pseudo-priors etc. */
3914 fprintf(pgd->pfileOut, "IndexT\t");
3915 for (i = 0; i < pgd->nPerks; i++)
3916 fprintf(pgd->pfileOut, "LnPseudoPrior(%ld)\t",i+1);
3917 }
3918 fprintf(pgd->pfileOut, "LnPrior\tLnData\tLnPosterior\n");
3919 fflush(pgd->pfileOut);
3920
3921 } /* WriteHeader */
3922
3923
3924 /* Function -------------------------------------------------------------------
3925 WriteParameterNames
3926
3927 Called from Traverse Levels
3928 Write the names of the sampled parameters to output file header
3929 */
3930 void WriteParameterNames (PLEVEL plevel, char **args)
3931 {
3932 PANALYSIS panal = (PANALYSIS)args[0];
3933 PFILE outFile = (FILE*)args[1];
3934 long n, m;
3935
3936 panal->iInstance[plevel->iDepth] = plevel->iSequence;
3937
3938 for (n = 0; n < plevel->nMCVars; n++) {
3939 fprintf (outFile, "%s(", plevel->rgpMCVars[n]->pszName);
3940 for (m = 1; m < plevel->iDepth; m++)
3941 fprintf (outFile, "%d.", panal->iInstance[m]);
3942 fprintf (outFile, "%d)\t", panal->iInstance[plevel->iDepth]);
3943 }
3944
3945 } /* WriteParameterNames */
3946
3947
3948 /* Function -------------------------------------------------------------------
3949
3950 WriteMCVars
3951
3952 Write the values of MC vars for one level to output file
3953 */
3954 void WriteMCVars (PLEVEL plevel, char **args)
3955 {
3956 PFILE pOutFile = (PFILE)args[0];
3957 long n;
3958 PMCVAR pMCVar;
3959
3960 for (n = 0; n < plevel->nMCVars; n++) {
3961 pMCVar = plevel->rgpMCVars[n];
3962 fprintf(pOutFile, "%5g\t", pMCVar->dVal);
3963 }
3964
3965 } /* WriteMCVars */
3966
3967 /* End */
3968