tstmid.c

雷达工具箱

    if(13 != iwhich) goto S290;
/*
     Multinomial outcomes
*/
    puts(" Enter (int) number of observations: ");
    scanf("%ld",&ntry);
    printf(" Enter (int) num. of categories: <= %ld\n",MXNCAT);
    scanf("%ld",&ncat);
    while(ncat > MXNCAT) {
      printf(" number of categories must be <= %ld\n",MXNCAT);
      puts(" Try again ... ");
      printf(" Enter (int) num. of categories: <= %ld\n",MXNCAT);
      scanf("%ld",&ncat);
    }
    printf(" Enter (real) prob. vector of length %ld\n",ncat-1L);
    for(i=1L;i<ncat;i++) scanf("%f",prob + i - 1);
    genmul(ntry,prob,ncat,iarray);
    ntot = 0L;
    if (ntry > 0L) {
      rtry = (float) ntry;
      for(i=0L;i<ncat;i++) {
	ntot += *(iarray + i);
	*(array + i) = *(iarray + i) / rtry;
      }
    } else {
      for(i=0L;i<ncat;i++) {
	ntot += *(iarray + i);
	*(array + i) = 0.0;
      }
    }
    psum = 0.0;
    for(i=1L;i<ncat;i++) psum += *(prob + i - 1);
    *(prob + ncat - 1) = 1.0 - psum;
    printf(" Total number of observations: %ld\n",ntot);
    puts(" Total observations by category: ");
    pr_ln_ar(iarray,ncat);
    puts(" True probabilities by category: ");
    pr_fl_ar(prob,ncat);
    puts(" Observed proportions by category: ");
    pr_fl_ar(array,ncat);
    goto S320;

S290:
    if(14 != iwhich) goto S300;
/*
     Normal deviates
*/
    strcpy( pr_type, "norm");
    puts(" Enter (real) AV, (real) SD for Normal");
    scanf("%f%f",param,(param+1));
    for(i=1; i<=1000; i++) *(array+i-1) = gennor(*param, *(param+1));
    statistics(array,1000L,&av,&var,&xmin,&xmax);
    trstatistics(pr_type,param,&avtr,&vartr);
    pr_ans( av, avtr, var, vartr, xmin, xmax );
    goto S320;

S300:
    if(15 != iwhich) goto S310;
/*
     Negative Binomial outcomes
*/
    strcpy( pr_type, "nbin");
    puts(" Enter required (int) Number of events, then");
    puts(" (real) Prob event for negative binomial outcomes");
    scanf("%ld%f",&ntry,&pevt);
    for(i=1; i<=1000; i++) *(iarray+i-1) = ignnbn(ntry,pevt);
    for(i=1; i<=1000; i++) *(array+i-1) = *(iarray+i-1);
    statistics(array,1000L,&av,&var,&xmin,&xmax);
    *param = ntry;
    *(param+1) = pevt;
    trstatistics(pr_type,param,&avtr,&vartr);
    pr_ans( av, avtr, var, vartr, xmin, xmax );
    goto S320;

S310:
S320:
    goto S10;
#undef mxwh
#undef MXNCAT
}

