/* * * * Simple Good-Turing Frequency Estimator * * * Geoffrey Sampson, with help from Miles Dennis * * Department of Informatics * Sussex University * * www.grsampson.net * * * First release: 27 June 1995 * Revised release: 24 July 2000 * This header information revised: 23 March 2005 * Further revised release: 8 April 2008 * Further revised release: 11 July 2008 * * * Takes a set of (frequency, frequency-of-frequency) pairs, and * applies the "Simple Good-Turing" technique for estimating * the probabilities corresponding to the observed frequencies, * and P.0, the joint probability of all unobserved species. * The Simple Good-Turing technique was devised by the late William * A. Gale of AT&T Bell Labs, and described in Gale & Sampson, * "Good-Turing Frequency Estimation Without Tears" (JOURNAL * OF QUANTITATIVE LINGUISTICS, vol. 2, pp. 217-37 -- reprinted in * Geoffrey Sampson, EMPIRICAL LINGUISTICS, Continuum, 2001). * * Anyone is welcome to take copies of this program and use it * for any purpose, at his or her own risk. If it is used in * connexion with published work, acknowledgment of Sampson and * the University of Sussex would be a welcome courtesy. * * The program is written to take input from "stdin" and send output * to "stdout"; redirection can be used to take input from and * send output to permanent files. The code is in ANSI standard C. * * The input file should be a series of lines separated by newline * characters, where all nonblank lines contain two positive integers * (an observed frequency, followed by the frequency of that frequency) * separated by whitespace. (Blank lines are ignored.) * The lines should be in ascending order of frequency, and must * begin with frequency 1. * * No checks are made for linearity; the program simply assumes that the * requirements for using the SGT estimator are met. * * The output is a series of lines each containing an integer followed * by a probability (a real number between zero and one), separated by a * tab. In the first line, the integer is 0 and the real number is the * estimate for P.0. In subsequent lines, the integers are the * successive observed frequencies, and the reals are the estimated * probabilities corresponding to those frequencies. * * Later releases cure bugs to which my attention has kindly been * drawn at different times by Martin Jansche of Ohio State University * and Steve Arons of New York City. No warranty is given * as to absence of further bugs. * * Fan Yang of Next IT Inc., Spokane, Washington, has suggested to me * that in the light of his experience with the SGT technique, for some * data-sets it could be preferable to use the 0.1 significance criterion * actually used in the experiments reported in the Gale & Sampson * paper, rather than the 0.05 criterion suggested in that paper * for the sake of conformity with standard statistical convention. * (See note 8 of the paper.) Neither Fan Yang nor I have pursued * this far enough to formulate a definite recommendation; but, in * order to make it easier for users of the software to experiment * with alternative confidence levels, the July 2008 release moves * the relevant "magic number" out of the middle of the program into * a #define line near the beginning where it is given the constant * name CONFID_FACTOR. The value 1.96 corresponds to the p < 0.05 * criterion; in order to use the p < 0.1 criterion, 1.96 in the * #define line should be changed to 1.65. * * */ #include #include #include #include #include #define TRUE 1 #define FALSE 0 #define MAX_LINE 100 #define MAX_ROWS 200 #define MIN_INPUT 5 #define CONFID_FACTOR 1.96 int r[MAX_ROWS], n[MAX_ROWS]; double Z[MAX_ROWS], log_r[MAX_ROWS], log_Z[MAX_ROWS], rStar[MAX_ROWS], p[MAX_ROWS]; int rows, bigN; double PZero, bigNprime, slope, intercept; int main(void) { int readValidInput(void); void analyseInput(void); if ((rows = readValidInput()) >= 0) { if (rows < MIN_INPUT) printf("\nFewer than %d input value-pairs\n", MIN_INPUT); else analyseInput(); } return(TRUE); } double sq(double x) { return(x * x); } int readValidInput(void) /* * returns number of rows if input file is valid, else -1 * NB: number of rows is one more than index of last row * */ { char line[MAX_LINE]; const char* whiteSpace = " \t\n\v\f\r"; int lineNumber = 0; int rowNumber = 0; const int error = -1; while (fgets(line, MAX_LINE, stdin) != NULL && rowNumber < MAX_ROWS) { char* ptr = line; char* integer; int i; ++lineNumber; while (isspace(*ptr)) ++ptr; /* skip white space at the start of a line */ if (*ptr == '\0') continue; if ((integer = strtok(ptr, whiteSpace)) == NULL || (i = atoi(integer)) < 1) { fprintf(stderr, "Invalid field 1, line %d\n", lineNumber); return(error); } if (rowNumber > 0 && i <= r[rowNumber - 1]) { fprintf(stderr, "Frequency not in ascending order, line %d\n", lineNumber); return(error); } r[rowNumber] = i; if ((integer = strtok(NULL, whiteSpace)) == NULL || (i = atoi(integer)) < 1) { fprintf(stderr, "Invalid field 2, line %d\n", lineNumber); return(error); } n[rowNumber] = i; if (strtok(NULL, whiteSpace) != NULL) { fprintf(stderr, "Invalid extra field, line %d\n", lineNumber); return(error); } ++rowNumber; } if (rowNumber >= MAX_ROWS) { fprintf(stderr, "\nInsufficient memory reserved for input\ values\nYou need to change the definition of\ MAX_ROWS\n"); return(error); } return(rowNumber); } void findBestFit(void) { double XYs, Xsquares, meanX, meanY; double sq(double); int i; XYs = Xsquares = meanX = meanY = 0.0; for (i = 0; i < rows; ++i) { meanX += log_r[i]; meanY += log_Z[i]; } meanX /= rows; meanY /= rows; for (i = 0; i < rows; ++i) { XYs += (log_r[i] - meanX) * (log_Z[i] - meanY); Xsquares += sq(log_r[i] - meanX); } slope = XYs / Xsquares; intercept = meanY - slope * meanX; } double smoothed(int i) { return(exp(intercept + slope * log(i))); } int row(int i) { int j = 0; while (j < rows && r[j] < i) ++j; return((j < rows && r[j] == i) ? j : -1); } void showEstimates(void) { int i; printf("0\t%.4g\n", PZero); for (i = 0; i < rows; ++i) printf("%d\t%.4g\n", r[i], p[i]); } void analyseInput(void) { int i, j, next_n; double k, x, y; int indiffValsSeen = FALSE; int row(int); void findBestFit(void); double smoothed(int); double sq(double); void showEstimates(void); bigN = 0; for (j = 0; j < rows; ++j) bigN += r[j] * n[j]; next_n = row(1); PZero = (next_n < 0) ? 0 : n[next_n] / (double) bigN; for (j = 0; j < rows; ++j) { i = (j == 0 ? 0 : r[j - 1]); if (j == rows - 1) k = (double) (2 * r[j] - i); else k = (double) r[j + 1]; Z[j] = 2 * n[j] / (k - i); log_r[j] = log(r[j]); log_Z[j] = log(Z[j]); } findBestFit(); for (j = 0; j < rows; ++j) { y = (r[j] + 1) * smoothed(r[j] + 1) / smoothed(r[j]); if (row(r[j] + 1) < 0) indiffValsSeen = TRUE; if (! indiffValsSeen) { x = (r[j] + 1) * (next_n = n[row(r[j] + 1)]) / (double) n[j]; if (fabs(x - y) <= CONFID_FACTOR * sqrt(sq(r[j] + 1.0) * next_n / (sq((double) n[j])) * (1 + next_n / (double) n[j]))) indiffValsSeen = TRUE; else rStar[j] = x; } if (indiffValsSeen) rStar[j] = y; } bigNprime = 0.0; for (j = 0; j < rows; ++j) bigNprime += n[j] * rStar[j]; for (j = 0; j < rows; ++j) p[j] = (1 - PZero) * rStar[j] / bigNprime; showEstimates(); }