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

// This is the energy routine, it needs to process the whole vector Xdata
//params is a simple struct containing 4 doubles.
double U_pot_GM(const params& theta, const double sig1, const double sig2, const double a1, const double a2,
const vector <double>& Xdata, const double sig0){
double U = 0;
for(int i=0; i<Xdata.size(); ++i){
U += log( likelihood_GM(theta, sig1, sig2, a1, a2, Xdata[i]) ); // sum over loglikelihoods
}
U = (theta.mu1*theta.mu1 + theta.mu2*theta.mu2)/(2*sig0*sig0) + log(2*PI*sig0*sig0); // logprior part.
return 1*U;
}
// This is the force routine, it needs to process B random elements of the vector Xdata
//forces and params are structs that contain 2 and 4 doubles, resp. .
forces get_noisy_force_GM(const params& theta, const double sig1, const double sig2, const double a1, const double a2,
vector <double>& Xdata, const size_t B, vector <int>& idx_arr, const double sig0){
forces F{0,0};
double P, e1, e2, x, scale;
scale = Xdata.size()/double(B);
int help_int, idx;
int size_minus_B = Xdata.size()B;
if(Xdata.size() != B){ // in case of subsampling...
// idx_arr stores the possible indices of the vector Xdata
// this loop randomly chooses B of them and stores them
// at the end of idx_arr.
for(int i=Xdata.size()1; i>=size_minus_B; i){
uniform_int_distribution<> distrib(0, i); // recreates this in every iter... is there a better way?
idx = distrib(twister);
help_int = idx_arr[i];
idx_arr[i] = idx_arr[idx];
idx_arr[idx] = help_int;
}
}
// actual force evaluation.
// the B data points to be considered are given
// by the B last indices stored in idx_arr.
for(int i=idx_arr.size()1; i>=size_minus_B; i){
x = Xdata[ idx_arr[i] ];
P = likelihood_GM(theta, sig1, sig2, a1, a2, x); // likelihood of a single data point
e1 = exp( 1*(xtheta.mu1)*(xtheta.mu1)/(2*sig1*sig1) );
e2 = exp( 1*(xtheta.mu2)*(xtheta.mu2)/(2*sig2*sig2) );
F.fmu1 += 1/P * e1 * (xtheta.mu1);
F.fmu2 += 1/P * e2 * (xtheta.mu2);
}
F.fmu1 *= a1/(sqrt(2*PI)*sig1*sig1*sig1) * scale;
F.fmu2 *= a2/(sqrt(2*PI)*sig2*sig2*sig2) * scale;
F.fmu1 = theta.mu1/(sig0*sig0); // prior part
F.fmu2 = theta.mu2/(sig0*sig0);
return F;
}
// This is the likelihood evaluation of a given data point x. It is used by both routines above.
double likelihood_GM(const params& theta, const double sig1, const double sig2, const double a1, const double a2, const double x){
double e1 = exp( 1*(xtheta.mu1)*(xtheta.mu1)/(2*sig1*sig1) );
double e2 = exp( 1*(xtheta.mu2)*(xtheta.mu2)/(2*sig2*sig2) );
double p = a1/(sqrt(2*PI)*sig1) * e1 + a2/(sqrt(2*PI)*sig2) * e2;
return p;
}
