jtc_em.h - Header file classes for EM analysis of "Beads-in-a-Jar" experiments



/* File with methods to perform Expectation-minimisation and such like, inheriting the experiment stuff in jtc_exp.h */ #ifndef JTC_EM_H #define JTC_EM_H #include "jtc_exp.h" // For multiple integration : #include #include #include #include #include // for Gamma distro, incl. dummulative dist fns. : #include #include // for special functions, esp digamma, trigamma: #include /* D.a. ---------------------------- first simple C declarations ----------- */ // pdf of a Gamma distribution at q given its mean Gmean and its variance Gvar. // To rely on gsl necessities being already initialised in a ctor ... double gammaPDF(double q, double Gmean, double Gvar) ; // Basic classes assembling bits for MC integration: -------------------------- class MCIntegr { protected: Array parArr; int parNum; // will be from parArr size (i.e. rows() ...) // various bizarre bits that GSL random no generators etc. need: const gsl_rng_type *T; gsl_rng *r; int dim, calls; double cur_result, cur_error, cur_chisq; // outputs of gsl, incl. chi square check double (*integrand)(double *q, size_t dim, void *param); // fn. itself HAS TO be // given externally. public: MCIntegr(); // ~MCIntegr(); void SetUpMCI( double (*new_integrand)(double *q, size_t dim, void *param) , size_t new_dim , size_t new_calls = 100000 ) ; double vegasIntegrate( double *xl // calls the integrand with param = NULL ! , double *xu, int auxCalls = 10000 ) ; double miserIntegrate( double *xl // similar, for rougher work. , double *xu, int auxCalls = 10000 ) ; }; /* Classes : ------------------------------------------------------------- */ class Beads_Exp_II : public Beads_Exp_I // See jtc_exp.h { protected: static const char *const ERROR; int currModelType; // model type, 0 for control, 1 for the 'real', basic // Bayesian, one, 2 for SPRT,3 for basic Mackintosh_El-Deredy ... int maxMicroParNum ; // THIS IS 3, BECAUSE Bayesian-based model can have a max. // of 3 micro-params, CS,CW and T, also covering poissonoid // control (1). SOME METHODS RETURN ERROR IF IT'S NOT 3 ! int wrkDim; // if 1, we are only working with 1 param - can go up to maxMicroParNum // for Monte-Carlo integration methods float ApproxInf; // can ignore contribution of Gamma pdf over this many times its SD. // for storage ... : int gridSizePerDim ; int subjN; double min_curr_Par; // Don't start scannig from zero, quite ! Array AuxPV; // AuxPV: sequence_Num x Probability_vector e.g: 3 // (for Wellcome) x fN - usu. holds Py(d2d=drawNum) Array GdPV; // GdPV: sequence_Num x Probability_vector like // above but usu. holding corresp. Py(decide G) Array PVecSlice; // Another auxiliary, where to copy a slice of PVec. Array p_l; // yet another working auxiliary, will hold 2 rows of llr pdf int p_l_N; // copy of col. no. of above. Array p_l0; // separate working auxiliary, as above but for t=0 only. int p_l_N0; // copy of col. no. of above. double outL0, outL, outH0, outH, d0, d1, mu_hi, mu_lo; // auxiliaries for convolutions int Nlh0, Nlh, Nt0, Nt; // integrations and related .. Array PDumD_L; // will have range 1 to subject_number - for DumDat work. double sum_ln_PDumD_L; // (sum of ln of elements of above) public: // ctor initialises various bizarre bits that GSL random no generators etc. need, i.e. // const gsl_rng_type *T; gsl_rng *r; T = gsl_rng_default; r = gsl_rng_alloc (T); : Beads_Exp_II(); ~Beads_Exp_II(); MCIntegr mcInt; // Integration engine // copies of 'parameters' that will be needed for estimation of various probabili- // ty densities etc., esp. incl. recognition distro : // Data point currently considered, by participant number and draw type. It // has slightly different meaning in different types of expt., Welcome vs. ME1 // etc : float curr_pt, curr_gp, curr_draw_type; // integers, really, but float_ed to use w. menus. // 'macro' parameters currentl considered: for control poissonoid and SPRT // and for Bayesian distr.: double currPoidPar_m, currPoidPar_v ; double currHiTh_m, currTh_v ; // also to use currTm and currTv for 'noise' double currCSm, currCSv, currTm, currTv ; double currCGBm, currCGBv, currCBGm, currCBGv ; // Boundaries for integrations over hypercube : double *xl ; // lower limit NOTE THIS GOES FROM ZERO INDEX: ZERO *NOT* IGNORED double *xu ; // upper limit NOTE THIS GOES FROM ZERO INDEX: ZERO *NOT* IGNORED // Look-ups etc: int SubjN(){ return subjN; }; int GridSizePerDim(){ return gridSizePerDim; } ; int WDim(){ return wrkDim; }; double wrkP(int trial, int step_to_decl){ return AuxPV(trial, step_to_decl); }; double PdG(int trial, int step_to_decl){ return GdPV(trial, step_to_decl); }; double quantAuxPV(int seqN, int k, int q=1000); // quantile d2d for f=k/q from AuxPV double PdumD_L(int i){ return PDumD_L(i); } ; double sum_ln_PdumD_L(){ return sum_ln_PDumD_L; }; int currModType(){ return currModelType; }; void writeAuxPV(const char *fname); // Adjustments: void setCurrModType(double newModT){ currModelType=lround(newModT); }; // The following don't just copy data from array MS etc. into AuxPV, but it also // converts from P[Sample if in State] to P[Sample AND State-is-reached]. // First, for default, Welcome study: int StoreAuxPV( RBWelcomeDat *RBWptr, int validZero=0 ); // Second, for Makintosh - El-Deredy protocol. This is different // from StoreAuxPV above in that it only stores for ONE COST CONDITION // and needs to know which sequence was used for that cost. cond. It has // to use a fixed matrix MS of policy probabilities to sample again! It stores // in AuxPV ACCORDING TO COST CONDITION, NOT SEQUENCE IDENTITY AuxPV(cond,step) ! int StoreME0AuxPV( RBWelcomeDat *RBWptr, int cond, int seq_n, int validZero=0 ); // the above, generalized / emphasizing that it's not just for ME0 : int StoreME_AuxPV( RBWelcomeDat *RBWpt,int cnd,int sq_n,int valZero){ return StoreME0AuxPV( RBWpt, cnd, sq_n, valZero ); } ; double cProd3P_2(int d1, int d2, int d3, double L=-666.0) { return CPDistr1(d1,L)*CPDistr1(d2,L)*CPDistr1(d3,L); } ; // Various bits for dummy data (can be not so dummy at all, but this one, DumDat, // contains mainly decisions-to-draw data. Array DumDat ; // has ranges starting from 1. // fill DumDat with a ready-made 2D array from file 'fname'. Rem it is usually: // ptNum param1 data1 data2 data3 param2 ... int fillDumDat_w_blitz_file_arr(const char *fname); // The next two use currXX values to create artificial data & fill DumDat : void fillDumDatPoid(int subjN =30, int datN=3); // the following produces output formatted as // ptNum CS data1 data2 ... data_datN kT , with datN usu. = 3 // this does NOT allow for decision to be made at zero draws. void fillDumDatBays(RBWelcomeDat *RBWptr, int subjN=30, long int rndSeed = 0); // then, similar for SPRT. The first col. is CDPar1; Also produces // actual decisions as well, so no separate ...DumDecSPRT... needed: void fillDumDatSPRT(RBWelcomeDat *RBWptr, int subjN=30, long int rndSeed = 0); // Makintosh-El-Deredy MDP version, producing output formatted as // ptNum CS data1 data2 ... data_datN kT , with datN usu. = 3, // and also allowing choice at zero draws (or not!). // M0 is the 'dirty' version relying on CS. W.r.t. sequences used for // each cost level ME0 will usse 4:1:1, like full dataset. void fillDumME0DatBays(RBWelcomeDat *RBWptr, int subjN=60, int zeroDraw=0, long int rndSeed = 0); // ... and this is a cleaner one, again using 4:1:1 which is either the // full dataset or excluding 4:1:1 people w. d2d=20 (giving 36:9:9 pts) ! // Note also that it explicitly updates CBG = CGB = gsl_ran_gamma(r, K_CW, ThCW) // at crucial point. void fillDumME1DatBays(RBWelcomeDat *RBWptr, int subjN=54, int zeroDraw=1, long int rndSeed = 0); // ... and this is for ME2 version: void fillDumME2DatBays(RBWelcomeDat *RBWptr, int subjN=54, int zeroDraw=1, long int rndSeed = 0); // load one set of experimental results into DumDat (after appropr. resizing) // and record the group and draw_type (beads or words) in curr_gp & curr_draw_type: void loadExpToDumDat(RBWelcomeDat *RBWptr, int group, int exp_type=1); // ... and for Mackintosh - El-Deredy, first expt. 1 : void loadME1ExpToDumDat(RBWelcomeDat *RBWptr, int data_set); // SPRT / log-likelihood test related // First, zero but don't resize SPRT machinery. Not needed if fill_SPRT_P used. void SPRT_reset(){ p_l=0; AuxPV=0; GdPV=0; }; // Given the ptNo and trialNo (i.e. seqN), look up PG appropriately and // return the absolute magnitude of the log-likelihood ratio. Error-free version. double absLLRDumDatDecIfInState( RBWelcomeDat *RBWptr, int ptNo, int trialNo); // Propagate noisy SPRT llr pdf by one step - assumes that p_l contains values // updated at t-1. Uses kT as the sd of the noise, CDPar1 as upper threshold and // -rCW*CDPar1 as the lower threshold. validZero is whether decision before seeing // any draws is allowed. details_out is for check / debug output. void propNoisyLLRpdf( int t, int Seq, RBWelcomeDat *RBWptr, double *PS, double *PG=NULL, double *PB=NULL, int validZero=0, int details_out=0 ); // the next is a slow version to test against: void propNoisyLLRpdf0( int t, int Seq, RBWelcomeDat *RBWptr, double *PS, double *PG=NULL, double *PB=NULL, int validZero=0, int details_out=0 ); // This will fill AuxPV and GdPV for all 3 sequences using propNoisyLLRpdf, like // StoreAuxPV does but actually calculates probabilities, doesn't use ready-made: void fill_SPRT_P( RBWelcomeDat *RBWptr, int validZero=0 ); // The following is like StoreAllPVecs of base class, but uses fill_SPRT-like // code to store at the place correspoinding to the current 'iter' iteration: int StoreSPRT_PVecs( RBWelcomeDat *RBWptr, int validZero=0, int details=0 ); Array DumDec ; // As DumDat, but containing the decisions taken i.e. 1 or 0 // (1=Green/Good in Young & Bentall, Yellow in Welcome etc.) // fill the DumDec with generated decisions based on presumed-ready probability // matrices, given that the decision was taken in the step stored in DumDat: void fillDumDecBays(RBWelcomeDat *RBWptr, long int rndSeed = 0); // **NO** fillDumDecSPRT - it's incorporated in fillDumDatSPRT !! // Set ... : // set xl to zero and xu to respective mean+ApproxInf*SD acc. to the model type, // 0 for control, 1 for the 'real', basic Bayesian, one, 2 for SPRT, // 3 for basic Mackintosh_El-Deredy ... void SetApproxInf(int model_type = 0); // Set curr.. variables, then run SetApproxInf. model_types as above. In the case // of more complicated parameter structures, arrays can be passed with the pointers // e.g. *newCWx can be (CGBx, CBGx) : void SetCurrAndBounds( double *newTm, double *newTv, double *newV2m, double *newV2v ,double *newCWm=NULL, double *newCWv=NULL , double *newPParm=NULL ,double *newPParv=NULL, int new_model_type = 1 ); // reset no. of subjects & what depends on it: void ResetAndSubjN(int new_subjN, int datN=3){subjN=new_subjN; PDumD_L.resize(Range(1,new_subjN)); PDumD_L=0; DumDat.resize(Range(1,new_subjN),Range(1,6)); DumDat = 0; } ; void ResetPDumD_L_etc2(int new_subjN); // set subjN and rezize & zero PDumD_L int ResetPDumD_L_by_DumDat(); // set subjN and rezize & zero PDumD_L acc. to DumDat // and return the new value of subjN. // write to a disk file a set of artificial d2d data, created with poissoinoid micro // dynamics and Gamma macro statistics. void writeDummyPoidDat( char *fname, double Gmean, double Gvar, int subjN=30, int datN=3); // write to file plain gamma distro, using x coord relevant to dimension chosen void write1DGamma(char *fname, int chosen_dim, double mean, double var); // from xl to xu void write1DGamma2(char *fname, int chosen_dim, double mean, double var); // from 0 to xl/xu // a grid array for storing/working on m/pdf's. Assume no more than 3D param space !: Array pdfGrid ; // Clear the pdfGrid array and fill it with the current Prior (i.e. the gamma fns. // corresponding to currTm, currTv, curr CSm etc. writeCurrPriorSlice1 entrry is for // writing first 2D 'slice' of the 3D array, at i,j,1, to a file currPri.csv. Last // entry is for storing the whole 3D array - usually not !!! void fill_pdfGrid_curr(int writeCurrPriorSlice1 = 1, int verbose=1, int threeD=0); // in the next one, a very crude (hyper)volume estimate for the fn. stored in // pdfGrid. griddimUsed can be 1 to maxMicroParNum: double crudeGridFnVol(int gridDimUsed=1); // Set up the probability mass function on the grid for ME type expts (Grid is // kT x CGB), usu. w. 3 returns patterns and 3 sequences, acc. to // SoftME1_M_PG_V_QG_QB_QS_AS_MS etc. First, for model ME1: void set_ME_PMF_Grid(RBWelcomeDat *RBWptr, int zero_allowed, int verbose=1); // ... next, for model ME2: void set_ME2_PMF_Grid(RBWelcomeDat *RBWptr, int zero_allowed, int verbose=1); // Given the denom. of the recodistr., P[Data; macrop] add to each point of the grid // the appropr. recodistr. pdf value, divided by totPtNum if that is not zero, or by // DumDat.rows() if it is. Defaults also mean that it uses the data in the // DumDat arrray. Assumes data are in triplets-per-participant, like in RBWelcome. void addRecoDistPDF( int ptNum, double PD_L, int datDim=1 , int totPtNum=0, RBWelcomeDat *RBWptr = NULL , int useRBWdata = 0); // Developped from addRecoDist, used to output to disk a set of posterior, // 'experimental Bayesian' recognition distros ... // Uses DumDat, curr_draw_type, curr_group and PDumD_L(i) which must be ready! // if datDim is 1, it uses a constant CS (usu = currCSm) : void BayesExpDumDs( RBWelcomeDat *RBWptr , char *fNameStem, int verbose=1, int datDim=2); // Rough estimate, based on Simpson rule integration, of mean and var of the variable // asked for, calculating from pdfGrid. The variables signify which one(s) are to be // integrated over if only one or two are, acc. to wrkDim. If wrkDim is 3 (or more), // arguments are to be ignored // and integration to be over the whole grid. Defaults signify that if wrkDim // is 1, only look at CS; if it is 2, look at CS (2) and kT (1). void SimpsonGridStats( int variable, double *mean, double *variance,int var2D=1); // Give some idea of how good or bad the integration over the grid is - should // add up to 1. Rest of explanations similar to above. double SimpsonGridPCheck(int var1D=2, int var2D=1); // laborious stuff ! 1st, update PD_L for data in DumDat, and return the // updated sum_ln_PDumD_L. Modeltype=1 is Bayesian, 2 is SPRT: double update_PDumD_L( double (*new_PDF)(double *q, size_t dim, void *param), int mcCalls=10000, int model_type=1, int verbose= 1); // similar for constant CS (usu. =currCSm) using Simpson's rule: double update_PDumD_L0( double (*new_PDF)(double *q, size_t dim, void *param), int SimpsPtN, /*int model_type,*/ int verbose= 1); // Estimate (return new mean but don't update) currTm and currCSm : // RELY ON UPDATED PdumD_L ! double estCurrDumTm( double (*new_PDF)(double *q, size_t dim, void *param), int mcCalls=10000, int model_type=1, int verbose = 1); double estCurrDumCSm( double (*new_PDF)(double *q, size_t dim, void *param), int mcCalls=10000, int verbose = 1); double estCurrDumHiTh_m(double (*new_PDF)(double *q, size_t dim, void *param), int mcCalls=10000, int verbose = 1) ; // for SPRT double estCurrDumCGBm( double (*new_PDF)(double *q, size_t dim, void *param), int mcCalls=10000, int verbose = 1); // for ME1 + // Estimate (return vals. but don't update) on way to currTv and currCSv : // RELY ON UPDATED PdumD_L AND ON UPDATED currTm AND currCSm ! double estCurrDum_QlnT(double (*new_PDF)(double *q, size_t dim, void *param), int mcCalls=10000, int model_type=1, int verbose = 1); double estCurrDum_QlnCS(double (*new_PDF)(double *q, size_t dim, void *param), int mcCalls=10000, int verbose = 1); double estCurrDum_QlnTh(double (*new_PDF)(double *q, size_t dim, void *param), int mcCalls=10000, int verbose = 1); double estCurrDum_QlnCGB(double (*new_PDF)(double *q, size_t dim, void *param), int mcCalls=10000, int verbose = 1); // for ME1 + // THE FOLLOWING TWO ARE SUBOPTIMAL W.R.T. EM THEORY ! double estCurrDumTv( double (*new_PDF)(double *q, size_t dim, void *param), int mcCalls=10000, int verbose = 1); double estCurrDumCSv( double (*new_PDF)(double *q, size_t dim, void *param), int mcCalls=10000, int verbose = 1); // Interpolation machinery : // These two prepare PVecSlice and call setRegInp2D_x and setInp2D_y : void prepIntBasics(int mod_type); void prepPVecSlice(int d2d, int seq_num, int mod_type=1); // double wrkP2(int trial, int step_to_decl){ // Costed Bayesian version return IP2.PolynInt(kT,CS, PVecs(Range::all(),Range::all(),trial,step_to_decl)); }; double wrkP2s(int trial, int step_to_decl){ // s for SPRT version; note zero floor. double p= IP2.PolynInt(kT,CDPar1, PVecs(Range::all(),Range::all(),trial,step_to_decl)); return ((p > 0) ? p : 0 ); }; double wrkP3(int trial, int step_to_decl){ // 'basic' Costed Bayesian version w. CS return IP2.RationInt(kT,CS, PVecs(Range::all(),Range::all(),trial,step_to_decl)); }; // 'ME' type, e.g. Costed Bayesian version w. CGB==CBG; // ret_n is the returns pattern number, seq_n is the sequence used for that cost pattern. // see e.g. Beads_Exp_I::CalcMEPVecs in jtc_exp.cc double wrkP4(int ret_n, int seq_n, int step_to_decl){ return IP2.RationInt( kT,CGB, PVecs( Range::all(),Range::all() , seq_n+retNum*(ret_n-1), step_to_decl)); }; double wrkP5(int ret_n, int seq_n, int step_to_decl){ return IP2.RationInt_w_Poly_rescue( kT,CGB, PVecs( Range::all(),Range::all() , seq_n+retNum*(ret_n-1), step_to_decl)); }; // for MODEL ME2 - changed from wrkP* for clarity: double P_T_CS2(int ret_n, int seq_n, int step_to_decl); /* { return IP2.RationInt_w_Poly_rescue( kT,CS, PVecs( Range::all(),Range::all() , seq_n+retNum*(ret_n-1), step_to_decl)); }; */ // M phase related methods: // 1. Use averaged integration results, obtained separately and including // new estimate of mean of gamma distr., to derive new estimate of variance // of gamma distr. Param. NRiterNum is number of Newton-Raphson iterations // used in final phase of calculation; frLastCorr is a rough error estimate // (see code) : double xVarEst( double IntQx, double IntQlnx, int NRiterNum=10 , double *frLastCorr=NULL, int verbose = 1); }; // Don't forget the semilcolon at the end of class declarations ! // = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = #endif // JTC_EM_H // eof