A3.2. Analysis of continuous travel time data

This first section contains the code to undertake the analysis of travel time data performed in chapter 4. The code is contained in 5 files, and a data file must be supplied. Before providing the code, I will briefly discuss the structure of the files and the routines contained within. I will also provide a sample data file and the output generated from the code with the sample data file.

main.c - contains the routine main() that controls the program by calling other routines. This routine also reads in the data file.

mle.c - contains the routines r_mle() and sig_mle() that compute maximum likelihood estimates of r and .

conf_int.c - contains the routines r_mle() and sig_mle() that determine (1-) · 100 percent confident intervals for the parameters and prints them out.

cumulative.c - contains the routine cumulative() that converts the travel time data to cumulative travel time based on the model and estimates parameters.

chi.c - contains the routine chi() that performs the goodness-of-fit test and prints out the results. This routine must be provided with a vector of cumulative travel times.

main.c

#include <stdio.h>

float   r_mle();
float   sig_mle();
void    confidence_intervals();
float   cumulative();
void    chi();

void    main()
{
 float  length;         /* length of the reach */
 int            num;    /* number of individuals */
 char   junk[12];       /* place to put header info from data file */
 float  *tt_vec;        /* vector to store travel time data */
 float  *cum_tt_vec; /* vector for values from tt cdf */
 int            i;              /* increment for tt vector */
 float  r, sig;         /* model parameters */
 float  alpha;          /* 1 - alpha is width of conf. intervals */
 int            pars;           /* # of parameters used in the model */
 FILE   *data;          /* pointer to data file */

 data = fopen("tt.data", "r");

 fscanf(data, "%s%s%f", junk, junk, &length);
 fscanf(data, "%s%s%d", junk, junk, &num);
 fscanf(data, "%s%s", junk, junk);

 /* allocate memory for the travel times vectors */
 tt_vec = (float *) malloc(num*sizeof(float));
 cum_tt_vec = (float *) malloc(num*sizeof(float));
 
 /* read in travel times from data file */
 for(i = 0; i < num; ++i)
        fscanf(data, "%f", &tt_vec[i]);

 /* compute maximum likelihood estimates */
 r = r_mle(tt_vec, length, num);
 sig = sig_mle(tt_vec, length, num);
 printf("mle r = %6.3f/n", r);
 printf("mle sig = %6.3f/n", sig);

 /* 95% confidence interval */  
 alpha = 0.05; 
 confidence_intervals(r, sig, length, num, alpha);
 
 /* X-squared goodness-of-fit test */
 /* the test needs values from the cumulative dist. func. */
 for(i = 0; i < num; ++i)
        cum_tt_vec[i] = cumulative(r, sig, length, tt_vec[i]);
 
 pars = 2; /* number of parameters used by the model */
 chi(cum_tt_vec, pars, num);

 fclose(data);
}

#include <math.h>

/* computes maximum likelihood estimate for the parameter r */
/* based on the travel time data                         */
float           r_mle(tt_vec, pool_length, num)
        float           *tt_vec;                        /* vector of travel times for group */
        float           pool_length;
        int             num;                    /* number of individuals in group */
{
        float           tt_bar = 0;      /* average travel time */
        int             i;

        for (i = 0; i < num; ++i){
                tt_bar += tt_vec[i];
 }
 
        tt_bar = tt_bar/num;

        return(pool_length/tt_bar);
}

/* computes maximum likelihood estimate for the parameter sigma */
/* based on the travel time data                         */
float           sig_mle(tt_vec, pool_length, num)
        float           *tt_vec;
        float           pool_length;
        int             num;
{
        float           tt_bar = 0;                             /* arithmetic mean */
        float           tt_recip = 0;                   /* harmonic mean */
        int              i;
 
        for (i = 0; i < num; ++i){
                tt_bar += tt_vec[i];
                tt_recip += 1.0/tt_vec[i];
        }
        tt_bar = tt_bar/num;
        tt_recip = tt_recip/num;

        return(pool_length * sqrt(tt_recip - (1.0/tt_bar)));
}

#include <math.h>
#include <stdio.h>

/* This routine is passed the maximum likelihood estimates for r */
/* and sigma, reach length, number of fish and alpha. It prints */
/* 100*(1-alpha) percent confidence intervals for the parameters */
/* r and sigma. The appropriate quantiles of the Student's t and */
/* chi-square are obtained from S-plus, which is provided with */
/* the degrees of freedom.                               */

void            confidence_intervals(r, s, L, num, alpha)
        float           r, s; /* mles of r and sigma */
        float           L; /* reach length */
        int             num; /* number of individuals */
        float           alpha; /* 1-alpha is length of C.I. */
{
        float           r_min, r_max; /* min and max of r C.I. */
        float           sig_min, sig_max; /* min and max of sigma C.I. */
        float           a, b; /* quantiles used in C.I. calc. */
        char            junk[10];        /* junk from input file */
        FILE            *iptr, *optr;    /* input and output files */

        /* provide degrees of freedom and alpha for S-plus routine */
        optr = fopen(".quant_info", "w");
        fprintf(optr, "%d/t%f/n", num - 1, alpha);
        fclose(optr);

        /* execute S-plus routine that prints quantiles to a file */
        system("S < quantile.s > /dev/null");

        /* open file and read in (1.0-alpha/2)th quantile of t dist. */
        iptr = fopen(".quantile", "r");
        fscanf(iptr,"%s%f", junk, &a);

        /* compute max and min values of r C.I. */
        r_min = r*(1.0 - a*sqrt((s*s)/(r*L*(num-1))));
        r_max = r*(1.0 + a*sqrt((s*s)/(r*L*(num-1))));

        /* print r C.I. */
        printf("%4.1f percent confidence interval for r:/n",
                100*(1.0-alpha));
        printf("(%6.2f,%6.2f)/n", r_min, r_max);

        /* read (alpha/2)th and (1-alpha/2)th quantiles of chi-sq. dist. */
        fscanf(iptr,"%s%f", junk, &a);
        fscanf(iptr,"%s%f", junk, &b);

        /* compute max and min values of sig C.I. */
        sig_min = s*sqrt((float)(num)/a);
        sig_max = s*sqrt((float)(num)/b);

        /* print sigma C.I. */
        printf("%4.1f percent confidence interval for sigma:/n",
                100*(1.0-alpha));
        printf("(%6.2f,%6.2f)/n", sig_min, sig_max);
}

