timeseries.h

最大似然估计算法

#include <math.h>
#include <stdio.h> 
#include <stdlib.h> 
#include <string.h>
#include <pwd.h>
#include <sys/types.h>
#include <sys/utsname.h>
#include <time.h> 
#include <unistd.h>
#include "f2c.h"
#include "blaswrap.h"

typedef struct {
	int n_data;
	int n_full;
	int n_series;
	int n_offsets;
	int n_psd;

	int got_errors;
	int current_series;
	int count;

	int *index;
	int *off_code;
	int *c_id;

	double fs;
	double time_span;
	double *f;
	double *P;
	double *t;
	double *data;
	double *formal_error;
	double *offsets;

} time_series;

typedef struct {
	int n_par;
	int n_data;

	int *name_id;   /* id for name of parameter           */
	int *alt_id;    /* id for alternate name of parameter */

	double *MLE;

	double *A;      /* data kernel         */
	double *d;      /* data vector         */
	double *dd;     /* formal error vector */

	double *params; /* estimated parameter vector     */
	double *covar;  /* estimated parameter covariance */
	
} data_kernel;

typedef struct {
                               /************************************************/
	char   model;          /* single character code for model              */
	int    c_id;           /* id for series names                          */
                               /*                                              */
	int    sigma_id;       /* id for name of sigma e.g. pl, wh, fogm       */
	int    *sigma_flag;    /* flag : 0 free, 1 fixed, 2 estimated          */
	double *sigma;         /* vector of sigmas                             */
	double *dsigma;        /* vector of sigma uncertainties                */ 
                               /*                                              */
	int     n_pvec;        /* size of parameter vector                     */
	int    *names;         /* vector of index for names of pars.           */
	int    *pvec_flag;     /* flag : 0 free, 1 fixed, 2 estimated          */
	double *pvec;          /* parameter vectors                            */
	double *pvec_sigma;    /* vector of pvec uncertainties                 */
                               /************************************************/

} noise_model;

typedef struct {
	int n_data;
	double *C;
} covariance;

typedef struct {

	int est_method;
	int cov_method;
	int verbose;
	int speed;
	int input_method;

	double *tol;
	double delta;

	FILE *fpout;
	FILE *fpsd;
	FILE *fpkernel;
} options;


static char *par_names[] = { 
	"INTER  : ", "SLOPE  : ", "SIN    : ", "COS    : ", 
	"OFFSET : ", "WH     : ", "PL     : ", "FOGM   : ", 
	"BETA   : ", "BP     : ", "INDEX  : ", "FCEN   : ",
	"WIDTH  : ", "NPOLE  : ", "VWH_S  : ", "BTIME  : ",
	"GGM    : ",
	"Spectral_index",
	"Beta",
	"Fcentral",
	"width",
	"number_poles",
	"vwh_scale",
	"begin_time",
	"VWH_T  : ",
	"vwh_time",
	"ETIME  : ",
	"end_time",
	"NULL"};

static char *comp_names[] = {"NORT", "EAST", "VERT","RLR ","PSMSL"};
static char *version = {"3.1.2"};

static double sec_per_year = 31556925.216;

void        mle_wrapper(noise_model *, int, time_series, data_kernel, options);
void        emp_wrapper(noise_model *, int, time_series, data_kernel, options);
void        wls_wrapper(noise_model *, int, time_series, data_kernel, options);
void        psd_wrapper(noise_model *, int, time_series, data_kernel, options, double *, double *, int);
time_series read_cats_series(char *, options, double, int);
time_series read_rlrdata(char *, options, double);
noise_model decode_model_string(char *, int);
void        print_model(noise_model, options, double);
data_kernel create_A_and_d(time_series, double *, int, int, int);
int         sort_offsets(time_series, int, int *);
int         create_spectrum(double *, double *, time_series, data_kernel, options);
double      white_only(time_series, data_kernel, noise_model, options, int);
double      color_only(time_series, data_kernel, noise_model, options, int);
double      brent_color_only(time_series, data_kernel, noise_model, options);
double      simplex_color_only(time_series, data_kernel, noise_model, options, int);
double      brent_color_white(time_series, data_kernel, noise_model *, int, options, int);
double      general_inner_simplex(time_series, data_kernel, noise_model *, int, options, int);
double      general_outer_simplex(time_series, data_kernel, noise_model *, int, options);
void        cats_empirical(time_series, data_kernel, noise_model, options, double *, double *, double *);
void        create_covariance(noise_model, time_series, int, double *, int);
void        fractional_diff(double *, int, double);
void        generate_FOGM_cov(double *, int, double, double, double, double *);
void        power_law_cov(double *, int, double, double, double, double *);
double      general_MLE(data_kernel, double *, options, double *, double *, double *, int);
double      whittle_inner_simplex(double *, double *, int, double, data_kernel, noise_model *, int, options, int);
double      whittle_outer_simplex(double *, double *, int, time_series, data_kernel, noise_model *, int, options);
void        band_pass_tv(double *, int, double, double, int);
void        calculate_fisher(double **, data_kernel, noise_model *, int, double *, double *);
double      calculate_fs(double *, int, int *, int *);
double      calculate_fs2(double *, int, int *, int *);
double      color_only_mle(double *, int, double *, double *, double *, int);
double      color_white_mle(double *, int, double *, double *, double *);
double      color_white_mle_res(double *, int, int, int, double *, double *, double *, double *);
void        cov_scale(double *, int, double *, double *);
void        create_power_spectrum(noise_model, double *, double, int, double *, int);
double      est_line_WH_only(data_kernel, options,  double *, double *);
double      exact_radius(int, double *, double *, double *, double, double *);
void        fft1d(double *, int, int);
void        fogm_tv(double *, int, double, double);
double      linefit_all_fast(double *, double *, double *, int, int, double *, double *, double *);
double      linefit_all_fast3(data_kernel, double *, int, double *, double *, double *, options);
void        linefit_fast(double *, double *, double *, int, int, double *, double *, double *);
double      MLE_withline_WH(double, int, double *);
void        power_law_tv(double *, int, double);
void        qr(double *, int, int, double *, double *); 
void        scargle(double *, double *, int, double, double, double *, double *, int, int *, int *, double *);
void        sigma_from_rphi(int, double *, double *);
void        sort_params_to_models(noise_model *, int, double *, int, int);
double      whittle_MLE(double *, double *, int, double *, int);
void        avevar(double *, int, double *, double *);
int         comp_double_asc(const void *p_1, const void *p_2);
double      MLE_fixed_param(double **, double *, int, data_kernel);
void        realft(double *, int , int);
void        four1(double *, int, int);
void        psd(double *, double *, double *, int, double);