C++ Parm::change方法代码示例

void Mcmc::runChain(void) {

	// get initial likelihood
	modelPtr->updateAllFlags();
	double oldLnL = modelPtr->lnLikelihood();
	
	for (int n=1; n<=numCycles; n++)
		{
		// pick a parameter to change
		Parm* parm = modelPtr->pickParmAtRandom();
		
		// propose a new state for the parameter
		double lnProposalRatio = parm->change();
		
		// calculate the prior ratio 
		double lnPriorRatio = parm->lnPriorRatio();
		
		// calculate the likelihood
		modelPtr->updateQ(parm);
		double lnL = modelPtr->lnLikelihood();
		double lnLikelihoodRatio = lnL - oldLnL;
		
		// calculate the acceptance probability of the proposed state
		double r = safeExp( lnLikelihoodRatio + lnPriorRatio + lnProposalRatio );
		
		// accept or reject the proposed state as the next state of the chain
		bool acceptState = false;
		if (ranPtr->uniformRv() < r)
			acceptState = true;
			
		// print information to the screen
		if ( n % printFrequency == 0 )
			{
			std::cout << std::setw(5) << n << " -- ";
			std::cout << std::fixed << std::setprecision(8) << oldLnL << " -> " << lnL;
			if (acceptState == true)
				std::cout << " -- Accepted change to " << parm->getName() << std::endl;
			else
				std::cout << " -- Rejected change to " << parm->getName() << std::endl;
			}
		
		// update the state of the chain
		if (acceptState == true)
			{
			parm->incrementNumAccepted();
			oldLnL = lnL;
			modelPtr->keep(parm);
			}
		else 
			{
			modelPtr->restore(parm);
			modelPtr->restoreQ(parm);
			}

		// print out current chain state
		if ( n == 1 || n % sampleFrequency == 0 || n == numCycles )
			printChainState(n, oldLnL);
			
		// print summary of chain
		if ( n % summarizeFrequency == 0 )
			treeSummary->printSummary();
		}
		
	treeSummary->print();
	modelPtr->printAcceptanceInfo();
}

Mcmc::Mcmc(Model *m, MbRandom *r, Settings *s, IoManager *outputPtr) {

	/* initialize some pointers and variables */
	modelPtr    = m;
	rnd         = r;
	chainLength = s->getChainLength();
	printFreq   = s->getPrintFrequency();
	sampleFreq  = s->getSampleFrequency();

	/* open files for output */
	string fileName = "";
	if (outputPtr->getFileName() == "")
		fileName = "mcmc_out";
	else
		fileName = outputPtr->getFileName();
	string filePathAndName = outputPtr->getFilePath() + PATH_SEPERATOR + fileName;
	string parmFile = filePathAndName + ".parm";
	string treeFile = filePathAndName + ".tree";
	ofstream parmOut( parmFile.c_str(), ios::out );
	ofstream treeOut( treeFile.c_str(), ios::out );
	if (!parmOut) 
		{
		cerr << "Cannot open file \"" + parmFile + "\"" << endl;
		exit(1);
		}
	if (!treeOut) 
		{
		cerr << "Cannot open file \"" + treeFile + "\"" << endl;
		exit(1);
		}
	
	cout << "   Running Markov chain" << endl;
	cout << "      Number of cycles      = " << chainLength << endl;
	cout << "      Print frequency       = " << printFreq << endl;
	cout << "      Sample frequency      = " << sampleFreq << endl;
	cout << "      Parameter file name   = " << fileName + ".parm" << endl;
	cout << "      Tree file name        = " << fileName + ".tree" << endl << endl;
	cout << "   Chain" << endl;
	
	/* run the Markov chain */
	int timeOn = clock();
	for (int n=1; n<=chainLength; n++)
		{
		/* pick a parameter to change */
		int whichParm = modelPtr->pickParm();
		Parm *parmOld = modelPtr->getCurParameter(whichParm);
		modelPtr->flipActiveParm();
		Parm *parmNew = modelPtr->getCurParameter(whichParm);

		/* change the parameter */
		double lnProposalRatio = parmNew->change();
		
		/* calculate likelihood for new state */
		modelPtr->upDateQ(parmNew);
		modelPtr->upDateTiProbs();
		double lnLOld = modelPtr->getLnL();
		double lnLNew = modelPtr->lnLikelihood();

		if (n % printFreq == 0)
			cout << setw(8) << n << " -- " << fixed << setprecision(2) << lnLOld << " -> " << fixed << setprecision(2) << lnLNew;
		
		/* calculate the acceptance probability */
		double lnLikelihoodRatio = lnLNew - lnLOld;
		double lnPriorRatio = parmNew->getLnPriorProbability() - parmOld->getLnPriorProbability();
		double lnR = lnLikelihoodRatio + lnPriorRatio + lnProposalRatio;
		double r = 0.0;
		if (lnR > 0.0)
			r = 1.0;
		else if (lnR < -100.0)
			r = 0.0;
		else
			r = exp(lnR);
		
		/* accept or reject the move */
		bool acceptMove = false;
		if (rnd->uniformRv() < r)
			acceptMove = true;
			
		/* update state of chain */
		if (acceptMove == true)
			{
			modelPtr->setLnL(lnLNew);
			int from = modelPtr->getActiveParm();
			int to   = flip(from);
			modelPtr->copyParms( from, to );
			if (n % printFreq == 0)
				cout << " -- Accepted change to " << parmNew->getParmName() << endl;
			}
		else
			{
			modelPtr->restoreQ(parmNew);
			int to = modelPtr->getActiveParm();
			int from   = flip(to);
			modelPtr->copyParms( from, to );
			modelPtr->flipActiveParm();
			if (n % printFreq == 0)
				cout << " -- Rejected change to " << parmNew->getParmName() << endl;
			}
			
		/* print state of the chain to a file */
//.........这里部分代码省略.........