jtc_exp.h: Header file with experiment structure for "Beads-in-a-Jar" simulations



// File with experimental setups for jumping to conclusions / // beads-in-a-jar and related systems. #ifndef JTC_EXP_H #define JTC_EXP_H #include #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 using namespace std; #include "maths.h" #include "nr_ut.h" #include "menu.h" #include "jtc.h" #include "blitz_array_aux.h" #include #include BZ_USING_NAMESPACE(blitz) /* A.a. Summary data from Beads and related experiments : a. */ /* = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = */ /* B.a. Summary data from Fear & Healy 1997 */ class FearHealy97dat { public: FearHealy97dat(); ~FearHealy97dat(); Array Cond1Tab2 ; /* This is: // Control OCD DD Mixed 0.54, 0.53, 0.51, 0.52, // Init. Certainty 0.10, 0.18, 0.06, 0.17, // " " SD 2.6, 3.4, 1.5, 2.7, // Draws to decision 1.3, 2.5, 0.9, 1.8, // " " SD 0.967, 0.923, 0.773, 0.933, // participant fraction deci- // ding on 85%-black jar, A. 30, 26, 22, 15 ); // and N per group. */ float PyBlackBallGivenJarA85(){ return 0.85; } ; }; /* = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = */ /* B.b. from Garety_et_al_91 */ class Garety91dat { protected: static const char *const ERROR; int gLab, bLab; // conventional labels, gLab=1, bLab=-1 void CheckInd(int ind){ if (ind<1 || ind>20) Error("%s%s \n",Garety91dat::ERROR,"Index off range" ); }; public: int SamplMaxN(){ return 20; }; // allow free access - we trust our users ! This is the standard bead // sequence used in Young & Bentall 97 and various others, and its inverse. float *gseq ; // e.g. green-starting or good-trait-starting in Young & Bentall: 11101... float *bseq ; // e.g. yellow-starting or bad-trait-starging " " " : 00010... // These are the number-of-green/good as fn. of no. sampled : float *ng_gseq; float *ng_bseq; Garety91dat(); ~Garety91dat(); // the following return integers: // Sequence draws : int gSeq(int ind){ CheckInd(ind); return roundfl(gseq[ind]); }; int bSeq(int ind){ CheckInd(ind); return roundfl(bseq[ind]); }; // ng (number of green/good results) out of nd draws of ... int ngGSeq(int nd){ if (nd<0 || nd>20) // ... the gseq Error("%s%s \n",Garety91dat::ERROR,"ngGSeq index off range"); return roundfl(ng_gseq[nd]); } ; int ngBSeq(int nd){ return nd - ngGSeq(nd) ; } ; // ... the bseq }; /* = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = */ /* B.b. from Bentall_et_al_06 ('Welcome study'- - The Phenomenology and Cognitive Structure...) and also the Mackintosh - El-Deredy variant */ class RBWelcomeDat { protected: static const char *const ERROR; static const int sN = 3; // auxiliary for No. of sequences presented int gLab, bLab; // conventional labels, gLab=1, bLab=-1 // Here I have called Yellow=g=1, black=b=0. int *validN; // holds the valid N for each group in Wellcome study by group // label, e.g. validN[0]=0, validN[1]=36 etc. int meRetN; // Number of different reward/cost pattern in M-E type expts. void CheckInd(int ind){ if (ind<1 || ind>20) Error("%s%s \n",RBWelcomeDat::ERROR,"Index off range" ); }; Array centiles; // working space for e.g. quartile values of // expt. data. Is in form (without the header numerical labels) : // \ group: 1 2 4 ... 9 1 2 4 ... 9 // centile \ exp beads1 beads2 ... // ------------------------------------------------------------------------ // 0 (i.e. min val.) g1min g2min ... // 25 // ... // 100 (i.e. max val) public: int SamplMaxN(){ return 20; }; int ptN(){ return 229; } ; // Tot. num. of participants with valid data for // paranoid processing. // allow free access - we trust our users ! This is the standard bead / word // sequences used in Young & Bentall 97 and various others, and its inverse. float **seq; // will point to the ones below. [0] wasted. float *seq1 ; // BYBBBBYBBBYBYYYYBYYY , here 01000010001011110111 float *seq2 ; // BYBBBYBBYBYBBBBYYBBY , here 01000100101000011001 float *seq3 ; // YBYYYYBYYYBYBBBBYBBB , here 10111101110100001000 // These are the number-of-yellow/good as fn. of no. sampled for Welcome: float **ng_tot; // will point to the running totals below float *ng_seq1; float *ng_seq2; float *ng_seq3; // Makintosh - El-Deredy param : float CS1,R1,CS2,R2 ; // usu. -5.0, 100.0, -10, 200 as for lo-cost cond and hi-cost // for Mackintosh - El-Deredy version (M.. E-D.. Seqences -> initials): float **MEseq; // will point to the ones Makintosh - El-Deredy used. [0] wasted. float *mes1 ; // BRBBBRBBRBRBBBBRRBBR, here 01000100101000011001 (like seq2) float *mes2 ; // RBRRRRBRRRBRBBBBRBBB, here 10111101110100001000 (like seq3) float *mes3 ; // BBRBBRBRRBRBBBBRRBBR, here 00100101101000011001 (new) // These are the number-of-yellow/good as fn. of no. sampled for Makintosh - El-Deredy: float **ng_met; // will point to the running totals below float *ng_mes1; float *ng_mes2; float *ng_mes3; /* See seqByCostGpME1 below ! So not yet: sequences used for each cost pt group: int **costgp_seq; // starts from costgp_seq(1,1) { gp1->112,2->123,2->231 */ // Big array with all the VALID JTC data. Some secondary variables have value // 9999 when invalid. This will be a 229 cases x 15 data-columns array. Array WJTC ; // Makintosh - El-Deredy data; ME1 is 1st costed exp. : Array ME1JTC ; RBWelcomeDat(); ~RBWelcomeDat(); // --------------- Methods of general applicability in this class --- // General Lookups : int seqN(){ return sN; }; // Either for Wellcome or ME variants float PyG_Given_Jar_G_90(){ return 0.90; } ; float PyG_Given_Jar_G_75(){ return 0.75; } ; float PyG_Given_Jar_G_60(){ return 0.60; } ; // .------------------------------------------------------------------. // | Methods concerning Welcome I - RPB data | // |__________________________________________________________________| int ValidN(int k){ return validN[k]; } ; // the following return integers: int Seq1(int ind){ CheckInd(ind); return roundfl(seq1[ind]); }; int Seq2(int ind){ CheckInd(ind); return roundfl(seq2[ind]); }; int Seq3(int ind){ CheckInd(ind); return roundfl(seq3[ind]); }; // ng (number of yellow/good results) out of nd draws of ... int ngSeq1(int nd){ if (nd<0 || nd>20) // ... the seq1 Error("%s%s \n",RBWelcomeDat::ERROR,"ngSeq1 index off range"); return roundfl(ng_seq1[nd]); } ; int ngSeq2(int nd){ if (nd<0 || nd>20) // ... the seq2 Error("%s%s \n",RBWelcomeDat::ERROR,"ngSeq2 index off range"); return roundfl(ng_seq2[nd]); } ; int ngSeq3(int nd){ if (nd<0 || nd>20) // ... the seq3 Error("%s%s \n",RBWelcomeDat::ERROR,"ngSeq3 index off range"); return roundfl(ng_seq3[nd]); } ; // reset bead drawN in sequence seqN to newVal, and make other members of the class // consistent. Of course if change off initial sequences, can't analyse decisions-to // draw data in the same way but can use for simulated experiments etc. void beadSet(int seqN, int drawN, int newVal); // Centile work: // Lookup, with all sort of checks. expType is beads=1, words=2; // expNum is 1 to 3 (the seq. num. used). Returns integer d2d value : int CentVal( int centile, int groupLabel, int expType, int expNum); // fill in the centiles array for all the groups. // partNum=4 would mean quartiles, etc.: void centileCalc(int partNum); // .------------------------------------------------------------------. // | Methods concerning Mackintosh - El-Deredy type data | // |__________________________________________________________________| int returnPatN(){return meRetN; }; // the following gives the sequence used for each cost level for each of the // three expt. groups of expt. ME1. Note it can be used with variable subject // number so that full or excl-20 datasets, or simulated datasets, can be used. void seqByCostGpME1(int this_pt, int *S1, int *S2, int *S3, int PtN) { if (this_pt<= PtN*2/3 ) { *S1=*S2=1; *S3=2; } else if (this_pt<= PtN*5/6 ) { *S1=*S3=3; *S2=2; } else if (this_pt<= PtN ) { *S1=2; *S2=3; *S3=1; } else { cerr<<"/n pt. number mismatch in RBWelcomeDat::seqByCostGp. /n"; exit(1); } }; }; /* = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = */ /* C.a. Classes with the 'Beads' machinery that can do more */ /* First, can also do just iteration: */ class Beads_Exp_I : public Beads // See attenmod.h -> class Exp_Dual_At_K_II : public DualAtK_II ... { protected: static const char *const ERROR; // The following are used by RefreshIter: float **iterParPtr, **iterIterNo, **iterEndPtr ; Iter iter ; // auxiliary stuff for file output : static const int fN = 21; // auxiliary for No. of P components / output files. ofstream **ofs ; // for files that carry Py(draws-to-dec.=n) */ ofstream Pd0, Pd1, Pd2, Pd3, Pd4, Pd5, Pd6, Pd7, Pd8, Pd9, Pd10; ofstream Pd11,Pd12,Pd13,Pd14,Pd15,Pd16,Pd17,Pd18,Pd19,Pd20 ; ofstream l3beads, l3words ; // for likelihood products for one participant // using the Bayesian ideal obsrver reasoning etc. ofstream cl3beads, cl3words ; // As immediately above for control distribution. ofstream results; // For summary and comparative results over all participants ! // will also be used for debugging until definitive prog. version ready. // 'results' will be aka ofs[fN+4] . int seqNum, seqLen; // local copies of number and length of draw-sequences // that the experiment uses - def. 20 (-> P00 to P20 !), 3 int retNum; // [local copy of] the number of different patterns of returns that // the expt. uses - e.g. 3 for ME1 that has zero, low and high. public: float llrSD0; // Basically an auxiliary to give the prior distibution of // log-lik-ratios (before any evidence) a small nonzero width. float CGB_e, CGB_N; // endpoint and number of CGB float kT_e, kT_N; float CS_e, CS_N; float rCW, rCW_e, rCW_N; // This is the ratio CGB/CBG, putatively similar // to false-neg.cost/false-pos.cost float q_e, q_N; // in case we need to iterate the P(g|G) float PoG_e, PoG_N; // in case we need to iterate the PoG. float CDPar1, CDPar1_e, CDPar1_N; // to handle & iterate contr. distr. param. 1 float CDPar2, CDPar2_e, CDPar2_N; // ... and 2. Array PVecs; // 4-D array: param1 x param2 x sequence_Num x Probability_vector // e.g: CS_N x CGB_N x 3 (for Wellcome) x fN // If trials have different Array Res; // for final output results - 15 columns or so (see ctr etc.): // BayesianLik3Beads BayesianLik3Words ContrLik3Beads ContrLik3Words BestParBayesBeads1 2 3 BestBayesianParamWords1 2 3 ControlParBeads1 2 ContolParWords1 2 // OptL3B OptPar1B OptPar2B OptPar3B \ // OptL3W OptPar1W OptPar2W OptPar3W \ // ContrL3B ContrPar1B \ // ContrL3W ContrPar1W // for file output, plan is to add the pt. numbers, // [the control distribution best-likelihoods and corresponding params] // and the two derivative ratios, Opt/Contr_L3B and Opt/Contr_L3W, // as additional cols. Interpolation2D IP2; Beads_Exp_I(); ~Beads_Exp_I(); // paramIter access fns. Initial 'Set'ting in ctor void RefreshIter() ; // for after each menu-dep input int DoIter() { return iter.Do() ; } ; int StepIter() { return iter.Step() ; } ; float StartIter(int j) { return iter.Start(j) ; } ; float EndIter(int j) { return iter.End(j) ;}; // determine total no. of iterations to do; also assigns result to totRunN unsigned long TotIterN() { return iter.TotN() ; } ; unsigned long ItersDone() { return iter.ItersDone() ; } ; int AtStepIter(int k){ return iter.AtStep(k); }; float FractRunsDone(){ return iter.FractDone() ; } ; /* void ShowProgress(); // pars, fract done etc -> cerr */ // resetting etc. // PreparePVecs1() etc. use existing params. Inline for clarity rather than speed. // PreparePVecs1() is for the Bayesian Costed model (kT/CS), PreparePVecs2() for the // SPRT (kT/CDPar1), PreparePVecsME1() for 1st Mackintosh-El-Deredy Bayesian model etc: void PreparePVecs1(){ PVecs.resize(Range(1,roundfl(*iterIterNo[1])),Range(1,roundfl(*iterIterNo[2])),Range(1,seqNum),Range(0,seqLen)); PVecs = 0; } ; void PreparePVecs2(){ PVecs.resize(Range(1,roundfl(*iterIterNo[1])),Range(1,roundfl(*iterIterNo[6])),Range(1,seqNum),Range(0,seqLen)); PVecs = 0; } ; void PreparePVecsME1(){ //rem iterParPtr[1]= &kT; iterParPtr[3]= &CGB; int d3 = seqNum*retNum; PVecs.resize(Range(1,roundfl(*iterIterNo[1])),Range(1,roundfl(*iterIterNo[3])),Range(1,d3),Range(0,seqLen)); PVecs = 0; } ; // Code for next one is in jtc_exp.cc, similar to above. void PreparePVecsME2(); // rem iterParPtr[1]= &kT; iterParPtr[2]= &CS; // ** Methods for writing results (mostly 2D arrays) to disk ** // Open output files, jst before itern. over params.: int OpenOutp() ; // Next write to the open s'th Pd2d file stream, from within param. i // itern., loop, assuming first iter. par. is CS (or, acc. to Peter, beta) // and second iter. par is CW (alpha acc to Peter). // Non-trivial !! int WritePd2d(int s, double result, char sepChar = '\t', int check=1 ); // Next one calculates the whole P(d2d=0, =1, ...) vector and writes to // all the appropriate streams. Assumes : int WriteAllPd2d( float *seqAr, int seqArEnd=20, int seqArStart=1 , char sepChar = '\t', int check=1 ); // Method for writing just two likelihood products for the current values of // the parameters, in same style as StoreAllPVecs, for just 1 participant : int Write2Prod3P( int part_i, RBWelcomeDat *RBWptr, int check=1, char sepChar = '\t') ; // Another version, similar to above, but designed so that first iter par // is kT and second is CS : int WriteAllPd2d_2( float *seqAr, int seqArEnd=20, int seqArStart=1 , char sepChar = '\t', int check=1 ); // Method for writing just two likelihood products for the current values of // the parameters, in same style as StoreAllPVecs, for just 1 participant: int Write2Prod3P_2( int part_i, RBWelcomeDat *RBWptr, int check=1, char sepChar = '\t') ; int CloseOutp(); // Close open file streams at end of param. itern. // Method for writing 21-D Bayesian-der. Prob. vectors to 4-D array of the form // PVecs.resize(Range(1,roundfl(GB.CS_N)), Range(1,roundfl(GB.CGB_N)), Range(1,RBWelc.seqN()), Range(0,20)) // validZero var says if decision at 0 draws is allowed. // Version using RBWelcomeDat : int StoreAllPVecs( RBWelcomeDat *RBWptr, int validZero=1 ) ; // Mackintosh-El-Deredy like protocols can't just store PVecs for a given // set of microparam; they have to recalc. SoftME$_M_PG_V_QG_QB_QS_AS_MS(zeroAllowed) // for the different cost/reward patterns too. All included in the following, // which sets PVecs(AtStepIter(1),AtStepIter(3), seq_n+retTot*(ret_n-1), s)=... int CalcMEPVecs( RBWelcomeDat *RBWptr, int validZero=1, int version = 1) ; // Return the product of the probabilities for the 3 trials of a certain // type, type 1 for beads, 2 for words, for the current combination of // parameters (CS, CGB, CBG, kT ...), for participant w. index part_i: double Prod3P( int part_i, RBWelcomeDat *RBWptr, int drawType=1, int check=1); // Control distributions & Py products : // 1. Poisson, except stage 20 which mops up the remainder of all the rest. // If param is <0, it uses CDPar1 as the lambda (or lambda*t) of the Poisson // distr. double CPDistr1(int d2d, double param = -666 ) ; // Likelihood products for one pt. using Poisson-like above, using // parameter for CDPar1 : double cProd3P(RBWelcomeDat *RBWptr, int participant, double param, int drawType=1 ); // find max-likelihood CDPar1 param within range CDPar1:CDPar1_e, CDPar1_N : double bestP3PPoissonoidParam( RBWelcomeDat *RBWptr, int participant , int drawType=1, double *Prod3Val = NULL ); // 2. Sequential Probability Ratio Test using CDPar1 // a. Directly calculate the error-free log-likelihood-ratio for a particular // step of particular sequence. double SPRT(int step, int seq, RBWelcomeDat *RBWptr); // zero max likelihoods and other results : void zero_B_W_3L(){ Res = 0; } ; // track max likelihoods for current pt - if current best than stored, update !: void update_B_W_3L( double L3B, double L3W, int ptInd ); // calculate and store in memory (in Res col. 9 to 12) results of control distr 1: void contr1Res(RBWelcomeDat *RBWptr) ; // Writing main results for pts. If participant=0 write for all pt, if ti=0 // ommit the title line. void WritePtRes(RBWelcomeDat *RBWptr, int participant=0, int ti=1 ); }; // Don't forget the semilcolon at the end of class declarations ! // = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = #endif // JTC_EXP_H // eof