home *** CD-ROM | disk | FTP | other *** search
/ Tricks of the Mac Game Programming Gurus / TricksOfTheMacGameProgrammingGurus.iso / More Source / Libraries / VideoToolbox 95.04.18 / Utilities / Quick3 / PsychometricFit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-09-19  |  6.2 KB  |  169 lines  |  [TEXT/MMCC]
  1. /* 
  2. PsychometricFit.c
  3. Copyright 1990 (c) Denis G. Pelli
  4. A general-purpose function that does a maximum likelihood fit of any
  5. psychometric function to psychometric data. The returned value is the level of
  6. significance at which the fit can be rejected. The degreesOfFreedom may be zero,
  7. in which case no parameters will be adjusted, but you'll get the log likelihood
  8. and significance of the supplied parameter values.
  9.  
  10. The psychometric function (which you supply as an argument) takes two arguments:
  11. a contrast and a pointer to a paramRecord. The function Weibull() is provided in
  12. Weibull.c. Others may be added as desired. It is assumed that the first
  13. parameter is the log of threshold. No assumptions are made about the other
  14. parameters, except that any that are to be iteratively fit are assumed to be of
  15. type "double".
  16.  
  17. Quick3.c is a stand-alone program that uses PsychometricFit() to do the real work.
  18. I suggest you read the introductory comments at the beginning of Quick3.c
  19.  
  20. HISTORY:
  21. 4/7/90        dgp wrote it
  22. 10/29/90    dgp    tidied up the comments
  23. 11/17/92    dgp "
  24. 1/25/93 dgp removed obsolete support for THINK C 4.
  25. 9/5/94 dgp removed assumption in printf's that int==short.
  26.  
  27. SOURCES:
  28. Quick3.h
  29. LogLikelihood.c
  30. MonotonicFit.c
  31. PsychometricFit.c
  32. SortAndMergeContrasts.c
  33. Weibull.c
  34. #From Denis Pelli's VideoToolbox:
  35. VideoToolbox.h
  36. Binomial.c
  37. ChiSquare.c
  38. Normal.c
  39. SetFileInfo.c        # Used only on the Macintosh
  40. #From Numerical Recipes in C:
  41. nr.h
  42. NRUTIL.h
  43. BRENT.C
  44. F1DIM.C
  45. LINMIN.C
  46. MNBRAK.C
  47. NRUTIL.C
  48. POWELL.C
  49.  
  50. LIMITATIONS
  51.  
  52. This program uses routines from Numerical Recipes in C. They're copyrighted, 
  53. so I can't distribute them. You can order the software:
  54.     Numerical Recipes C Diskette for Macintosh $29.95
  55. and book:
  56.     Numerical Recipes in C: The Art of Scientific Computing $44.50
  57. from:
  58.     Cambridge University Press
  59.     Order Department
  60.     510 North Avenue
  61.     New Rochelle, NY 10801
  62.  
  63. Note that I have made several changes to the Numerical Recipes in C routines: 
  64. 1.In every file I changed "float" to "FLOAT", and #included nr.h. I inserted the 
  65. statement "typedef double FLOAT;" in the file nr.h. This is because the 
  66. Macintosh computes doubles much faster than floats. If you'd rather run
  67. slowly than modify your Numerical Recipes in C files, then you will need to insert 
  68. "typedef float FLOAT;"
  69. in CalibrateLuminance.c & PsychometricFit.c in order to compile those files.
  70. The rest of the VideoToolbox doesn't care.
  71. */
  72. #include "VideoToolbox.h"
  73. #include "Quick3.h"
  74. #include "nr.h"
  75. #include "nrutil.h"
  76.  
  77. #define TOLERANCE 0.001    /* fractional tolerance of log likelihood. Not critical. */
  78.  
  79. /*
  80. I hate global variables because they hide the flow of information. However,
  81. some sort of cludge is necessary to pass the extra arguments to Error(), bypassing
  82. the Numerical Recipes routines that call it, since they only pass the 
  83. parameters they know about. These static declarations at least restrict the scope of
  84. these "globals" to this file.
  85. */
  86. static double Error(FLOAT *p);
  87. static dataRecord *myDataPtr;                /* for Error() */
  88. static PsychometricFunctionPtr MyPsychFun;    /* for Error() */
  89. static paramRecord myParams;                /* for Error() */
  90. static int myDegreesOfFreedom;                /* for Error() */
  91. static int iter;                            /* for Error() */
  92.  
  93. double PsychometricFit(paramRecord *paramPtr,PsychometricFunctionPtr PsychFun
  94.     ,dataRecord *dataPtr,double *logLikelihoodPtr,int degreesOfFreedom
  95.     ,double *chiSquarePtr,int *chiSquareDFPtr)
  96. {
  97.     int i,j;
  98.     FLOAT *p,**direction,ftol,fret;
  99.     dataRecord monotonicData;
  100.     double monotonicLL;
  101.     int monotonicDF;
  102.     double P;
  103.     
  104.     myDataPtr=dataPtr;        /* copy these for use by Error() */
  105.     MyPsychFun=PsychFun;
  106.     myParams=*paramPtr;
  107.     myDegreesOfFreedom=degreesOfFreedom;
  108.     
  109.     p=vector(1,degreesOfFreedom);
  110.     direction=matrix(1,degreesOfFreedom,1,degreesOfFreedom);    /* initial directions */
  111.     if(p==NULL || direction == NULL)
  112.         PrintfExit("PsychometricFit: not enough room for arrays.\007\n");
  113.     for(i=1;i<=degreesOfFreedom;i++) p[i]=((double *)paramPtr)[i-1];
  114.     for(i=1;i<=degreesOfFreedom;i++)for(j=1;j<=degreesOfFreedom;j++)direction[i][j]=0.0;
  115.     for(i=1;i<=degreesOfFreedom;i++)direction[i][i]=0.03;    /* initial step size */
  116.     ftol=TOLERANCE;    /* fractional tolerance on Error value when done */
  117.     iter=0;
  118.     
  119.     /* do it. The psychometric function is passed to Error by the global MyPsychFun */
  120.     if(degreesOfFreedom==0)fret=Error(p);
  121.     else powell(p,direction,degreesOfFreedom,ftol,&iter,&fret,&Error);
  122.  
  123.     for(i=1;i<=degreesOfFreedom;i++) ((double *)paramPtr)[i-1]=p[i];
  124.     free_matrix(direction,1,degreesOfFreedom,1,degreesOfFreedom);
  125.     free_vector(p,1,degreesOfFreedom);
  126.  
  127.     *logLikelihoodPtr=-fret;
  128.     
  129.     /* Now compute the degree of significance at which the fit can be rejected */
  130.     monotonicData= *dataPtr;
  131.     MonotonicFit(&monotonicData,&monotonicLL,&monotonicDF);    /* overwrites data with fit */
  132.     *chiSquarePtr= -2.0*(*logLikelihoodPtr-monotonicLL);    /* -2 log likelihood ratio of hypotheses */
  133.     *chiSquareDFPtr=monotonicDF-degreesOfFreedom;            /* difference in degrees of freedom */
  134.     P=PChiSquare(*chiSquarePtr,*chiSquareDFPtr);            /* significance */
  135.     return P;
  136. }
  137.  
  138. /* There is a subtlety here. I thought that I could use Powell with the whole
  139. paramRecord, yet ask Powell to only twiddle the first few parameters, figuring that
  140. even when I was asking Powell to fit only the first few parameters the other
  141. parameters would still be there in the array when the pointer to the array was
  142. passed to Error(). Alas, Powell() and its subroutines make COPIES of the array,
  143. and naturally fail to copy the non-twiddled parameters, since they don't know about
  144. them. The solution is for Error() to have its own complete copy of the paramRecord.
  145. Each time Error() is called it updates the twiddled parameters before calling
  146. LogLikelihood(), which calls the psychometric function (*MyPsychFun)().
  147. */
  148.  
  149. double Error(FLOAT *p)
  150. {
  151.     double error;
  152.     int i;
  153.     static int lastIter=0;
  154.     
  155.     for(i=1;i<=myDegreesOfFreedom;i++) ((double *)&myParams)[i-1]=p[i];
  156.     
  157.     error=-LogLikelihood(myDataPtr,&myParams,MyPsychFun);
  158.     
  159.     /* Diagnostic printout for difficult cases */
  160.     if(iter>0 && iter%50 == 0 && iter!=lastIter){
  161.         printf("Error(): Warning, %d iterations:\n",(int)iter);
  162.         printf("logAlpha %5.2f, beta %5.1f, gamma %5.2f, delta %6.3f -log likelihood %9.0g\n"
  163.             ,myParams.logAlpha,myParams.beta,myParams.gamma,myParams.delta,error);
  164.         lastIter=iter;
  165.     }
  166.     return error;
  167. }
  168.  
  169.