56 REAL8 h0range = 0, h0step = 0;
57 INT4 h0steps = 0,
i = 0;
60 ProcessParamsTable *
ppt;
61 REAL8 scale = 1., scalemin = 0., tmpscale = 0., tmpmin = 0., tmpgridval = 0.;
63 ProcessParamsTable *commandLine = runState->
commandLine;
71 CHAR *parname = NULL, parscale[256], parmin[256], outputgrid[256];
77 ProcessParamsTable *ppt2;
85 sprintf( parscale,
"%s_scale", parname );
86 sprintf( parmin,
"%s_scale_min", parname );
95 h0min = atof( ppt2->value );
103 h0max = atof( ppt2->value );
111 h0steps = atoi( ppt2->value );
120 fprintf( stderr,
"Calculating posterior on %s over a grid from:\n", parname );
121 fprintf( stderr,
"\t%le --> %le in %d steps.\n", h0min, h0max, h0steps );
124 h0range = h0max - h0min;
125 h0step = h0range / (
REAL8 )( h0steps - 1. );
144 sprintf( outputgrid,
"%s_grid_posterior.txt", parname );
146 if ( (
fp = fopen( outputgrid,
"w" ) ) == NULL ) {
147 fprintf( stderr,
"Error... cannot open grid posterior file %s.\n",
152 for (
i = 0;
i < h0steps;
i++ ) {
153 REAL8 h0val = h0min +
i * h0step;
160 if ( logL->
data[
i] < minL ) {
161 minL = logL->
data[
i];
166 for (
i = 0;
i < h0steps - 1;
i++ ) {
167 sumPost += ( exp( logL->
data[
i] - minL ) + exp( logL->
data[
i + 1] - minL ) ) *
172 for (
i = 0;
i < h0steps;
i++ ) {
173 REAL8 h0val = h0min +
i * h0step;
174 fprintf(
fp,
"%le\t%le\n", h0val, exp( logL->
data[
i] - minL ) / sumPost );
219 fprintf( stderr,
"Error... testing Gaussian likelihood required the \"H0\" parameter to be set in the prior file.\n" );
226 loglike = -
log( sqrt( 2.*
LAL_PI ) * like_sigma );
227 loglike -= 0.5 * (
h0 - like_mean ) * (
h0 - like_mean ) / ( like_sigma * like_sigma );
230 data->likeli_counter += 1;
247 ProcessParamsTable *
ppt = NULL;
255 fprintf( stderr,
"Error... no prior range set for the test Gaussian likelihood" );
260 REAL8 h0mean = 0., h0sigma = 0., h0sigma2 = 0., h095 = 0.;
263 h0sigma2 = h0sigma * h0sigma;
274 REAL8 p_Z = exp( lnC - Z );
277 REAL8 G = ( 1. / ( sqrt(
LAL_TWOPI ) * h0sigma ) ) * ( (
min - h0mean ) * exp( -0.5 * (
min - h0mean ) * (
min - h0mean ) / h0sigma2 ) - (
max - h0mean ) * exp( -0.5 * (
max - h0mean ) * (
max - h0mean ) / h0sigma2 ) );
278 REAL8 KLdiv = -0.25 * p_Z * (
D + 2.*G );
292 if ( (
fp = fopen(
outfile,
"w" ) ) != NULL ) {
293 fprintf(
fp,
"%.12le\t%.12le\t%.12le\n", Z, h095, KLdiv );
296 fprintf( stderr,
"Warning... could not open test Gaussian output file '%s'.",
outfile );
301typedef struct tagul_params {
314 double sigma =
p->sigma;
316 double area =
p->area;
318 double C = 1. / ( sqrt( 2. ) * sigma );
320 return ( 0.5 * ( erf(
C * (
mu -
min ) ) - erf(
C * (
mu -
x ) ) ) / area ) -
ul;
331 int iter = 0, max_iter = 100;
332 const gsl_root_fsolver_type *
T;
335 double epsrel = 1
e-4;
337 double C = 1. / ( sqrt( 2. ) * sigma );
338 double area = 0.5 * ( erf(
C * (
mu -
min ) ) - erf(
C * (
mu -
max ) ) );
345 x_lo =
mu - 10.*sigma;
346 x_hi =
mu + 10.*sigma;
351 T = gsl_root_fsolver_brent;
352 s = gsl_root_fsolver_alloc(
T );
353 gsl_root_fsolver_set(
s, &F, x_lo, x_hi );
357 gslstatus = gsl_root_fsolver_iterate(
s );
358 x = gsl_root_fsolver_root(
s );
359 x_lo = gsl_root_fsolver_x_lower(
s );
360 x_hi = gsl_root_fsolver_x_upper(
s );
363 gslstatus = gsl_root_test_interval( x_lo, x_hi, 0., epsrel );
364 }
while ( gslstatus == GSL_CONTINUE && iter < max_iter );
366 if ( gslstatus != GSL_SUCCESS ) {
371 gsl_root_fsolver_free(
s );
398 if ( (
fp = fopen(
"prior_samples.txt",
"w" ) ) == NULL ) {
399 fprintf( stderr,
"Error... could not open prior samples file\n" );
452 ProcessParamsTable *
ppt = NULL;
453 CHAR *sampfile = NULL, *namefile = NULL;
454 CHAR *parambuf = NULL, *databuf = NULL;
458 REAL8 logLnew = 0., logL = 0.;
467 namefile =
ppt->value;
481 sampfile =
ppt->value;
492 fp = fopen(
"likelihoodComp.txt",
"w" );
519 fprintf( stderr,
"%.16le\t%.16le\t%.16le\n", logL, logLnew, logL - logLnew );
520 fprintf(
fp,
"%.16le\t%.16le\t%.16le\n", logL, logLnew, logL - logLnew );
#define __func__
log an I/O error, i.e.
const double scale
multiplicative scaling factor of the coordinate
char * XLALFileLoad(const char *path)
void LALInferenceAddVariable(LALInferenceVariables *vars, const char *name, const void *value, LALInferenceVariableType type, LALInferenceParamVaryType vary)
REAL8 LALInferenceGetREAL8Variable(LALInferenceVariables *vars, const char *name)
void LALInferenceRemoveVariable(LALInferenceVariables *vars, const char *name)
void LALInferenceCopyVariables(LALInferenceVariables *origin, LALInferenceVariables *target)
void * LALInferenceGetVariable(const LALInferenceVariables *vars, const char *name)
ProcessParamsTable * LALInferenceGetProcParamVal(ProcessParamsTable *procparams, const char *name)
int LALInferenceCheckVariable(LALInferenceVariables *vars, const char *name)
void LALInferenceSetVariable(LALInferenceVariables *vars, const char *name, const void *value)
LALINFERENCE_PARAM_CIRCULAR
LALINFERENCE_PARAM_LINEAR
int LALInferenceCheckMinMaxPrior(LALInferenceVariables *priorArgs, const char *name)
void LALInferenceGetMinMaxPrior(LALInferenceVariables *priorArgs, const char *name, REAL8 *min, REAL8 *max)
void * XLALCalloc(size_t m, size_t n)
char char * XLALStringDuplicate(const char *s)
int char * XLALStringAppend(char *s, const char *append)
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
#define XLAL_ERROR_VOID(...)
#define XLAL_ERROR_REAL8(...)
int XLALPrintError(const char *fmt,...) _LAL_GCC_PRINTF_FORMAT_(1
p
RUN ANALYSIS SCRIPT ###.
void test_gaussian_output(LALInferenceRunState *runState)
Output the analytic evidence for the test Gaussian likelihood and a 95% upper limit.
REAL8 test_gaussian_log_likelihood(LALInferenceVariables *vars, LALInferenceIFOData *data, LALInferenceModel *get_model)
Test the sampler using a Gaussian likelihood.
void compare_likelihoods(LALInferenceRunState *rs)
Read in an ascii text file of nested samples and compare the log likelihoods.
void outputPriorSamples(LALInferenceRunState *runState)
Output a number of prior samples based on the initial live points.
static double ul_gauss_cdf_function(double x, void *params)
void gridOutput(LALInferenceRunState *runState)
A test function to calculate a 1D posterior on a grid.
static double ul_gauss_CDFRoot(double mu, double sigma, double min, double max)
Find the root of the Gaussian CDF to give a 95% upper limit.
Header file for functions used in testing the parameter estimation code for known pulsar searches usi...
#define USAGEGRID
The usage format for the test case of performing the analysis on a one-dimensional grid.
LALInferenceVariables * params
REAL8 * ifo_loglikelihoods
LALInferenceIFOModel * ifo
ProcessParamsTable * commandLine
LALInferenceVariables * algorithmParams
struct tagLALInferenceIFOData * data
LALInferenceVariables * priorArgs
LALInferenceThreadState * threads
LALInferenceVariables ** livePoints
LALInferenceLikelihoodFunction likelihood
LALInferenceVariables * currentParams
LALInferenceModel * model
struct tagVariableItem * next
LALInferenceParamVaryType vary
LALInferenceVariableItem * head