em_bootstraps_02.cc - 'main' file for Bootstrap CI project



/* Like em_rb_beads_07, but objective is to produce a random sample of group pseudodata, e.g. several x 36 simulated participants, and then em the result. Version em_bootstraps_2 has menued initial condition mean values. 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 float nIter = 24; // 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 = 2; // 1 for beads, 2 for words float bootIterN = 13; // This is the no. of pseudodata groups that // will be produced and 'bootstrap-analysed' by EM float startID = 1.0; // To label outpouts to distinguish different runs. // 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 / counters long int n; double aux0, aux1, aux2, aux3, aux4, aux5, aux6, aux7; // useful even in debugging float rCSm, rCSsd, rTm, rTsd; // 'reference' storage values. float iCSm, iCSsd, iTm, iTsd; // central tendencies of init. val. distros. double eCSm, eCSsd, eTm, eTsd; // point estimates for generated data double tempCSm, tempCSv, tempTm, tempTv, IntQlnCS, IntQlnT; //temporary values. double eps = 0.25; // 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 =36; Array sumRes(Range(0,roundfl((nIter+1)*bootIterN)-1),Range(1,8)); sumRes =0; GB.ResetAndSubjN(ptN); Array keyRes(Range(1,roundfl(bootIterN)),Range(1,8)); keyRes =0; Array inCond(Range(1,roundfl(bootIterN)),Range(1,5)); inCond =0; // Initialize default reference values : GB.currCSm = rCSm = -3.04999; // -0.045845 ; the paranoid EM fit values ; rCSsd = sqrt(21.6348); // 0.001036); GB.currCSv =rCSsd*rCSsd; GB.currTm = rTm = 5.21107; // 12.069600 ; rTsd = sqrt( 3.90749); // 51.478100); GB.currTv =rTsd*rTsd; GB.kT_N = 80; GB.CS_N= 120. ; // default central tendencies for initial value distributions: // iCSm= -0.044397; iCSsd= sqrt(0.00165247); iTm= 8.68103; iTsd= sqrt(44.4617) ; // From /home/mm/michael.s/sci/simdata/Bayesian_Beads/Using_Wellcome-RPB_setup/group1words_02 : iCSm = -2.60816; iCSsd = sqrt(31.0668) ; iTm = 8.50498 ; iTsd = sqrt(151.794) ; // *********************************** 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]=&rCSsd; menuNo[5]=&rTm ; menuNo[6]=&rTsd ; menuNo[7]=&nIter ; menuNo[8]=&GB.kT_N ; menuNo[9]=&GB.CS_N ; menuNo[10]=&startID ; menuNo[11]=&bootIterN ; menuNo[12]=&iCSm ; menuNo[13]=&iCSsd ; menuNo[14]=&iTm ; menuNo[15]=&iTsd ; MakeMenu(); // Construct & Run menu char menuFlow = menu.GetChoices(menuNo, menuSt); // might need to ensure all params, inc. derivative & non-public ones, are updated ! // GB.RefreshParams() ; // elementary random no. generator seed for this level srand( ( time( NULL ) % 604800) ) ; // -------------------- 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 = lround(mcalls); k = roundfl(dataGroup); // Set pseudodata pt N etc : if ((k >0) && (k != 3) && (k <=9)) { // Read up number of participants: ptN = RBWelc.ValidN(k); // later ... GB.ResetPDumD_L_etc2(ptN); } else if ((k >= 20) && (k<=999)) // consider this to be No. of 'participants' ptN = k; else { cerr<<"Uprepared for pt Num choice given, k="<=20,<=999) : " ,'N',menuNo, menuSt, 0); menu.SetOption( menOff+3,"ref CSm : " ,'N',menuNo, menuSt, 0); menu.SetOption( menOff+4,"ref CSsd : " ,'N',menuNo, menuSt, 0); menu.SetOption( menOff+5,"ref Tm : " ,'N',menuNo, menuSt, 0); menu.SetOption( menOff+6,"ref Tsd : " ,'N',menuNo, menuSt, 0); menu.SetOption( menOff+7,"no. of EM iterations per loop : " ,'N',menuNo, menuSt, 0); menu.SetOption( menOff+8,"no. T grid points : " ,'N',menuNo, menuSt, 0); menu.SetOption( menOff+9,"no. CS grid points : " ,'N',menuNo, menuSt, 0); menu.SetOption( menOff+10,"run label starting, startID : " ,'N',menuNo, menuSt, 0); menu.SetOption( menOff+11,"no. of EM-estimation points : " ,'N',menuNo, menuSt, 0); menu.SetOption( menOff+12,"initial mean CSm : " ,'N',menuNo, menuSt, 0); menu.SetOption( menOff+13,"initial mean CSsd : " ,'N',menuNo, menuSt, 0); menu.SetOption( menOff+14,"initial mean Tm : " ,'N',menuNo, menuSt, 0); menu.SetOption( menOff+15,"initial mean Tsd : " ,'N',menuNo, menuSt, 0); menu.Refresh(); } // end of MakeMenu // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - // 2a. fn to integrate w. MC for denom. of Bayesian model recog. distr. // using Interpolated results. double baysPDat_if_macropPDF2( double *q, size_t dim, void *param) { // use of *param not needed as fn. has access to global objects double aux1, aux2, aux3; // double keep_kT; keep_kT=GB.kT; int d1, d2, d3 ; // d1 = lround(GB.DumDat(roundfl(GB.curr_pt),3)); d2 = lround(GB.DumDat(roundfl(GB.curr_pt),4)); d3 = lround(GB.DumDat(roundfl(GB.curr_pt),5)); GB.kT = q[0]; GB.CS = q[1]; // NEEDED by wrkP2 ... // cout <<" kT:"<