diff options
| -rw-r--r-- | poiss.c | 41 |
1 files changed, 25 insertions, 16 deletions
@@ -8,6 +8,9 @@ void usage() { Uses actual Poisson probability (i.e. not Gaussians, and therefore not least squares). +Typical usage: +:; ./poiss -if my_data.csv -fit cu.csv + Options include: -h # print this message (and immediately exit) -ifile $fn # input file, date for -fit @@ -38,7 +41,9 @@ See ./doit for a working example, including the SS --meld stanzas. #include "poissonier.h" // Used for equal weighting: -vector<double> const u5(5, 1.0); +vector<double> uvec(int const ddd) { + return vector<double>(ddd, 1.0); +} // Class to describe a sitch, i.e. a situation, i.e. // everything the objective function needs to know: @@ -257,12 +262,14 @@ struct guesser { << endl; vector<double> slow_model{slow_r, slow_a, bl_r}; + int nprm = slow_model.size(); double fast_t0 = evt[0].t; double fast_n0 = evt[0].n; double fast_t1 = evt[fast_i].t; double fast_n1 = evt[fast_i].n; - double model_inty = inty_multi_exp(slow_model, fast_t0, fast_t1, u5); + double model_inty = inty_multi_exp(slow_model, + fast_t0, fast_t1, uvec(nprm)); double fast_tot = fast_n1 - fast_n0; double fast_dn = fast_tot - model_inty; double fast_r = fast_dn * fast_a; @@ -296,21 +303,20 @@ struct guesser { } ostream& xout(ofile != "-" ? xxout : cout); - size_t ss = sitch.t1.size(); - xout << ss << endl; - xout << sitch.t1[ss-1] - << " " << sitch.t2[ss-1] - << " " << sitch.obs[ss-1] + size_t nprm = sitch.t1.size(); + xout << nprm << endl; + xout << sitch.t1[nprm-1] + << " " << sitch.t2[nprm-1] + << " " << sitch.obs[nprm-1] << endl; - vector<double> prm(u5); + vector<double> prm(uvec(nprm)); xout << dump(model) << endl; xout << "bl/norm,bl,limp" << endl; double bl0 = prm[4]; void* context = (void*)(&sitch); - int nprm(prm.size()); for (double bl = bl0*.9; bl <= bl0*1.1; bl += bl/100.) { prm[4] = bl; @@ -377,9 +383,9 @@ void do_fit(poissonier& poi, string const ifile, ggg.more(); ggg.sitch.norm = ggg.model; ggg.show(gfile); + int nprm = ggg.model.size(); - vector<double> prm = u5; - int nprm(prm.size()); + vector<double> prm = uvec(nprm); vector<double> const lower(nprm, 0.5); vector<double> const upper(nprm, 1.5); @@ -402,6 +408,7 @@ void do_fit(poissonier& poi, string const ifile, } if (OK && ofile != "") { + string meta = "fast.amp, fast.dk, slow.amp, slow.dk, bl"; ofstream xxout; if (ofile != "-") { xxout.open(ofile); @@ -421,7 +428,7 @@ void do_fit(poissonier& poi, string const ifile, << ", perdof:, " << limp_end / DoF << endl; xout << endl; - xout << ",fast.amp, fast.dk, slow.amp, slow.dk, bl" << endl; + xout << "," << meta << endl; xout << "Normed:, " << dump(found) << endl; vector<double> combi(nprm); for (int ii = 0; ii < nprm; ii++) { @@ -437,12 +444,13 @@ void do_fit(poissonier& poi, string const ifile, // In particular, calculate the Mahalanobis metric // i.e. the second derivative (Hessian) of the log improbability. vector<vector<double> > covar(nprm, vector<double>(nprm)); - vector<double> probe(u5); + vector<double> probe(uvec(nprm)); sitcher sitch(ggg.sitch); sitch.norm = combi; void* context = (void*)(&sitch); double delta(0.01); xout << endl; + xout << meta << endl; for (int ii = 0; ii < nprm; ii++) { string sep = ""; for (int jj = 0; jj < nprm; jj++) { @@ -486,12 +494,13 @@ void do_dk(poissonier& poi, string const ofile) { double step(1); vector<double> prm( {30, .1, 10, .01, 1} ); + int nprm(prm.size()); int dt(5); for (int ii = 0; ii < 500; ii+= dt) { double tt1(step*ii); double tt2(tt1 + dt); - double intensity(inty_multi_exp(prm, tt1, tt2, u5)); + double intensity(inty_multi_exp(prm, tt1, tt2, uvec(nprm))); int obs = poi.sample(intensity); sitch.t1.push_back(tt1); sitch.t2.push_back(tt2); @@ -509,13 +518,13 @@ void do_dk(poissonier& poi, string const ofile) { // and the observations (to the extent they exist) for (int ii = 0; ii < 500; ii++) { double tt1(smooth.t1[ii]); - double sminty(inst_multi_exp(prm, tt1, u5)); + double sminty(inst_multi_exp(prm, tt1, uvec(nprm))); xout << tt1 << "," << sminty; if (ii < (int) sitch.t1.size()) { double fake(max(0.+sitch.obs[ii], 0.2)); // to facilitate log axes double tt1(sitch.t1[ii]); double tt2(sitch.t2[ii]); - double minty(inty_multi_exp(prm, tt1, tt2, u5)); + double minty(inty_multi_exp(prm, tt1, tt2, uvec(nprm))); xout << "," << tt1 << "," << tt2 << ",x" |