#include <math.h>
#include "input.h"

#define pi      3.1415

/* phi is the cumulative distribution for a standard normal */
float phi(x)
float x;
{
        return(0.5 + erf(x/sqrt(2.0))/2);
}

/* this routine returns a value from the cumulative distribution */
/* of the basic travel time model. The routine must be passed /x11*/
/* the model parameters and the travel time. The procedure for */
/* generating the value is described in appendix 4.a /x11/x11/x11/x11/x11/x11/x11/x11/x11/x11*/
float   cumulative(r, sig, L, t)
        float                   r, sig, L; /* model parameters */
        float                   t; /* travel time */    
{
        float                   mu, lam;                /* reparameterization */
        float                   first, second;
        double                  y, z;
        double                  z2, z4;
        double                  fac1, fac2;
        float                   d0 = 0.2316419;
        float                   d1 = 0.319381530;
        float                   d2 = -0.356563782;
        float                   d3 = 1.781477937;
        float                   d4 = -1.821255978;
        float                   d5 = 1.330274429;
        float                   qz, qz2, qz4;

        mu = L/r;
        lam = L/sig;

        y = lam*(t-mu)/(mu*sqrt(t));
        z = lam*(t+mu)/(mu*sqrt(t));
 
        if (z<4){
                qz = 1/(1+d0*z);
                qz2 = qz*qz;
                qz4 = qz2*qz2;
 
                fac1 = (exp(-(y*y)/2))/(sqrt(2*pi));
                fac2 = (d1*qz + d2*qz2 + d3*qz2*qz
                        + d4*qz4 + d5*qz4*qz);
        }
        if (z>=4){
                z2 = z*z;
                z4 = z2*z2;
 
                fac1 = (exp(-(y*y)/2))/(sqrt(2*pi))/z;
                fac2 = 1 - 1/z2 + 3/z4 - 3*5/(z2*z4) + 3*5*7/(z4*z4)
                        - 3*5*7*9/(z4*z4*z2) + 3*5*7*9*11/(z4*z4*z4)
                        - 3*5*7*9*11*13/(z4*z4*z4*z2);
        }
        second = fac1*fac2;
        first = phi(y);

        return(first + second);
}

#include <math.h>
#include <stdio.h>

float           gammp(); /* This routine is from Numerical Recipes in C */
                                /* Press, et al. 1988 */

void            chi(cum_vec, params, num_fish)
        float           *cum_vec;                       /* vector cdf values of tt dist. */
        int             params;                 /* number of parameters estimated */
        int             num_fish;                       /* number of individuals */
{
        int             num_bins;                       /* number of bins */
        float           bin_width;                      /* width of each bin */
        float           expect;                 /* expected individuals per bin */
        int             *obs;                   /* vector of oberved individuals */
        int             i;                      /* counter */
        float           X = 0.0;                        /* chi square statistic */
        float           prob;                   /* chi square probability */
        int             df;                     /* degrees of freedom */
 
         /* num_bins is determined by Mann-Wald calculation */
        num_bins = (int)(3.76*pow((float)(num_fish), 0.4) );

        bin_width = 1.0/num_bins;
        expect = bin_width*num_fish;

        /* allocate memmory for vector of observed values and set */
        /* each element to zero */
        obs = (int *) malloc(num_bins*sizeof(int));
        bzero((char *) obs, num_bins*sizeof(float));

        /* determine which bin each individual falls into */
        for ( i = 0; i < num_fish; ++i)
                ++obs[(int)(cum_vec[i]/bin_width)];

        /* compute chi square statistic */
        for (i = 0; i < num_bins; ++i)
                X += ((expect-obs[i])*(expect-obs[i]))/expect;

        df = num_bins - params - 1;
        
        /* compute percentile of chi-square distribution */
        prob = gammp((float)(df)/2.0, X/2.0);
 
        printf("/nX-squared goodness-of-fit test/n");
        printf("X-squared = %7.3f/n", X);
        printf("degrees of freedom = %3d/n", df);
        printf("p = %7.3f/n", 1.0 - prob);
}


sample data file

reach length: 52.0
num fish: 57
travel times:
15.74    12.06   30.63   24.79   17.54   
26.39    14.87   9.25    4.83    12.32   
14.61    9.08    20.34   7.74    16.98   
3.99     10.69   23.38   20.02   19.74   
22.66    24.62   20.62   18.24   22.48   
10.76    12.01   9.99    6.34    21.47   
18.09    22.25   15.74   13.68   5.11    
10.35    10.41   22.41   8.21    36.66   
21.45    13.17   18.64   18.69   11.85   
20.57    34.52   15.73   9.46    37.39   
21.53    92.03   33.30   21.67   21.94   
21.45    8.23   

        program output

mle r = 2.773
mle sig = 7.251

95.0 percent confidence interval for r:
( 2.33, 3.22)
95.0 percent confidence interval for sigma:
( 6.18, 8.97)

X-squared goodness-of-fit test
X-squared = 22.263
degrees of freedom = 15
p = 0.101