em_rb_beads07.cc - 'main' file for EM fitting of Corcoran et al 08 data



/* Like 06, but version for Words or Beads. Using interpolation for faster MC integrations to implement better M phase formulae, as of Apr 08. Again two 'jars', B and G, which have a higher proportion of b and g type of beads. There are actions choose B (chB), choose G (chG) and sample (S). Action values - Q(X,Ng,N) = value of taking action x when there have been N draws so far, from which Ng are of type g. Basic procedure: 1. Estimate 'ideal observer' probabilities re. Jar-as-Cause prob.s 2. Take into account possible costs and derive Q values for chosing actions 3. Fit with one set of results, and apply to another ... Version 05rbw uses data from Wellcome 07 Richard Bentall et al study, and concentrates on exploration of noise-param. vs. cost-of-slowness slice of the parametre space. */ #include #include #include #include // input & output streams #include // i-o parametrized manipulators #include // in-memory string manipulations for i/o, new standard #include #include #include "jtc.h" // Beads etc. setup classes. #include "jtc_exp.h" // Experimental procedures/setups based on jtc.h #include "jtc_em.h" // EM and such like operating on stuff in jtc_exp.h; // contains random distro, cummul distro & multiple integrn. stuff. #include "maths.h" #include "nr_ut.h" #include "menu.h" #include "blitz_array_aux.h" // ******************** Global Variables & Objects ***************************** // model Objects & Variables --------------------------------------------------- Beads_Exp_II GB ; // General Beads with iteration class RBWelcomeDat RBWelc ; // provides sequence & exp. data from Wellcome study // Various model Parameters ---------------------------------------------------- int groupN =1; // would be 3 for RB expts, 4 for Fear & Hayley expt. simulation int deprLab=0; // For RPB simulations. int cntrLab=1; int paranLab =2; // for all int ocdLab=3; int mixedLab =4; // for F&H etc. float nIter = 10; // Default no. of EM iterations float mcalls = 60000; // default calls to Vegas for integr float dataGroup = 1; // <=0 numbers - read from various files; esp. -2. // >0 numbers - load respective RBW data float expType = 1; // 1 for beads, 2 for words // Accessory Objects ----------------------------------------------------------- Menu25 menu; float **menuNo; char **menuSt; // to connect to menu params. // ************************ In-file Functions ********************************** double poidPDat_if_macropPDF( double *q, size_t dim, void *param); double baysPDat_if_macropPDF( double *q, size_t dim, void *param); // to be used to calculate means and variances of kT and CS: double Tm_baysPDat_if_macropPDF( double *q, size_t dim, void *param); double Tv_baysPDat_if_macropPDF( double *q, size_t dim, void *param); double CSm_baysPDat_if_macropPDF( double *q, size_t dim, void *param); double CSv_baysPDat_if_macropPDF( double *q, size_t dim, void *param); // versions using interpolation: double baysPDat_if_macropPDF2( double *q, size_t dim, void *param); // to be used to calculate means and variances of kT and CS: double Tm_baysPDat_if_macropPDF2( double *q, size_t dim, void *param); double QlnT_baysPDat_if_macropPDF2( double *q, size_t dim, void *param); double CSm_baysPDat_if_macropPDF2( double *q, size_t dim, void *param); double QlnCS_baysPDat_if_macropPDF2( double *q, size_t dim, void *param); void MakeMenu(); // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~ MAIN ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ int main() { int i, j, k, m, s; // various auxiliary variables long int n; double aux0, aux1, aux2, aux3; // useful even in debugging float rCSm, rCSv, rTm, rTv; // 'reference' storage values. double tempCSm, tempCSv, tempTm, tempTv, IntQlnCS, IntQlnT; //temporary values. // could have 'expType' for Beads, Words or combined, but not no ! // Parametrest to stay const (not all needed ...): GB.CGB_N= GB.kT_N = GB.CS_N= GB.rCW_N= GB.q_N= GB.PoG_N= GB.CDPar1_N= 1.0; GB.CGB =GB.CBG = -100.0 ; int ptN =33; // 69; // 36; Array sumRes(Range(0,roundfl(nIter)),Range(1,8)); sumRes =0; GB.ResetAndSubjN(ptN); // Initialize : GB.currCSm = rCSm = -1.5 ; GB.currCSv = rCSv = 17.0 ; GB.currTm = rTm = 4.75 ; GB.currTv = rTv = 11.0 ; GB.kT_N = 80; GB.CS_N = 120. ; // *********************************** Menu *********************************** // to tell menu which parametres to adjust : menuNo = new float* [15+1]; menuSt = new char* [5] ; // see in_out.h menuNo[1]=&mcalls ; menuNo[2]=&dataGroup ; menuNo[3]=&rCSm ; menuNo[4]=&rCSv; menuNo[5]=&rTm ; menuNo[6]=&rTv ; menuNo[7]=&nIter ; menuNo[8]=&GB.kT_N ; menuNo[9]=&GB.CS_N ; menuNo[10]=&expType ; MakeMenu(); // Construct & Run menu char menuFlow = menu.GetChoices(menuNo, menuSt); // might need to make sure all params, inc. derivative & non-public ones, // are updated ! // GB.RefreshParams() ; // -------------------- Menu dependent loop : ------------------------------ while (menuFlow != 'Q') { // rem menu choices run before loop and at its end // convert menu item re. calls of MC routine etc. to integers n = 1000*(roundfl(mcalls/1000.0)); k = roundfl(dataGroup); // Find / load appropriate data set(s) : if ((k >0) && (k != 3) && (k <=9)) { // Load valid data from RBW experiment : GB.loadExpToDumDat(&RBWelc,k,roundfl(expType)); // note trial type last entry ptN = GB.DumDat.rows(); GB.ResetPDumD_L_etc2(ptN); } else if (k == -2) { // load from basic combined results, "con-par.blar" ifstream inpStr("con-par.blar"); if (inpStr.bad()) { cerr<<"Arr. read error... /n"; exit(1); } GB.DumDat.resize(Range(1,33+36),Range(1,6)); GB.DumDat=0; inpStr >> GB.DumDat; // has to be 2D ! ptN = GB.DumDat.rows(); GB.ResetPDumD_L_etc2(ptN); } else if (k == -1) { // load from paranoid-without-outlier results, "par2.blar" ifstream inpStr("par2.blar"); if (inpStr.bad()) { cerr<<"Arr. read error... /n"; exit(1); } GB.DumDat.resize(Range(1,35),Range(1,6)); GB.DumDat=0; inpStr >> GB.DumDat; // has to be 2D ! ptN = GB.DumDat.rows(); GB.ResetPDumD_L_etc2(ptN); } else { cerr<<"Uprepared for dataGroup choice given ... /n"; exit(1); } // set various bits influenced by menu: sumRes.resize(Range(0,roundfl(nIter)),Range(1,8)); sumRes =0; GB.ResetPDumD_L_etc2(ptN); // set starting guess : GB.currCSm = rCSm ; GB.currCSv = rCSv ; GB.currTm = rTm ; GB.currTv = rTv ; GB.SetApproxInf(1); write2DArray_csv("data_used.csv",GB.DumDat); cout << GB.DumDat ; cout <<"\nPriors: CS_mean: "<