#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 () |
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 }