analyzeMCMC.c File Reference

#include <TROOT.h>
#include <TFile.h>
#include <TTree.h>
#include <TH2D.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <math.h>
#include <iostream.h>

Go to the source code of this file.

Functions

int main ()


Function Documentation

int main (  ) 

Definition at line 11 of file analyzeMCMC.c.

00012 {
00013 
00014    // open file containing the MCMC tree 
00015 
00016    TFile * file = new TFile("./MCMC.root", "READ"); 
00017 
00018    // get the MCMC tree from the file 
00019 
00020    TTree * tree = (TTree*) file -> Get("MarkovChainTree"); 
00021 
00022    // get the number of parameters 
00023 
00024    int nparameters; 
00025 
00026    tree -> SetBranchAddress("fNParameters", &nparameters); 
00027 
00028    tree -> GetEntry(0); 
00029 
00030    // get probability density 
00031 
00032    double logpdf; 
00033 
00034    tree -> SetBranchAddress("fLogLikelihood", &logpdf); 
00035 
00036    // get the parameters 
00037    
00038    std::vector<double> parameter; 
00039 
00040    for (int i = 0; i < nparameters; ++i)
00041       parameter.push_back(0.0); 
00042 
00043    // loop over parameters 
00044 
00045    for (int i = 0; i < nparameters; ++i)
00046       tree -> SetBranchAddress(Form("fParameter%i", i), &(parameter.at(i))); 
00047                                                                
00048    // loop over iterations and calculate limits 
00049 
00050    int niterations = tree -> GetEntries(); 
00051 
00052    // debug
00053    niterations = 100001; 
00054 
00055    double minlogpdf; 
00056    double maxlogpdf; 
00057 
00058    double meanlogpdf = 0; 
00059 
00060    std::vector<double> minparameter; 
00061    std::vector<double> maxparameter; 
00062 
00063    std::vector<double> meanparameter; 
00064 
00065    for (int i = 0; i < nparameters; ++i)
00066       {
00067          minparameter.push_back(0.0); 
00068          maxparameter.push_back(0.0); 
00069          meanparameter.push_back(0.0); 
00070       }
00071 
00072    for (int iiteration = 0; iiteration < niterations; ++iiteration)
00073       {
00074          tree -> GetEntry(iiteration); 
00075          
00076          if (logpdf < minlogpdf || iiteration == 0) 
00077             minlogpdf = logpdf; 
00078 
00079          if (logpdf > maxlogpdf || iiteration == 0) 
00080             maxlogpdf = logpdf; 
00081          
00082          for (int i = 0; i < nparameters; ++i)
00083             {
00084                if (parameter.at(i) < minparameter.at(i) ||  iiteration == 0)
00085                   minparameter[i] = parameter.at(i); 
00086 
00087                if (parameter.at(i) > maxparameter.at(i) ||  iiteration == 0)
00088                   maxparameter[i] = parameter.at(i); 
00089 
00090                meanparameter[i] += parameter.at(i); 
00091             }
00092 
00093          meanlogpdf += logpdf; 
00094       }
00095 
00096    // calculate mean values 
00097 
00098    meanlogpdf = meanlogpdf / double(niterations); 
00099    
00100    for (int i = 0; i < nparameters; ++i)
00101       meanparameter[i] = meanparameter.at(i) / double(niterations); 
00102 
00103    // define histograms 
00104 
00105    double minlimit = minlogpdf; 
00106    double maxlimit = maxlogpdf; 
00107 
00108    if (minlimit < 0)
00109       minlimit = 1.1 * minlimit; 
00110 
00111    if (minlimit > 0)
00112       minlimit = 0.9 * minlimit; 
00113 
00114    if (maxlimit < 0)
00115       maxlimit = 0.9 * maxlimit; 
00116 
00117    if (maxlimit > 0)
00118       maxlimit = 1.1 * maxlimit; 
00119 
00120    TH2D * hist_pdf_vs_iteration = new TH2D("hist_pdf_vs_iteration", "", 
00121                                                                100, 0.0, double(niterations), 
00122                                                                100, minlimit, maxlimit); 
00123    hist_pdf_vs_iteration -> SetStats(kFALSE); 
00124    hist_pdf_vs_iteration -> SetXTitle("Iteration"); 
00125    hist_pdf_vs_iteration -> SetYTitle("log(pdf)");
00126 
00127    // loop over iterations 
00128 
00129    int order = 0; 
00130 
00131    for (int iiteration = 0; iiteration < niterations; ++iiteration)
00132       {
00133          tree -> GetEntry(iiteration); 
00134 
00135          // fill histograms 
00136 
00137          hist_pdf_vs_iteration -> Fill(iiteration, logpdf); 
00138       }
00139 
00140    // define autocorrelation functions 
00141 
00142    std::vector<double> aclogpdf; 
00143 
00144    int norders = int(floor(log10(niterations))); 
00145 
00146    TGraph * graph_aclogpdf = new TGraph(norders); 
00147    graph_aclogpdf -> SetMarkerStyle(20); 
00148 
00149    std::vector <TGraph *> graph_acparameter; 
00150 
00151    for (int i = 0; i < nparameters; ++i)
00152       {
00153          TGraph * g = new TGraph(); 
00154          g -> SetMarkerStyle(20+i+1); 
00155          graph_acparameter.push_back(g); 
00156       }
00157 
00158    std::vector <double> Aparameter; 
00159    std::vector <double> Bparameter; 
00160    std::vector <double> tempAparameter; 
00161 
00162    for (int i = 0; i < nparameters; ++i)
00163       {  
00164          Aparameter.push_back(0.0); 
00165          Bparameter.push_back(0.0); 
00166          tempAparameter.push_back(0.0); 
00167       }
00168 
00169    for (int iorder = 0; iorder < norders; ++iorder)
00170       {
00171          double Apdf = 0; 
00172          double Bpdf = 0; 
00173 
00174          int k = int(pow(10, double(iorder))); 
00175          
00176          for (int i = 0; i < niterations - k; ++i) 
00177             {
00178                tree -> GetEntry(i); 
00179 
00180                double tempApdf = logpdf - meanlogpdf; 
00181                
00182                for (int iparameter = 0; iparameter < nparameters; ++iparameter)
00183                   tempAparameter[iparameter] = parameter.at(iparameter) - meanparameter.at(iparameter); 
00184 
00185                tree -> GetEntry(i + k); 
00186                      
00187                Apdf += tempApdf * (logpdf - meanlogpdf); 
00188 
00189                for (int iparameter = 0; iparameter < nparameters; ++iparameter)
00190                   Aparameter[iparameter] += tempAparameter.at(iparameter) * 
00191                      (parameter.at(iparameter) - meanparameter.at(iparameter)); 
00192                
00193             }
00194 
00195          for (int i = 0; i < niterations; ++i) 
00196             {
00197                tree -> GetEntry(i); 
00198                
00199                Bpdf += ((logpdf - meanlogpdf) * 
00200                             (logpdf - meanlogpdf)); 
00201 
00202                for (int iparameter = 0; iparameter < nparameters; ++iparameter)
00203                   Bparameter[iparameter] += ((parameter.at(iparameter) - meanparameter.at(iparameter)) * 
00204                                                           (parameter.at(iparameter) - meanparameter.at(iparameter))); 
00205             }
00206          
00207          graph_aclogpdf -> SetPoint(iorder, double(iorder), Apdf/Bpdf); 
00208 
00209          for (int iparameter = 0; iparameter < nparameters; ++iparameter)
00210             {
00211                TGraph * g = graph_acparameter.at(iparameter); 
00212 
00213                g -> SetPoint(iorder, 
00214                                     double(iorder), 
00215                                     Aparameter.at(iparameter) / Bparameter.at(iparameter)); 
00216 
00217             }
00218       }
00219 
00220    // print to screen 
00221 
00222    TCanvas * canvas1 = new TCanvas("canvas1"); 
00223    canvas1 -> cd(); 
00224 
00225    hist_pdf_vs_iteration -> Draw("COLZ"); 
00226    
00227    canvas1 -> Print("canvas1.ps"); 
00228 
00229 
00230    TCanvas * canvas_aclogpdf = new TCanvas("canvas_aclogpdf"); 
00231    canvas_aclogpdf -> cd(); 
00232 
00233    TH2D * hist_acaxes = new TH2D("hist_acaxes", "", 
00234                                                 1, -0.5, double(norders)+0.5, 
00235                                                 1, -0.1, 1.1); 
00236    hist_acaxes -> SetStats(kFALSE); 
00237    hist_acaxes -> Draw(); 
00238 
00239    graph_aclogpdf -> Draw("SAMEPL"); 
00240 
00241    for (int iparameter = 0; iparameter < nparameters; ++iparameter)
00242       graph_acparameter.at(iparameter) -> Draw("SAMEPL"); 
00243 
00244    canvas_aclogpdf -> Print("canvas_aclogpdf.ps"); 
00245 
00246    // close root file 
00247 
00248    // file -> Close(); 
00249 
00250    return 0; 
00251 
00252 }


Generated on Thu Feb 11 11:29:30 2010 for BayesianAnalysisToolkit by  doxygen 1.5.1