forked from basvanopheusden/BMD
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbmd.h
142 lines (133 loc) · 3.71 KB
/
bmd.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
#ifndef CHANGEPOINTS_H_INCLUDED
#define CHANGEPOINTS_H_INCLUDED
#include <vector>
#include <random>
#include <cmath>
#include <fstream>
#include <string>
using namespace std;
struct integraltable{
double* table;
double* a_upper;
double* a_lower;
int Nd;
int Na;
double d_step_size;
integraltable(char*);
~integraltable();
int get_a_ind(const double);
int get_d_ind(const double);
double get_a(const int);
double get_d(const int);
double operator()(const double,const double);
};
struct params{
double lambda0=0.004, lambda1=0.1; // parameters of the prior distributions (gamma) over durations
long double d0=1.0;
long double d1,sigma0,sigma1,sigmaz,sigmax;
double const_term_up;
double const_term_down;
double prefactor;
mt19937_64 engine;
params();
params(char*);
params(const double, const double, const double, const double, const double);
void write(char*);
void write();
void set_d_sigma_up(const double, const double);
void set_sigma_down(const double);
void set_sigmaz(const double);
void set_sigmax(const double);
};
struct changepoints;
struct data{
double* x[2];
double* z[2];
double* zf[2];
int T;
double sum_dz_squared;
data(char*);
~data();
void write_z(char*);
void kalman_filter(const params& );
double get_dz_squared(int,int);
void get_residual(changepoints&, double*[2]);
};
struct changepoints{
vector<int> t01;
vector<int> t10;
int n,N;
vector<bool> get_vec();
changepoints(int);
changepoints(data&);
changepoints(char*,int);
void write(char*);
void write();
void execute_step(int,int,int); //makes a step with a particular type, changepoint index k and position
};
struct MCMC_chain{
double beta;
data* dat;
params theta;
changepoints C;
integraltable* table;
double logposterior;
double logprior;
double loglik;
double logpostbest;
int nsteps;
int swaps_acc;
int swaps_tried;
mt19937_64 engine;
discrete_distribution<int> get_step_type;
MCMC_chain(double,data*,params&,changepoints&,integraltable*);
bool make_step(const params&); //generates a step type and k, accepts with Metropolis probability
//and returns true if successful
bool accept(double);//implements the metropolis acceptance probability
bool swap_chains(MCMC_chain*);
double get_logprior();
double get_logprior(const int,const int);
double loglik_up(const int,const int);
double loglik_down(const int,const int);
double get_loglik();
double get_logpost();
void change_params(const params&);
};
class parallel_tempered_chain{
public:
//static const int Nchains=82;
static const int Nchains=1;
MCMC_chain* chain[Nchains];
uniform_int_distribution<int> random_chain;
//static const int neo=55;
static const int neo=0;
mt19937_64 engine;
data dat;
params theta;
//changepoints C_real; //if we want the Creal input file
params theta_real;
integraltable table;
double logpostreal;
vector<changepoints> C_samples;
int nsamples;
//parallel_tempered_chain(double*,char*,char*,char*); //if we want the Creal input file
parallel_tempered_chain(double*,char*,char*);
~parallel_tempered_chain();
double get_loglik(params&,changepoints&);
double get_logprior(params&,changepoints&);
double get_logpost(params&,changepoints&);
void update_params();
void sweep();
void sweep(int);
void estimate_allthethings(char*,char*);
void sample_changepoints(int);
void initialize_changepoints();
void estimate_z();
void estimate_sigmazx();
void estimate_dsigma_up();
void estimate_sigma_down();
void write_posteriors(int,char*, char*);
void write_changepoints(char*);
void write_samples(char*);
};
#endif // CHANGEPOINTS_H_INCLUDED