void trstatistics(char pr_type[4],float *parin,float *av,float *var)
/*
**********************************************************************
     SUBROUTINE TRSTATISTICS( TYPE, PARIN, AV, VAR )
               TRue STATistics
     Returns mean and variance for a number of statistical distribution
     as a function of their parameters.
                              Arguments
     TYPE --> Character string indicating type of distribution
             'chis' chisquare
             'ncch' noncentral chisquare
             'f   ' F (variance ratio)
             'ncf ' noncentral f
             'unif' uniform
             'beta' beta distribution
             'bin ' binomial
             'pois' poisson
             'exp ' exponential
             'gamm' gamma
             'norm' normal
             'nbin' negative binomial
                         CHARACTER*(4) TYPE
     PARIN --> Array containing parameters of distribution
              chisquare
               PARIN(1) is df
              noncentral chisquare
               PARIN(1) is df
               PARIN(2) is noncentrality parameter
              F (variance ratio)
               PARIN(1) is df numerator
               PARIN(2) is df denominator
              noncentral F
               PARIN(1) is df numerator
               PARIN(2) is df denominator
               PARIN(3) is noncentrality parameter
              uniform
               PARIN(1) is LOW bound
               PARIN(2) is HIGH bound
              beta
               PARIN(1) is A
               PARIN(2) is B
              binomial
               PARIN(1) is Number of trials
               PARIN(2) is Prob Event at Each Trial
              poisson
               PARIN(1) is Mean
              exponential
               PARIN(1) is Mean
              gamma
               PARIN(1) is A
               PARIN(2) is R
              normal
               PARIN(1) is Mean
               PARIN(2) is Standard Deviation
              negative binomial
               PARIN(1) is required Number of events
               PARIN(2) is Probability of event
                         REAL PARIN(*)

     AV <-- Mean of specified distribution with specified parameters
                         REAL AV
     VAR <-- Variance of specified distribution with specified paramete
                         REAL VAR
                              Note
     AV and Var will be returned -1 if mean or variance is infinite
*/
{
static float a,b,range;

    if(strcmp("chis",pr_type) !=0 ) goto S10;
    *av = *parin;
    *var = 2.0**parin;
    goto S210;
S10:
    if(strcmp("ncch",pr_type) != 0) goto S20;
    a = *parin+*(parin+1);
    b = *(parin+1)/a;
    *av = a;
    *var = 2.0*a*(1.0+b);
    goto S210;
S20:
    if(strcmp("f   ",pr_type) != 0) goto S70;
    if(!(*(parin+1) <= 2.0001)) goto S30;
    *av = -1.0;
    goto S40;
S30:
    *av = *(parin+1)/(*(parin+1)-2.0);
S40:
    if(!(*(parin+1) <= 4.0001)) goto S50;
    *var = -1.0;
    goto S60;
S50:
    *var = 2.0*pow(*(parin+1),2.0)*(*parin+*(parin+1)-2.0)/(*parin*pow(*(parin
      +1)-2.0,2.0)*(*(parin+1)-4.0));
S60:
    goto S210;
S70:
    if(strcmp("ncf ",pr_type) != 0) goto S120;
    if(!(*(parin+1) <= 2.0001)) goto S80;
    *av = -1.0;
    goto S90;
S80:
    *av = *(parin+1)*(*parin+*(parin+2))/((*(parin+1)-2.0)**parin);
S90:
    if(!(*(parin+1) <= 4.0001)) goto S100;
    *var = -1.0;
    goto S110;
S100:
    a = pow(*parin+*(parin+2),2.0)+(*parin+2.0**(parin+2))*(*(parin+1)-2.0);
    b = pow(*(parin+1)-2.0,2.0)*(*(parin+1)-4.0);
    *var = 2.0*pow(*(parin+1)/ *parin,2.0)*(a/b);
S110:
    goto S210;
S120:
    if(strcmp("unif",pr_type) != 0) goto S130;
    range = *(parin+1)-*parin;
    *av = *parin+range/2.0;
    *var = pow(range,2.0)/12.0;
    goto S210;
S130:
    if(strcmp("beta",pr_type) != 0) goto S140;
    *av = *parin/(*parin+*(parin+1));
    *var = *av**(parin+1)/((*parin+*(parin+1))*(*parin+*(parin+1)+1.0));
    goto S210;
S140:
    if(strcmp("bin ",pr_type) != 0) goto S150;
    *av = *parin**(parin+1);
    *var = *av*(1.0-*(parin+1));
    goto S210;
S150:
    if(strcmp("pois",pr_type) != 0) goto S160;
    *av = *var = *parin;
    goto S210;
S160:
    if(strcmp("exp ",pr_type) != 0) goto S170;
    *av = *parin;
    *var = *av * *av;
    goto S210;
S170:
    if(strcmp("gamm",pr_type) != 0) goto S180;
    *av = *(parin + 1) / *parin;
    *var = *av / *parin;
    goto S210;
S180:
    if(strcmp("norm",pr_type) != 0) goto S190;
    *av = *parin;
    *var = *(parin + 1) * *(parin + 1);
    goto S210;
S190:
    if(strcmp("nbin",pr_type) != 0) goto S200;
    *av = *parin * (1.0 - *(parin + 1)) / *(parin + 1);
    *var = *av / *(parin + 1);
    goto S210;
S200:
    puts("Unimplemented pr_type in TRSTATISTICS");
    exit(1);
S210:
    return;
}

void pr_ans( float av, float avtr, float var, float vartr,
             float xmin, float xmax) 
{

    printf(" Mean Generated          %15.7g     True          %15.7g\n",
           (double)av, (double)avtr);
    printf(" Variance Gener          %15.7g     True          %15.7g\n",
           (double)var, (double)vartr);
    printf(" Minimum                 %15.7g     Maximum       %15.7g\n",
           (double)xmin, (double)xmax);
}

void pr_ln_ar(long *array, long larray)
/*
  Print array of longs; length of array is larray.  Print 10 to a
  line using 7 characters for each field.
*/
{
    int lo, hi, i, up;
    up = larray-1;
    for (lo=0; lo <=up; lo +=10) {
        hi = min(lo+9,up);
    	for (i=lo; i<=hi; i++)
            printf(" %7ld",*(array+i));
    puts("");    
    }
}    

void pr_fl_ar(float *array, long larray)
/*
  Print array of floats; length of array is larray.  Print 8 to a
  line using 9 characters for each field.
*/
{
    int lo, hi, i, up;
    up = larray-1;
    for (lo=0; lo <=up; lo +=8) {
        hi = min(lo+7,up);
    	for (i=lo; i<=hi; i++)
            printf(" %.7f",*(array+i));
    puts("");    
    }
}