Revision 49cb69c5a18fdb262964fbfeb47ab2099eb32c5c authored by Wesley Tansey on 03 May 2018, 19:46:59 UTC, committed by Wesley Tansey on 03 May 2018, 19:46:59 UTC
braindt.h
#ifndef briandtH
#define braindtH
#include<iostream>
#include <fstream>
#include<vector>
#include"randomc.h"
#include"stocc.h"
#include<cstring>
using namespace std;
class braindt{
public:
braindt(int &LX, int &LY, int &LZ,int &Num,char *argv ,vector<int> &Loca,int &PN,vector<double>& est,int &L_hat,int sweep_b, int sweep_r,int sweep_lis,int iter_max);//constructor function of class braindt
double normpdf(double &x, double &mu, double &sigma_sq);// The probability density function of the normal distribution with mean=mu, variance=sigma_sq
bool IsFiniteNumber(double &x);// Judge whether x is a finite numebr
bool IsNumber(double &x);//Judge whether x is a numebr
double sigmasum(vector<double> &p, int begin, int end);//sum of elements of vector p
double sigmasum(vector<double> &p1, vector<double> &p2, int begin, int end);//sum of elements of the entrywise-product of vectors p1 and p2
double sigmasum(vector<double> &p1, vector<double> &p2, double p2minus, int begin, int end);//sum of elements of the entrywise-product of vectors p1 and p3 squared, where p3=p2-p2minus
void two_d_var_matrix_inv(ofstream &fileout,int splnumber);//compute the inverse matrix of a 2 by 2 matrix
void two_d_var_matrix_inv(int splnumber);
void gem(ofstream &fileout,char* filename);//The Generalized EM Algorithm
void gem(char* filename);
void backtrack(ofstream &fileout, StochasticLib1 &rnd2);//The backtracking line search method
void backtrack(StochasticLib1 &rnd2);
void func(ofstream &fileout,StochasticLib1 &rnd2);//update U and q2 in the Armijo condition
void func(StochasticLib1 &rnd2);
void computeLIS(char* filename);//computing LIS
private:
static const double TINY;//a very small number to avoid the zero denominator of a fraction
int Lx,Ly,Lz,N,YbyX;//N=Lx*Ly*Lz is the size of the image grid, and YbyX=Ly*Lx
vector<double> y;//observed data
vector<int>::iterator loca;// the location of the brain region of interest (ROI) on the Lx*Ly*Lz image grid
int pN;//the number of voxels in the ROI, i.e., the size of loca
///////////////////////////////////////////
int L;// the number of components in the normal mixture of non-null distribtion
//parameters in the null distribution; we assume the null is standard normal.
static const double mu0;//
static const double sigma0_sq;//
//parameters in the non-null distribution
vector<double> mu1_hat;
vector<double> sigma1_sq_hat;
vector<double> pEll_hat;//the weight for each component
//parameters in the Ising model
double beta_hat;//parameter for sum of x_s*x_t,where x is the unobservable state
double h_hat;//parameter for sum of x_s
///////////////////////////////////////////
vector<double> U;// the score function
vector<double> I;// the information matrix
vector<double> invI;// the inverse information matrix
int swp_b, swp_r,max_iter;//swp_b is the number of burn-in iterations, swp_r is samples from the Gibbs sampler, max_iter is the maximum limit of iterations in GEM.
vector<double> delta;//the regular Newton step
double beta_hat_old, h_hat_old;//previous values for m iteration, i.e., the interations in backtracking method
vector<int> initial;//the initial random field for the gibbs sampler
vector<double> H_cond_mean;//mean H conditional on observed data
vector<double> H_mean0;//Used to reset H_cond_mean or H_mean to be zeros
vector<vector<double> > H0;//Used to reset H to be zeros
vector<double> H_mean;// mean H
vector<vector<double> > H;//H
////////////////////
double log_plikelihood;// sum of exp(-phi^T*H(sample[i])), in the difference in Q_2
//previous values for t iteration, i.e., the interations in GEM
vector<double> mu1_hat_pre;
vector<double> sigma1_sq_hat_pre;
vector<double> pEll_hat_pre;
double alam;// the lambda in the Armijo condition
int stpover;//indicate whether over the max step
double q2,q2_old, delta_q2;// q_2 is log(sum of exp(-phi^T*H(sample[i]))), in the in the difference in Q_2
////////////////////////////////////
vector<double> LIS;
int swp_lis;// the gibbs sample-size for computing LIS
};
#endif
Computing file changes ...