aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--poiss.c41
1 files changed, 25 insertions, 16 deletions
diff --git a/poiss.c b/poiss.c
index a286184..dace05e 100644
--- a/poiss.c
+++ b/poiss.c
@@ -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"