Please cut the rest of this message into six files with the names: readme.txt, notice.txt, bootprob.c, test1.dat, test1b.dat, test2.dat. ========================readme.txt==CUT HERE============================= This readme file contains information about two programs, bootprob.c and bootprobdiff.c. -- bootprob.c DHF, WFB Ver 2.5 16-Oct-06 (output fitted curve for plotting; suppress delay; increase maximum levels and repetitions) WFB Ver 2.4 24-Nov-96 (replace all routines from Numerical Recipes) DHF Ver 2.3 18-Apr-96 (minor change to conform to ANSI C) Ver 2.2 11-Mar-95 (cosmetic change to chi-sqr output) 20-Apr-95 (update of timing data under Timing) Ver 2.1 20-Jan-92 This C program fits a cumulative Gaussian psychometric function to a set of binomial data and saves the fitted curve in a separate file called fitted_curve.txt. It gives the threshold (for a given criterion level of performance), the gradient (slope), and 1/gradient (spread) for the fitted curve. It then estimates the standard deviations (SDs) and confidence intervals (CIs) of the threshold, slope, and spread by a bootstrap procedure (the 'parametric' bootstrap). The gradient is of the underlying linear regression, and, for a criterion level of performance of p = 0.5, the value needs to be scaled by 1/sqrt(2*pi) to obtain the gradient of the cumulative Gaussian at this threshold; likewise with the SD and CIs of the gradient. The bootstrap procedure is similar to that given in Foster and Bischof (1991; see below), but a more robust procedure has been adopted here in that SDs are computed from centiles, assuming an approximately normal distribution. When that assumption is violated, CIs should be used anyway. The CIs are based on the basic percentile method rather than on the bias-corrected and accelerated (BCa) method; see Efron, B. (1987) Journal of the American Statistical Association, 82 (397), 171-185. Differences between the two kinds of estimates were found to be small. Since this is a resampling program, repeated applications will produce slightly different values of the estimated SDs and CIs, the size of the difference depending on the number of bootstrap replications. To run the program, type: bootprob 1000 where is the full name (including any extension) of the datafile and 1000 is the number of bootstrap replications. The file is formatted thus: 5 0.5 -2.0 5 0 -1.0 5 1 0.0 5 3 1.0 5 4 2.0 5 5 where the first row has the number of levels (5) and the criterion level of performance (0.5), and the second and subsequent rows give the stimulus level (e.g. -2.0), the number of trials at that level (e.g. 5) and the number of hits or successes (e.g. 0) at that level. Any text may be written after the last line of data. The maximum number of lines of data is 500, but this may be increased by changing the maximum value in the program. The number at the end of the command line (1000) is the number of bootstrap iterations, which can be reduced to e.g. 200 if just SDs are to be used. The maximum number of bootstrap replications is 5000, but this may be increased by changing the maximum value in the program. As an example on a Windows platform, open a command-prompt window in the directory containing both the executable program bootprob.exe and the data file test1.dat (or any other data files) and type: bootprob test1.dat 1000 This program has been tested against corresponding large-scale bootstrap analyses using the probit link function under GLIM (Numerical Algorithms Group Oxford). For the same number of bootstrap iterations (1000), differences in SDs were less than 1% with the data set test2.dat and less than 7% with the data set test1.dat. The latter set has only 25 trials in all, and the increase in error is due to differences in convergence when pathological samples are generated. In general, for data sets with few trials, SDs and CIs can be strongly influenced by these pathological samples. Although all reasonable effort has been invested in validating this program, it does not have the precision and traps for outliers described in Foster and Bischof (1991), nor has it been validated on as full a set of test data. Please read the notice.txt file. Language: The source code is in ANSI C. Timing: For a full 1000-iteration bootstrap, the program should take a few seconds on a 486 or less on a 586 PC. Source reference: Foster, D.H. and Bischof, W.F. (1991) Thresholds from psychometric functions: superiority of bootstrap to incremental and probit variance estimators. Psychological Bulletin, 109, 152-159. Introduction to software: Foster, D.H. and Bischof, W.F. (1997) Bootstrap estimates of the statistical accuracy of thresholds obtained from psychometric functions. Spatial Vision, 11, 135-139. Algorithms: Routines were developed by W.F. Bischof and D.H. Foster based on the following sources: D.J. Finney (1947, 1952, 1971) Probit Analysis, C.U.P., Cambridge, England. Odeh, R.E. and Evans, J.O. (1974) Algorithm AS 70: Percentage points of the normal distribution, Appl. Stat. 23, 96-97, given in Kennedy, W.J. and Gentle, J.E. (1980) Statistical Computing, Dekker, New York. Cody, W. J. (1969). Rational Chebyshev approximations for the error function, Math. Comp., 23, 631-638. Described in Kennedy, W. J. & Gentle, J. E. (1980) Statistical Computing, New York: Dekker, pp. 90ff. Kachitvichyanukul, V. and Schmeiser, B. W. Binomial Random Variate Generation. Communications of the ACM, 31, 2 (February, 1988) 216. Random number routines from BSD sources, copyright 1993 Martin Birgmeier. Further details on sources are given in bootprob.c Authors: David H. Foster and Walter F. Bischof -- bootprobdiff.c DHF, WFB Ver 1.0 15-Oct-06 (difference of two thresholds, based on bootprob, Ver 2.1) This C program fits cumulative Gaussian psychometric functions to two sets of binomial data and saves the fitted curves in separate files. It gives the threshold (for a given criterion level of performance), the slope, and spread for each of the fitted curves, and the differences in threshold, slope, and spread. It then estimates the standard deviations (SDs) and confidence intervals (CIs) of these differences by the bootstrap procedure. The CIs are based on the basic percentile method rather than on the bias-corrected and accelerated (BCa) method; see Efron, B. (1987) Journal of the American Statistical Association, 82 (397), 171-185. Differences between the two kinds of estimates were found to be small. The CIs can be used to decide whether the two thresholds (or two slopes or two spreads) are significantly different from each other. For example, if the 95% CIs for the threshold difference does not contain zero, then the thresholds may be said to be significantly different, p < 0.05, two-tailed test. To run the program, type: bootprobdiff 1000 where is the full name (including any extension) of the first datafile; is the full name (including any extension) of the second datafile; and 1000 is the number of bootstrap replications. For details of the format of the datafile, see the readme file. As an example on a Windows platform, open a command-prompt window in the directory containing both the executable program bootprobdiff.exe and the two data files test1.dat and test1b.dat (or any other pair of data files) and type: bootprobdiff test1.dat test1b.dat 1000 This program has been tested against an implementation in Splus (Ver. 7.0 for Windows, Insightful Corp.) by K. Zychaluk. For the same data sets and number of bootstrap iterations (1000), but with extreme values excluded, the estimates for a 90% CI for the difference in thresholds were as follows: Splus percentile: [9.2024, 10.7678] Splus BCa: [9.1785, 10.7246] bootprobdiff (percentile): [9.1284, 10.8129] The signs of the limits depend on whether test1.dat is used as datafile1 and test1b.dat as datafile2 or vice versa. Authors: David H. Foster and Walter F. Bischof -- Compilation hints: Here are a few examples for compiling bootprob and bootprobdiff. Windows (Borland C/C++) bcc32 -O2 bootprob.c bcc32 -O2 bootprobdiff.c Linux (gcc) cc -O2 -o bootprob bootprob.c -lm cc -O2 -o bootprobdiff bootprobdiff.c -lm Mac OSX (gcc) cc -O2 -Dunix -o bootprob bootprob.c -lm cc -O2 -Dunix -o bootprobdiff bootprobdiff.c -lm ========================notice.txt==CUT HERE============================= Notice: This program and procedures are provided free of charge. No warranty, express or implied, is made that this program and procedures are free of error, or are consistent with any particular standard of quality, or that they will meet your requirements for any particular application. They should not be relied on for solving a problem whose incorrect solution could result in injury to a person or loss of property. If you do use the program or procedures in such a manner, it is at your own risk. ========================bootprob.c==CUT HERE============================= /* bootprob.c DHF, WFB Ver 2.5 16-Oct-06 (output fitted curve for plotting; suppress delay; increase maximum levels and repetitions) WFB Ver 2.4 24-Nov-96 (replace all routines from Numerical Recipes) DHF Ver 2.3 18-Apr-96 (minor change to conform to ANSI C) Ver 2.2 11-Mar-95 (cosmetic change to chi-sqr output) 20-Apr-95 (update of timing data under Timing) Ver 2.1 20-Jan-92 This C program fits a cumulative Gaussian psychometric function to a set of binomial data and saves the fitted curve in a separate file called fitted_curve.txt. It gives the threshold (for a given criterion level of performance), the gradient (slope), and 1/gradient (spread) for the fitted curve. It then estimates the standard deviations (SDs) and confidence intervals (CIs) of the threshold, slope, and spread by a bootstrap procedure (the 'parametric' bootstrap). The gradient is of the underlying linear regression, and, for a criterion level of performance of p = 0.5, the value needs to be scaled by 1/sqrt(2*pi) to obtain the gradient of the cumulative Gaussian at this threshold. The bootstrap procedure is similar to that given in Foster and Bischof (1991; see later), but a more robust procedure has been adopted here in that SDs are estimated from centiles, assuming an approximately normal distribution. When that assumption is violated, CIs should be used anyway. The CIs are based on the basic percentile method rather than on the bias- corrected and accelerated (BCa) method; see Efron, B. (1987) Journal of the American Statistical Association, 82, 171-185. Differences between estimates by the two methods were found to be small; see the example calculation of the CIs of threshold differences based on bootprobdiff in the readme file. Since this is a resampling program, repeated applications will produce slightly different values of the estimated SDs and CIs, the size of the difference depending on the number of bootstrap replications. To run the program, type: bootprob 1000 where is the full name (including any extension) of the datafile and 1000 is the number of bootstrap replications (which can be reduced if just SDs are to be used). For details of the format of the datafile, see the readme file. As an example on a Windows platform, open a command-prompt window in the directory containing both the executable program bootprob.exe and the data file test1.dat (or any other data files) and type: bootprob test1.dat 1000 Source reference: Foster, D.H. and Bischof, W.F. (1991) Thresholds from psychometric functions: superiority of bootstrap to incremental and probit variance estimators. Psychological Bulletin, 109, 152-159. Introduction to software: Foster, D.H. and Bischof, W.F. (1997) Bootstrap estimates of the statistical accuracy of thresholds obtained from psychometric functions. Spatial Vision, 11, 135-139. Algorithms: Routines were drawn from Finney, D.J. (1947, 1952, 1971), Probit Analysis, C.U.P., Cambridge, England; and other sources as detailed later. Language: The source code is in ANSI C. Further details: For further information, see the readme.txt file. Notice: This program and procedures are provided free of charge. No warranty, express or implied, is made that this program and procedures are free of error, or are consistent with any particular standard of quality, or that they will meet your requirements for any particular application. They should not be relied on for solving a problem whose incorrect solution could result in injury to a person or loss of property. If you do use the program or procedures in such a manner, it is at your own risk. David H. Foster and Walter F. Bischof */ #include #ifndef unix #include #include #endif #include #include #include #ifndef unix #include #else #include #endif #define TRUE 1 #define FALSE 0 #define MAXPROB 5.0 /* maximum value for probit transformation */ #ifndef MAXFLOAT #define MAXFLOAT 3.37E+38 /* maximum absolute value of midpoint */ #endif #define EPS_C 0.0001 /* stopping criterion for convergence */ #define EPS_Q 1.19E-07 /* minimum divisor; so invcgauss(0.0) = -5.17 */ #ifndef PI #define PI 3.14159265358979323846 #endif #define SQRT2 1.41421356237309504880 #define NLEVEL 500 /* maximum number of stimulus levels */ #define NBOOT 5000 /* maximum number of bootstrap replications */ #define MBIG 1000000000L #define MSEED 161803398L #define MZ 0 #define FAC (1.0/MBIG) #define ABS(x) ((x) >= 0 ? (x) : -(x)) #define MIN(a,b) ((a) <= (b) ? (a) : (b)) float meanran=0.0, count=0.0; /* global, used for ran3 check */ int ifails=0; /* global, used for fit check */ float ord(float z); float fsqr(float x); FILE *strmdata; /* data file */ char datafile[100]; int probit(int nlevel, float slevel[], int trials[], int hits[], float pcrit, float *pa, float *pb, float *pchi, int *pdof, float param[]); void boot(int nlevel, float slevel[], int trials[], int hits[], float prcrit, float a, float b, int nboot, float param[], float thrlimits[], float spreadlimits[], float slopelimits[]); int regress(int n, float x[], float y[], float w[], float pcrit, float *pa, float *pb, float *pchi, int *pdof, float param[]); void sort(int n, float ra[]); float centiles(int nboot, float ra[], float limits[]); int degen(int nlevel, int trials[], int hits[]); float cgauss(float x); float invcgauss(float p); float erfcc(float x); int binom_rand(float p,int n); void srand48(long seed); double drand48(void); /* ----------------------------------------------------------------------- */ int main(int argc, char *argv[]) { float slevel[NLEVEL]; float thr, slope, spread; float thr0, slope0, spread0; float sdthr, sdslope, sdspread; float a, b, chi, param[6], thrlimits[4], spreadlimits[4], slopelimits[4]; float pcrit; float xbegin, xend, xinc, xfine, pfine; /* for pretty fitted curve */ int nlevel, hits[NLEVEL], trials[NLEVEL], dof; int i; int nboot; int overflow; char dump[100]; FILE *strmfitcrv; /* fitted curve */ #ifndef unix // void clrscr(void); system("cls"); #endif printf("bootprob.c Ver 2.4\n\n" "Please read the readme.txt and notice.txt files on disk.\n\n" "Reference:\n" "Foster, D.H. and Bischof, W.F. (1991) Thresholds from\n" "psychometric functions: superiority of bootstrap to\n" "incremental and probit variance estimators. Psychological\n" "Bulletin, 109, 152-159.\n"); srand48((long)19961996); if (argc <= 1) { /* test for command line args */ printf("Data file not given!\n"); exit(1); } if (argc > 1) { /* read data file */ strcpy(datafile,argv[1]); strmdata = fopen(datafile, "r"); if (strmdata == NULL) { printf("\nFailed to open data file \"%s\"!\n", datafile); exit(1); } if (strmdata == NULL) { printf("Failed to open data file\n"); exit(1); } fgets(dump, 100, strmdata); /* get first line of data set up to CR */ if (sscanf(dump, "%d %f", &nlevel, &pcrit) < 2) { printf("\nNumber of levels or criterion performance not given!\n"); exit(1); } printf("\nInput data:\n"); printf(" number of levels: %d, criterion level of " "performance: %4.2f\n", nlevel, pcrit); if (nlevel > NLEVEL) { printf("\nNumber of stimulus levels > %d. Increase source code limit NLEVEL.\n", NLEVEL); exit(1); } if (pcrit <= 0.0 || pcrit >= 1.0) { pcrit = 0.5; /* default */ printf("assuming that criterion p = 0.5\n"); } for(i = 0; i < nlevel; i++) { fscanf(strmdata, "%f %d %d", &slevel[i], &trials[i], &hits[i]); printf(" %8.4f %3d %3d\n",slevel[i],trials[i],hits[i]); } (void)fclose(strmdata); /* close input datafile */ } if (argc == 3) nboot = atoi(argv[2]); /* read nboot from command line */ else { nboot = 1000; /* default */ printf("Default number of bootstrap replications set to 1000.\n"); } if (nboot < 1000) { printf("\nNumber of bootstrap replications is less than 1000!\n" "Result will be approximate. If possible, use 1000 instead.\n"); } if (nboot > NBOOT) { printf("\nNumber of bootstrap replications > %d. Increase source code limit NBOOT.\n", NBOOT); exit(1); } if (degen(nlevel, trials, hits) > 0) printf("\nCAUTION: possibly degenerate data set!\n"); overflow = probit(nlevel,slevel,trials,hits,pcrit,&a,&b,&chi,&dof,param); if (overflow) printf("\nCAUTION: probit fit inadequate!\n"); thr0 = param[0]; spread0 = param[1]; slope0 = 1.0/spread0; sdthr = param[2]; sdslope = param[4]; /* printf("Probit results\n"); */ /* printf("a = %-.4f b = %-.4f", a, b); */ printf("\nRESULTS OF INITIAL MAXIMUM LIKELIHOOD FIT\n"); printf("For a criterion performance level of %2d%%:\n", (int)(pcrit*100.0)); printf(" Estimate\n"); printf("threshold: %9.4f\n", thr0); printf("gradient: %9.4f\n", slope0); printf("1/gradient: %9.4f\n", spread0); /* Create and save pretty fitted curve for plotting */ xbegin = (invcgauss(0.01)-a)/b; /* begin at p = 0.01 */ xend = (invcgauss(0.99)-a)/b; /* end at p = 0.99 */ if (xbegin > slevel[0]) xbegin = slevel[0]; /* extend if necessary */ if (xend < slevel[nlevel-1]) xend = slevel[nlevel-1]; xinc = (xend - xbegin)/50; /* 50 points per curve */ strmfitcrv = fopen("fitted_curve.txt", "w"); for (xfine = xbegin; xfine <= xend; xfine += xinc) { pfine = cgauss(a+b*xfine); fprintf(strmfitcrv, "%f %f\n", xfine, pfine); } fclose(strmfitcrv); printf("\nNote 1: The file \"fitted_curve.txt\" contains plot data for the" "\n psychometric function.\n"); printf("\nNote 2: The gradient refers to the slope of the linear regression" "\n underlying the psychometric function.\n"); printf( " The lack of fit of the psychometric function is given \n" " approximately by chi-squared = %-.1f with df = %-d.\n\n", chi, dof); printf("Bootstrap starting!\n"); boot(nlevel, slevel, trials, hits, pcrit, a, b, nboot, param, thrlimits, spreadlimits, slopelimits); thr = param[0]; spread = param[1]; slope = param[2]; sdthr = param[3]; sdspread = param[4]; sdslope = param[5]; printf("\nBOOTSTRAP ESTIMATES OF STANDARD DEVIATIONS AND CONFIDENCE " "INTERVALS\n"); printf("(based on %d bootstrap replications)\n",nboot); printf(" SD"); printf(" 95%% CIs 90%% CIs\n"); printf("threshold: %9.4f", sdthr); printf(" [%9.4f,%9.4f] [%9.4f,%9.4f]\n", thrlimits[0], thrlimits[1], thrlimits[2], thrlimits[3]); printf("gradient: %9.4f", sdslope); printf(" [%9.4f,%9.4f] [%9.4f,%9.4f]\n", slopelimits[0], slopelimits[1], slopelimits[2], slopelimits[3]); printf("1/gradient: %9.4f", sdspread); printf(" [%9.4f,%9.4f] [%9.4f,%9.4f]\n", spreadlimits[0], spreadlimits[1], spreadlimits[2], spreadlimits[3]); /* DIAGNOSTICS ONLY printf("\nPress any key to continue ..."); (void)getch(); printf("\n\nDiagnostics:\n"); printf(" Bias\n"); printf("thresh: %9.4f\n", thr-thr0); printf("spread: %9.4f\n", spread-spread0); printf("slope: %9.4f\n", slope-slope0); printf("Failed fits during bootstrap replications = %d\n", ifails); printf("Check mean random number (0.5) = %f\n", meanran/count); */ return 0; } /* ----------------------------------------------------------------------- */ int probit(int nlevel, float slevel[], int trials[], int hits[], float pcrit, float *pa, float *pb, float *pchi, int *pdof, float param[]) { float p[NLEVEL]; /* empirical probabilities */ float prob[NLEVEL]; /* empirical probits */ float weight[NLEVEL]; /* corresponding weights */ float ep[NLEVEL]; /* exp. probabilities */ float eprob[NLEVEL]; /* exp. probits */ float a, b, chi; float othr, ostd, sthr, std; int dof; int i; int iterations, cycles; float z; int overflow; othr = -MAXFLOAT; ostd = -MAXFLOAT; iterations = 50; /* normally need 3-4 */ for (i = 0; i < nlevel; i++) { weight[i] = 1.0; /* initially unweighted regression */ p[i] = (float)hits[i]/(float)trials[i]; /* prop 1 responses - p(1) */ prob[i] = invcgauss(p[i]); /* get probit of p(1) */ } overflow = regress(nlevel, slevel, prob, weight, pcrit, &a, &b, &chi, &dof, param); if( overflow ) { b = 0.0; /* If overflow, gradient = 0.0, and */ a = 0.0; /* intercept = 0.0 */ } for (cycles = 0; cycles < iterations; cycles++) { for (i = 0; i < nlevel; i++) { eprob[i] = a + b*slevel[i]; ep[i] = cgauss(eprob[i]); z = ord(eprob[i]); weight[i] = z*z*(float)trials[i]/(ep[i]*(1.0-ep[i])); eprob[i] = eprob[i] + (p[i]-ep[i])/z; if (eprob[i] > MAXPROB) eprob[i] = MAXPROB; /* trap excess values */ if (eprob[i] < -MAXPROB) eprob[i] = -MAXPROB; } overflow = regress(nlevel, slevel, eprob, weight, pcrit, &a, &b, &chi, &dof, param); if ( !overflow ) { /* if ok result, see if converging */ sthr = param[0]; std = param[1]; if (fabs(sthr-othr) < EPS_C && fabs(std-ostd) < EPS_C) { overflow = 0; /* converged flag */ break; } else { /* not converged, so update values */ othr = sthr; ostd = std; overflow = 1; /* unconverged flag */ } } else { /* printf(" **OVERFLOW**\n"); suppress regress overflow message */ /* printf("cycle=%d, chi=%f a=%f b=%f \n", cycles, chi, a, b); */ } } *pa = a; *pb = b; *pchi = chi; *pdof = dof; return(overflow); } /* ----------------------------------------------------------------------- */ float ord(float z) { if (fabs(z) < MAXPROB) return( exp(-z*z/2.0) / 2.506628 ); else return( exp(-fsqr(MAXPROB)/2.0) /2.506628 ); } /* ----------------------------------------------------------------------- */ float fsqr(float x) { return(x*x); } /* ----------------------------------------------------------------------- */ int regress(int n, float x[], float y[], float w[], float pcrit, float *pa, float *pb, float *pchi, int *pdof, float param[]) { /* This routine performs weighted linear regression on up to NLEVEL points. Arguments are identified below. */ /* n number of points to be fitted x[NLEVEL], y[NLEVEL]; x, y coordinates of each point w[NLEVEL]; weights on each point a intercept b coefficient chi chi-square value for fit dof degrees of freedom param[0] thr param[1] spread param[2] estimated sd of thr param[3] estimated variance of thr param[4] estiamted sd of gradient b param[5] estimated variance of gradient b */ int i; /* loop counter */ int dof; /* dof */ float sw = 0.0; /* sum of weights */ float swx = 0.0, swy = 0.0; /* sum of weighted x, y, values */ float swxx = 0.0, swyy = 0.0; /* sum of weighted x^2, y^2 values */ float swxy = 0.0; /* sum of weighted x*y values */ float sx, sy, sxx, sxy, syy; /* working variables (WFB) */ float a, b, chi, ww; for(i = 0; i < n; i++) { ww = w[i]; sw += ww; swx += x[i] * ww; swy += y[i] * ww; swxx += x[i] * x[i] * ww; swyy += y[i] * y[i] * ww; swxy += x[i] * y[i] * ww; } if(sw > EPS_Q) { sx = swx/sw; sy = swy/sw; sxx = swxx - swx*swx/sw; sxy = swxy - swx*swy/sw; syy = swyy - swy*swy/sw; b = sxy/sxx; a = sy - b*sx; chi = syy - sxy*sxy/sxx; dof = n - 2; if (b*b > EPS_Q) { param[0] = (invcgauss(pcrit)-a)/b; param[1] = 1.0/b; param[3] = (1.0/sw + fsqr(param[0]-sx)/sxx)/(b*b); param[2] = sqrt(fabs(param[3])); param[5] = 1.0/sxx; param[4] = sqrt(fabs(param[5])); *pa = a; *pb = b; *pchi = chi; *pdof = dof; return(0); /* successful regression */ } return(0); } else return(1); /* overflow in regression */ } /* ----------------------------------------------------------------------- */ void sort(int n, float *r) /* Sort array r of length n using shell sort. */ { float temp; int i,j,gap; for (gap=1; gap<=n; gap=3*gap+1); for (; gap>0; gap/=3) { for (i=gap; i=0) && (r[j] > temp); j-=gap) r[j+gap]=r[j]; r[j+gap]=temp; } } } /* ----------------------------------------------------------------------- */ float centiles(int nboot, float ra[], float limits[]) { /* Sorts bootstrap samples in array. Returns either standard estimate of the SD, which is susceptible to outliers, or robust but biased estimate based on 15.87% and 84.13% centiles. */ float lower, upper, sd, fboot, mb, ssb, i; fboot = (float)nboot; sort(nboot,ra); /* sort into ascending order for centiles */ lower = ra[(int)(0.1587*fboot)]; upper = ra[(int)(0.8413*fboot)]; /* Robust biased normal estimate of SD */ sd = (upper-lower)*0.5; /* Standard estimate of SD*/ //mb = 0; //ssb = 0; //for (i = 0; i < nboot; i++) { // mb = mb + ra[i]; // ssb = ssb + ra[i]*ra[i]; //} //mb = mb/fboot; //sd = sqrt(ssb/fboot - mb*mb); /* printf("\nlower = %f upper = %f sd = %f\n",lower,upper,sd); */ limits[0] = ra[(int)(0.025*fboot)]; limits[1] = ra[(int)(0.975*fboot)]; limits[2] = ra[(int)(0.05*fboot)]; limits[3] = ra[(int)(0.95*fboot)]; return(sd); } /* ----------------------------------------------------------------------- */ void boot(int nlevel, float slevel[], int trials[], int hits[], float pcrit, float a, float b, int nboot, float param[], float thrlimits[], float spreadlimits[], float slopelimits[]) { int hitsboot[NLEVEL]; float pfit[NLEVEL]; float bootthr[NBOOT], bootspread[NBOOT], bootslope[NBOOT]; float thr, sumthr, slope, sumslope, spread, sumspread; float fboot; float chi; int i; int iboot; int overflow; int dof; sumthr = sumspread = sumslope = 0.0; printf("Now bootstrapping ...\n"); for (i = 0; i < nlevel; i++) { pfit[i] = cgauss(a+b*slevel[i]); /* parametric bootstrap */ /* pfit[i] = (float)hits[i]/(float)trials[i]; nonparametric bootstrap */ } for (iboot = 0; iboot < nboot; iboot++) { for (i = 0; i < nlevel; i++) { hitsboot[i] = binom_rand(pfit[i], trials[i]); } overflow = probit(nlevel,slevel,trials,hitsboot,pcrit,&a,&b,&chi, &dof,param); if (chi > 10.0*(float)nlevel || overflow) { /* trap badly fitting data sets */ /* printf("%d %f %f\n", overflow, param[0], param[1]); */ ifails++; iboot--; continue; } /* if ((float)iboot/20.0 - (float)(iboot/20) < 0.00001) printf("\n"); */ thr = bootthr[iboot] = param[0]; spread = bootspread[iboot] = param[1]; slope = bootslope[iboot] = 1.0/spread; // CARE //printf("%f\n", slope); //if (slope > 7) { // for (i = 0; i < nlevel; i++) { // printf("%d %d\n", trials[i], hitsboot[i]); // } // getch(); // exit(1); //} sumthr += thr; /* as a check on bias */ sumspread += spread; sumslope += slope; } fboot = (float)nboot; param[0] = sumthr/fboot; /* mean of bootstrap sample */ param[1] = sumspread/fboot; param[2] = sumslope/fboot; param[3] = centiles(nboot, bootthr, thrlimits); /* robust estimator of sd */ param[4] = centiles(nboot, bootspread, spreadlimits); param[5] = centiles(nboot, bootslope, slopelimits); } /* ----------------------------------------------------------------------- */ int degen(int nlevel, int trials[], int hits[]) { /* Check for degenerate data */ float p[NLEVEL]; int reject, i, ih, il, m1; /* reject (0 0 0 0 0) and (0 0 0 0 p) */ il = 0; while (hits[il] <= 0 && il < nlevel) il++; il++; if (il >= nlevel) { if (il > nlevel) reject = 2; else reject = 1; return(reject); } /* reject (1 1 1 1 1) and (p 1 1 1 1) */ ih = nlevel-1; while (hits[ih] >= trials[ih] && ih > -1) ih--; ih++; if (ih <= 1) { if (ih < 1) reject = 2; else reject = 1; return(reject); } /* reject (0 0 p 1 1) and (0 0 1 1 1) */ if (il >= ih) { reject = 3; return(reject); } /* reject symmetric samples */ for (i = 0; i < nlevel; i++) p[i] = (float)hits[i]/(float)trials[i]; m1 = nlevel/2; il = 0; while (fabs(p[il] - p[nlevel-il]) < 0.0001 && il < m1) il++; if (il == m1) reject = 4; else reject = 0; return(reject); } /* ----------------------------------------------------------------------- */ float invcgauss(float p) { /* Inverse cumulative Gaussian transform. This routine takes a proportion p and returns the corresponding x-value of a cumulative Gaussian of zero mean and unit standard deviation. Odeh, R. E. & Evans, J. O. (1974) Algorithm AS 70: Percentage points of the normal distribution. Applied Statistics, 23, 96-97. Described in Kennedy, W. J. & Gentle, J. E. (1980) Statistical Computing, New York: Dekker, pp. 95-96. */ float lim=EPS_Q; float p0=-0.322232431088, p1=-1.0, p2=-0.342242088547; float p3=-0.0204231210245, p4=-0.453642210148E-4; float q0=0.0993484626060, q1=0.588581570495; float q2=0.531103462366, q3=0.103537752850, q4=0.38560700634E-2; float xp, y, ptrue; ptrue = p; /* store value */ if (p > 0.5) p = 1.0 - p; if (p < lim) p = lim; /* minimum value of p */ if (p == 0.5) return(0.0); y = sqrt(log(1.0/(p*p))); xp = y+((((y*p4+p3)*y+p2)*y+p1)*y+p0)/((((y*q4+q3)*y+q2)*y+q1)*y+q0); return(ptrue < 0.5 ? -xp : xp); } /* ----------------------------------------------------------------------- */ float cgauss(float x) { /* Cumulative Gaussian. This routine takes a number between -infinity and infinity and returns the corresponding p-value of a cumulative Gaussian, the underlying Gaussian having zero mean and unit standard deviation. */ if (x > MAXPROB) x = MAXPROB; /* trap excessive value */ if (x < -MAXPROB) x = -MAXPROB; if (x < 0.0) return(0.5*erfcc(-x/SQRT2)); else return(1.0 - 0.5*erfcc(x/SQRT2)); } /* ----------------------------------------------------------------------- */ float erfcc(float x) { /* Error function. Source: Cody, W. J. (1969). Rational Chebyshev approximations for the error function, Math. Comp., 23, 631-638. Described in Kennedy, W. J. & Gentle, J. E. (1980) Statistical Computing, New York: Dekker, pp. 90ff */ static float p1[4] = { 2.4266795523053175e2,2.1979261618294152e1,6.9963834886191355, -3.5609843701815385e-2}; static float q1[4] = { 2.1505887586986120e2,9.1164905404514901e1,1.5082797630407787e1, 1.0}; static float p2[8] = { 3.004592610201616005e2,4.519189537118729422e2,3.393208167343436870e2, 1.529892850469404039e2,4.316222722205673530e1,7.211758250883093659, 5.641955174789739711e-1,-1.368648573827167067e-7}; static float q2[8] = { 3.004592609569832933e2,7.909509253278980272e2,9.313540948506096211e2, 6.389802644656311665e2,2.775854447439876434e2,7.700015293522947295e1, 1.278272731962942351e1,1.0}; static float p3[5] = { -2.99610707703542174e-3,-4.94730910623250734e-2,-2.26956593539686930e-1, -2.78661308609647788e-1,-2.23192459734184686e-2}; static float q3[5] = { 1.06209230528467918e-2,1.91308926107829841e-1,1.05167510706793207, 1.98733201817135256,1.0}; float z,zz,zpow,p,q,r; int j; z=fabs(x); if (z < 0.46875) { zz=z*z; zpow=zz; p=p1[0]; q=q1[0]; for (j=1;j<=3;j++) { p+=p1[j]*zpow; q+=q1[j]*zpow; zpow*=zz; } r=1-z*p/q; } else if (z < 4.0) { zpow=z; p=p2[0]; q=q2[0]; for (j=1;j<=7;j++) { p+=p2[j]*zpow; q+=q2[j]*zpow; zpow*=z; } r=exp((double)(-z*z))*p/q; } else { zz=1.0/(z*z); zpow=zz; p=p3[0]; q=q3[0]; for (j=1;j<=4;j++) { p+=p3[j]*zpow; q+=q3[j]*zpow; zpow*=zz; } r=exp((double)(-zz))/z*(0.564189583548+zz*p/q); } return r; } /* ------------------------------------------------------------------------ */ int binom_rand(float pp,int nn) { /* Generates a single random deviate from a binomial distribution whose number of trials is N and whose probability of an event in each trial is P. This is algorithm BTPE from: Kachitvichyanukul, V. and Schmeiser, B. W. Binomial Random Variate Generation. Communications of the ACM, 31, 2 (February, 1988) 216. Source: ranlib.c at www.netlib.org */ static float psave = -1.0; static long nsave = -1; static long ignbin,i,ix,ix1,k,m,mp,n,T1; static float al,alv,amaxp,c,f,f1,f2,ffm,fm,g,p,p1,p2,p3,p4,q,qn,r,u,v,w,w2,x,x1, x2,xl,xll,xlr,xm,xnp,xnpq,xr,ynorm,z,z2; n=nn; if(pp != psave) goto S10; if(n != nsave) goto S20; if(xnp < 30.0) goto S150; goto S30; S10: /* *****SETUP, PERFORM ONLY WHEN PARAMETERS CHANGE */ psave = pp; p = MIN(psave,1.0-psave); q = 1.0-p; S20: xnp = n*p; nsave = n; if(xnp < 30.0) goto S140; ffm = xnp+p; m = ffm; fm = m; xnpq = xnp*q; p1 = (long) (2.195*sqrt(xnpq)-4.6*q)+0.5; xm = fm+0.5; xl = xm-p1; xr = xm+p1; c = 0.134+20.5/(15.3+fm); al = (ffm-xl)/(ffm-xl*p); xll = al*(1.0+0.5*al); al = (xr-ffm)/(xr*q); xlr = al*(1.0+0.5*al); p2 = p1*(1.0+c+c); p3 = p2+c/xll; p4 = p3+c/xlr; S30: /* *****GENERATE VARIATE */ u = drand48()*p4; v = drand48(); /* TRIANGULAR REGION */ if(u > p1) goto S40; ix = xm-p1*v+u; goto S170; S40: /* PARALLELOGRAM REGION */ if(u > p2) goto S50; x = xl+(u-p1)/c; v = v*c+1.0-ABS(xm-x)/p1; if(v > 1.0 || v <= 0.0) goto S30; ix = x; goto S70; S50: /* LEFT TAIL */ if(u > p3) goto S60; ix = xl+log(v)/xll; if(ix < 0) goto S30; v *= ((u-p2)*xll); goto S70; S60: /* RIGHT TAIL */ ix = xr-log(v)/xlr; if(ix > n) goto S30; v *= ((u-p3)*xlr); S70: /* *****DETERMINE APPROPRIATE WAY TO PERFORM ACCEPT/REJECT TEST */ k = ABS(ix-m); if(k > 20 && k < xnpq/2-1) goto S130; /* EXPLICIT EVALUATION */ f = 1.0; r = p/q; g = (n+1)*r; T1 = m-ix; if(T1 < 0) goto S80; else if(T1 == 0) goto S120; else goto S100; S80: mp = m+1; for(i=mp; i<=ix; i++) f *= (g/i-r); goto S120; S100: ix1 = ix+1; for(i=ix1; i<=m; i++) f /= (g/i-r); S120: if(v <= f) goto S170; goto S30; S130: /* SQUEEZING USING UPPER AND LOWER BOUNDS ON ALOG(F(X)) */ amaxp = k/xnpq*((k*(k/3.0+0.625)+0.1666666666666)/xnpq+0.5); ynorm = -(k*k/(2.0*xnpq)); alv = log(v); if(alv < ynorm-amaxp) goto S170; if(alv > ynorm+amaxp) goto S30; /* STIRLING'S FORMULA TO MACHINE ACCURACY FOR THE FINAL ACCEPTANCE/REJECTION TEST */ x1 = ix+1.0; f1 = fm+1.0; z = n+1.0-fm; w = n-ix+1.0; z2 = z*z; x2 = x1*x1; f2 = f1*f1; w2 = w*w; if(alv <= xm*log(f1/x1)+(n-m+0.5)*log(z/w)+(ix-m)*log(w*p/(x1*q))+(13860.0- (462.0-(132.0-(99.0-140.0/f2)/f2)/f2)/f2)/f1/166320.0+(13860.0-(462.0- (132.0-(99.0-140.0/z2)/z2)/z2)/z2)/z/166320.0+(13860.0-(462.0-(132.0- (99.0-140.0/x2)/x2)/x2)/x2)/x1/166320.0+(13860.0-(462.0-(132.0-(99.0 -140.0/w2)/w2)/w2)/w2)/w/166320.0) goto S170; goto S30; S140: /* INVERSE CDF LOGIC FOR MEAN LESS THAN 30 */ qn = pow(q,(double)n); r = p/q; g = r*(n+1); S150: ix = 0; f = qn; u = drand48(); S160: if(u < f) goto S170; if(ix > 110) goto S150; u -= f; ix += 1; f *= (g/ix-r); goto S160; S170: if(psave >= 0.5) ix = n-ix; ignbin = ix; return (int)ignbin; } /* -------------------------------------------------------------------------- * The random number routines were adapted from the bsd sources, that can be * obtained from e.g. * ftp://gatekeeper.dec.com/pub/BSD/FreeBSD/FreeBSD-current/src/lib/libc/gen * and come with the following copyright: * * Copyright (c) 1993 Martin Birgmeier * All rights reserved. * * You may redistribute unmodified or modified versions of this source * code provided that the above copyright notice and this and the * following conditions are retained. * * This software is provided ``as is'', and comes with no warranties * of any kind. I shall in no event be liable for anything that happens * to anyone/anything when using this software. * -------------------------------------------------------------------------- */ #define RAND48_SEED_0 (0x330e) #define RAND48_SEED_1 (0xabcd) #define RAND48_SEED_2 (0x1234) #define RAND48_MULT_0 (0xe66d) #define RAND48_MULT_1 (0xdeec) #define RAND48_MULT_2 (0x0005) #define RAND48_ADD (0x000b) unsigned short _rand48_seed[3] = { RAND48_SEED_0, RAND48_SEED_1, RAND48_SEED_2 }; unsigned short _rand48_mult[3] = { RAND48_MULT_0, RAND48_MULT_1, RAND48_MULT_2 }; unsigned short _rand48_add = RAND48_ADD; void _dorand48(unsigned short xseed[3]) { unsigned long accu; unsigned short temp[2]; accu = (unsigned long) _rand48_mult[0] * (unsigned long) xseed[0] + (unsigned long) _rand48_add; temp[0] = (unsigned short) accu; /* lower 16 bits */ accu >>= sizeof(unsigned short) * 8; accu += (unsigned long) _rand48_mult[0] * (unsigned long) xseed[1] + (unsigned long) _rand48_mult[1] * (unsigned long) xseed[0]; temp[1] = (unsigned short) accu; /* middle 16 bits */ accu >>= sizeof(unsigned short) * 8; accu += _rand48_mult[0] * xseed[2] + _rand48_mult[1] * xseed[1] + _rand48_mult[2] * xseed[0]; xseed[0] = temp[0]; xseed[1] = temp[1]; xseed[2] = (unsigned short) accu; } double erand48(unsigned short xseed[3]) { _dorand48(xseed); return ldexp((double) xseed[0], -48) + ldexp((double) xseed[1], -32) + ldexp((double) xseed[2], -16); } double drand48(void) { return erand48(_rand48_seed); } void srand48(long seed) { _rand48_seed[0] = RAND48_SEED_0; _rand48_seed[1] = (unsigned short) seed; _rand48_seed[2] = (unsigned short) (seed >> 16); _rand48_mult[0] = RAND48_MULT_0; _rand48_mult[1] = RAND48_MULT_1; _rand48_mult[2] = RAND48_MULT_2; _rand48_add = RAND48_ADD; } ========================test1.dat==CUT HERE============================== 5 0.5 -2.0 5 0 -1.0 5 1 0.0 5 3 1.0 5 4 2.0 5 5 test1.dat test data The first row has the number of levels (5) and the criterion level of performance (0.5), and the second and subsequent rows give the stimulus level (e.g. -2.0), the number of trials at that level (e.g. 5) and the number of hits or successes (e.g. 0) at that level. Any text may be written after the last line of data. bootprob yields (with 1000 replications): Estimate threshold: -0.0886 gradient: 1.0433 1/gradient: 0.9585 SD 95% CIs 90% CIs threshold: 0.3373 [ -0.8207, 0.6877] [ -0.7328, 0.5000] gradient: 1.7005 [ 0.6864, 4.8213] [ 0.7059, 4.8213] 1/gradient: 0.4244 [ 0.2074, 1.4569] [ 0.2074, 1.4167] ========================test1b.dat==CUT HERE============================= 5 0.5 8.0 5 0 9.0 5 1 10.0 5 3 11.0 5 4 12.0 5 5 test1b.dat test data The first row has the number of levels (5) and the criterion level of performance (0.5), and the second and subsequent rows give the stimulus level (e.g. -2.0), the number of trials at that level (e.g. 5) and the number of hits or successes (e.g. 0) at that level. Any text may be written after the last line of data. bootprob on test1b.dat yields (with 1000 replications): Estimate threshold: 9.9114 gradient: 1.0433 1/gradient: 0.9585 SD 95% CIs 90% CIs threshold: 0.3398 [ 9.1806, 10.5587] [ 9.2499, 10.5000] gradient: 1.6158 [ 0.6418, 4.7879] [ 0.7059, 4.7740] 1/gradient: 0.4402 [ 0.2089, 1.5582] [ 0.2095, 1.4167] bootprobdiff on test1.dat and test1b.dat yields (with 1000 replications): SD 95% CIs 90% CIs threshold difference: 0.5088 [ 8.9790, 10.9930] [ 9.1284, 10.8129] gradient difference: 1.8808 [ -3.8935, 4.0153] [ -3.5935, 3.5147] 1/gradient difference: 0.5224 [ -1.0132, 1.0414] [ -0.7496, 0.9035] ========================test2.dat==CUT HERE============================== 11 0.5 0.0 100 1 10.0 100 2 20.0 100 5 30.0 100 12 40.0 100 27 50.0 100 50 60.0 100 73 70.0 100 88 80.0 100 95 90.0 100 98 100.0 100 99 test2.dat test data The first row has the number of levels (5) and the criterion level of performance (0.5), and the second and subsequent rows give the stimulus level (e.g. -2.0), the number of trials at that level (e.g. 5) and the number of hits or successes (e.g. 0) at that level. Any text may be written after the last line of data. bootprob yields (with 1000 replications): Estimate threshold: 50.0000 gradient: 0.0538 1/gradient: 18.5877 SD 95% CIs 90% CIs threshold: 0.9788 [ 47.9460, 51.9792] [ 48.4116, 51.6328] gradient: 0.0027 [ 0.0490, 0.0599] [ 0.0498, 0.0587] 1/gradient: 0.9131 [ 16.6977, 20.3897] [ 17.0462, 20.0807]