home *** CD-ROM | disk | FTP | other *** search
/ Explore the World of Soft…e: Engineering & Science / Explore_the_World_of_Software_Engineering_and_Science_HRS_Software_1998.iso / programs / statistc / aiuq9308.exe / AIU-QLV.CXX next >
Text File  |  1993-08-01  |  14KB  |  412 lines
  1. /*
  2.  *    Acceptable Quality Level (AQL), Indifference Quality Level (IQL) and
  3.  *    Unacceptable Quality Level (UQL) for Single-Stage Binomial Sampling
  4.  *
  5.  *                   A Quality Assurance Training Tool:
  6.  *            Statistics Committee of the QA Section of the PMA
  7.  *
  8.  *               Bob Obenchain, CompuServe User [72007,467]
  9.  *
  10.  *                     Usage: DOSprompt>  aiu-qlv
  11.  *
  12.  *      Parameter settings are either...
  13.  *             specified interactively in response to prompts, or
  14.  *             redirected to come from a batch input file.
  15.  */
  16.  
  17. #include <math.h>
  18. #include <stdio.h>
  19. #include <signal.h>
  20. #include <stdlib.h>
  21. #include <string.h>
  22. #include <process.h>
  23.  
  24. #define abs(a)      (( (a) >= 0 ) ? (a) : (-a))  /* absolute value of a */
  25. #define sgn(a)      (( (a) >= 0 ) ?  1  :  -1 )  /* numerical sign of a */
  26. #define max(x,y)    (( (x) > (y) ) ?  (x) : (y) )
  27. #define min(x,y)    (( (x) > (y) ) ?  (y) : (x) )
  28. #define ITMAX   100
  29. #define EPS     3.0e-7
  30.  
  31. /* externals */
  32.  
  33. FILE    *inpfile, *outfile;
  34. char    buf[256], inpnam[13], outnam[13], dtlnam[13], yesno;
  35. double  alpha, omalpha, term, quan, facp1, fnmac;
  36. long    sn, ac, rj;
  37.  
  38. /* prototypes */
  39.  
  40. void   main( void );
  41. double qbeta( double p, double a, double b );
  42. double qnorms( double pr );
  43. double betai( double a, double b, double x );
  44. double betacf( double a, double b, double x );
  45. double gammln( double xx );
  46. extern int get1key( void );
  47.  
  48. double qbeta(p,a,b)
  49. double p, a, b; {
  50.       int right, i;
  51.       double quan, tol, xl, pl, xr, pr, xm, sd, xx, pp, pdif;
  52.       tol = 1.17577E-37;
  53.  
  54.       quan = p;
  55.       if( a == 1.0 && b == 1.0 )
  56.           return( quan );
  57.  
  58.       if( p < 0.05 ) {
  59.           xl = 0.0;
  60.           pl = 0.0;
  61.           right = 0;
  62.           }
  63.       else if( p > 0.95 ) {
  64.           xr = 1.0;
  65.           pr = 1.0;
  66.           right = 1;
  67.           }
  68.       else {
  69.           /* normal approximations... */
  70.           xm = a + b;
  71.           sd = sqrt((a*b)/(xm*xm*(xm+1.)));
  72.           xl = qnorms(p) * sd + a / xm;
  73.           if( xl < tol ) {
  74.               xl=0.0;
  75.               pl=0.0;
  76.               right=0;
  77.               }
  78.           else if( xl > (1.0-tol) ) {
  79.               xr=1.0;
  80.               pr=1.0;
  81.               right=1;
  82.               }
  83.           else {
  84.               pl=betai(a,b,xl);
  85.               right=0;
  86.               if( pl > p ) {
  87.                   right=1;
  88.                   pr=pl;
  89.                   xr=xl;
  90.                   }
  91.               }
  92.           }
  93.  
  94.       if( right == 1 ) {
  95.  
  96.           /* find left end... */
  97.           while( 1 ) {
  98.               xl = max( xr-0.05, 0.0 );
  99.               if( xl <= 0.0 ) {
  100.                   pl = 0.0;
  101.                   break;
  102.                   }
  103.               pl=betai(a,b,xl);
  104.               quan = xl;
  105.               if( pl == p )
  106.                   return( quan );
  107.               if( pl < p )
  108.                   break;
  109.               xr=xl;
  110.               pr=pl;
  111.               }
  112.           }
  113.  
  114.       if( right == 0 ) {
  115.  
  116.           /* find right end... */
  117.           while( 1 ) {
  118.               xr = min( xl+0.05, 1.0 );
  119.               if( xr >= 1.0 ) {
  120.                   pr = 1.0;
  121.                   break;
  122.                   }
  123.               pr=betai(a,b,xr);
  124.               quan = xr;
  125.               if( pr == p )
  126.                   return( quan );
  127.               if( pr > p )
  128.                   break;
  129.               xl=xr;
  130.               pl=pr;
  131.               }
  132.           }
  133.  
  134.       /* now quantile is between xl and xr; use bisection */
  135.       for( i=1; i<=6; i++ ) {
  136.           xx = (xl+xr)*.5;
  137.           pp = betai(a,b,xx);
  138.           pdif = pp-p;
  139.           quan = xx;
  140.           if( abs(pdif) < tol )
  141.               return( quan );
  142.           if( pdif > 0.0 ) {
  143.               xr=xx;
  144.               pr=pp;
  145.               }
  146.           else {
  147.               xl = xx;
  148.               pl = pp;
  149.               }
  150.           }
  151.  
  152.       /* next, try secant method... */
  153.       for( i=1; i<=6; i++ ) {
  154.           xx = xl+(p-pl)*(xr-xl)/(pr-pl);
  155.           pp = betai(a,b,xx);
  156.           pdif = pp-p;
  157.           quan =xx;
  158.           if( abs(pdif) < tol)
  159.               return( quan );
  160.           if( pdif > 0.0) {
  161.               xr = xx;
  162.               pr = pp;
  163.               }
  164.           else {
  165.               xl = xx;
  166.               pl = pp;
  167.               }
  168.           }
  169.       quan = xx;
  170.       /* failed to converge... */
  171.       return( quan );
  172.       }
  173.  
  174. double qnorms( pr )
  175. double pr; {
  176.       double quan, p, eta, f1, f2, f3, f4, f5, f6;
  177.       quan = 0.0;
  178.       if( pr <= 0.0 || pr == 0.5 || pr >= 1.0 )
  179.           return( quan );
  180.       f1 = .010328;
  181.       f2 = .802853;
  182.       f3 = 2.515517;
  183.       f4 = .001308;
  184.       f5 = .189269;
  185.       f6 = 1.432788;
  186.       p = pr;
  187.       if( pr > 0.5 )
  188.           p = 1.0 - pr;
  189.       eta = sqrt( -2.0 * log( p ) );
  190.       quan = eta - ((f1*eta+f2)*eta+f3)/(((f4*eta+f5)*eta+f6)*eta+1.0);
  191.       if( pr <= 0.5 )
  192.           quan *= -1.0;
  193.       return( quan );
  194.       }
  195.  
  196. void main() {
  197.  
  198.       char kb;
  199.  
  200.       printf("\n\n                 Acceptable Quality Level (AQL)," );
  201.       printf(  "\n               Indifference Quality Level (IQL) and" );
  202.       printf(  "\n               Unacceptable Quality Level (UQL) for" );
  203.       printf(  "\n                 Single-Stage Binomial Sampling" );
  204.       printf("\n\n                 A/I/U-QLV.EXE....Version 9308");
  205.       printf("\n\n                A Quality Assurance Training Tool:");
  206.       printf(  "\n         Statistics Committee of the QA Section of the PMA");
  207.       printf("\n\n            Bob Obenchain, CompuServe User [72007,467]\n\n" );
  208.       
  209.       printf(      "\tWill Parameter Settings be Input via K = Keyboard ? ...or" );
  210.       printf(    "\n\t                                     B = Batch File ?" );
  211.       printf(  "\n\n\tPress the   K   or   B   key now --> " );
  212.       while( 1 ) {
  213.               gets( buf );
  214.               sscanf( buf, "%c", &kb );
  215.               if( kb == 'k' )
  216.                   kb  = 'K';
  217.               if( kb == 'b' )
  218.                   kb  = 'B';
  219.               if( kb == 'K' || kb == 'B' )
  220.                   break;
  221.               printf( "\r\tENTER either K  or  B  ...Try Again --> " );
  222.               }
  223.  
  224.       strcpy( outnam, "aiu-qlv" );
  225.       strcpy( inpnam, "aiu-qlv" );
  226.       strcpy( dtlnam, "aiu-qlv" );
  227.  
  228.       if( kb == 'B' )
  229.               printf( "\n\n\tBatch Input of Parameter Settings Selected..." );
  230.       else
  231.               printf( "\n\n\tKeyboard Input of Parameter Settings Selected..." );
  232.       
  233.       printf( "\n\n\tAt colon Prompts : ...simply press ENTER to get the [default]." );
  234.  
  235.       if( kb == 'B' ) {
  236.               printf(  "\n\n\tSpecify filename of Batch Input File [%s] : ",
  237.                   inpnam );
  238.               gets( buf );
  239.               if( buf[0] != '\0' && buf[0] != '\r' && buf[0] != '\n' )
  240.                   sscanf( buf, "%s", inpnam );
  241.               inpnam[ min( 8, (int)strcspn( inpnam, ". " ) ) ] = '\0';
  242.               if( inpnam[0] == '\0' )
  243.                   strcpy( inpnam, "aiu-qlv" );
  244.               strcpy( outnam, inpnam );
  245.               strcpy( dtlnam, inpnam );
  246.               strcat( inpnam, ".inp" );
  247.               printf( "\n\tThe Batch Input file is to be: %s\n", inpnam );
  248.               if( ( inpfile = fopen( inpnam, "r" ) ) == NULL ) {
  249.                       printf(  "\tCannot read Batch Input file: %s\n", inpnam );
  250.                       printf(  "\t...using Keyboard Input from standard Infile, stdin.\n" );
  251.                       kb = 'K';
  252.                       }
  253.               }
  254.       
  255.       if( kb == 'K' )
  256.               inpfile = stdin;
  257.  
  258.       printf(  "\n\n\tSpecify filename to save primary A/I/U-QLV Output [%s] : ",
  259.           outnam );
  260.       fgets( buf, 80, inpfile );
  261.       if( buf[0] != '\0' && buf[0] != '\r' && buf[0] != '\n' )
  262.           sscanf( buf, "%s", outnam );
  263.       outnam[ min( 8, (int)strcspn( outnam, ". " ) ) ] = '\0';
  264.       if( outnam[0] == '\0' )
  265.           strcpy( outnam, "aiu-qlv" );
  266.       strcpy( dtlnam, outnam );
  267.       strcat( outnam, ".out" );
  268.       printf(  "\n\tThe A/I/U-QLV primary Output file is to be: %s\n",
  269.               outnam );
  270.       if( ( outfile = fopen( outnam, "w" ) ) == NULL ) {
  271.               printf(  "\tCannot write to Output filename : %s\n", outnam );
  272.               if( ( outfile = fopen( "aiu-qlv.out", "w" ) ) == NULL )
  273.                       printf(  "\tCannot write to Output filename : aiu-qlv.out\n" );
  274.               else {
  275.                       printf(  "\t...using default Outfile name : aiu-qlv.out\n" );
  276.                       strcpy( outnam, "aiu-qlv.out" );
  277.                       }
  278.               }
  279.  
  280.       if( outfile != NULL )
  281.           fprintf( outfile, "Output file for A/I/U-QLV..." );
  282.      
  283.       printf("\n\nThis software calculates lower and upper confidence limits on the" );
  284.       printf("\npercentage of nonconforming units in an infinite population of");
  285.       printf("\nunits given accept - reject numbers for random sampling.");
  286.       
  287.       sn = 100;
  288.       ac =   0;
  289.       rj =   1;
  290.       omalpha = 0.95;
  291.  
  292.       while( 1 ) {
  293.  
  294.           printf("\n\nWhat is the Sample Size ? [n=%ld] : ", sn );
  295.           fgets( buf, 80, inpfile );
  296.           if( buf[0] != '\0' && buf[0] != '\r' && buf[0] != '\n' )
  297.               sscanf( buf, "%ld", &sn );
  298.           if( sn < (long)1 )
  299.               sn = (long)1;
  300.           printf("\n\tSample Size: n =%7ld", sn);
  301.           if( outfile != NULL ) {
  302.               fprintf( outfile, "\n\n*****************\n" );
  303.               fprintf( outfile, "\nSample Size: n =%7ld", sn);
  304.               }
  305.  
  306.           printf("\n\nWhat is the Accept Number? [ac=%ld] : ",
  307.               ac );
  308.           fgets( buf, 80, inpfile );
  309.           if( buf[0] != '\0' && buf[0] != '\r' && buf[0] != '\n' )
  310.               sscanf( buf, "%ld", &ac );
  311.           if( ac > sn - (long)1 )
  312.               ac = sn - (long)1 ;
  313.           if( ac < (long)0 )
  314.               ac = (long)0;
  315.           printf("\n\tAccept Number: ac =%7ld", ac );
  316.           if( outfile != NULL )
  317.               fprintf( outfile, "\nAccept Number: ac =%7ld", ac );
  318.  
  319.           rj = ac + (long)1;
  320.           printf("\n\tReject Number: rj =%7ld", rj );
  321.           if( outfile != NULL )
  322.               fprintf( outfile, "\nReject Number: rj =%7ld", rj );
  323.  
  324.           printf("\n\nDesired   Yield   Fraction at AQL" );
  325.           printf(  "\n    and Rejection Fraction at UQL ? [1-alpha=%5.3lf] : ",
  326.               omalpha );
  327.           fgets( buf, 80, inpfile );
  328.           if( buf[0] != '\0' && buf[0] != '\r' && buf[0] != '\n' )
  329.               sscanf( buf, "%lf", &omalpha );
  330.           if( omalpha > 0.99 )
  331.               omalpha = 0.99;
  332.           if( omalpha < 0.75 ) {
  333.                  printf("\n\tWARNING: Yield/Rejection Fraction is usually large.");
  334.                  omalpha = 0.75;
  335.                  }
  336.           printf("\n\tYield/Rejection Fraction = %5.3lf", omalpha);
  337.           if( outfile != NULL )
  338.               fprintf( outfile, "\nYield/Rejection Fraction = %5.3lf", omalpha);
  339.  
  340.           alpha = 1.0 - omalpha;
  341.  
  342.           printf("\n\n\nAcceptable Percent Nonconforming: 100 * %ld / %ld = %7.3lf%%",
  343.               ac, sn, 100.0 * (double)ac / (double)sn );
  344.           printf("\n\nRejectable Percent Nonconforming: 100 * %ld / %ld = %7.3lf%%",
  345.               rj, sn, 100.0 * (double)rj / (double)sn );
  346.           if( outfile != NULL ) {
  347.               fprintf( outfile, "\nAcceptable Percent Nonconforming: 100 * %ld / %ld = %7.3lf%%",
  348.                   ac, sn, 100.0 * (double)ac / (double)sn );
  349.               fprintf( outfile, "\nRejectable Percent Nonconforming: 100 * %ld / %ld = %7.3lf%%",
  350.                   rj, sn, 100.0 * (double)rj / (double)sn );
  351.               }
  352.  
  353.           printf("\n\nAt Yield/Rejection Percentage %5.1lf%%,", 100.0 * omalpha );
  354.  
  355.           if( rj == sn )
  356.               quan = 100.0 * exp( log(alpha) / (double)sn );
  357.           else if( rj == (long)0 )
  358.               quan=0.0;
  359.           else {
  360.               facp1 = (double)(sn-rj+(long)1);
  361.               fnmac = (double)rj;
  362.               quan = 100.0 * (1.0 - qbeta( omalpha, facp1, fnmac ));
  363.               }
  364.  
  365.           printf(    "\n   ...the AQL is %7.3lf%%", quan );
  366.           if( outfile != NULL ) {
  367.               fprintf( outfile, "\n\nAt Yield/Rejection Percentage %5.1lf%%,", 100.0 * omalpha );
  368.               fprintf( outfile, "\n   ...the AQL is %7.3lf%%", quan );
  369.               }
  370.  
  371.           if( rj == sn )
  372.               quan = 100.0 * exp( log(0.5) / (double)sn );
  373.           else if( rj == (long)0 )
  374.               quan=0.0;
  375.           else {
  376.               facp1 = (double)(sn-rj+(long)1);
  377.               fnmac = (double)rj;
  378.               quan = 100.0 * (1.0 - qbeta( 0.5, facp1, fnmac ));
  379.               }
  380.  
  381.           printf(    "\n   ...the IQL is %7.3lf%%", quan );
  382.           if( outfile != NULL )
  383.               fprintf( outfile, "\n   ...the IQL is %7.3lf%%", quan );
  384.  
  385.           if( ac == (long)0 )
  386.               quan = 100.0 * ( 1.0 - exp( log(alpha) / (double)sn ));
  387.           else if( ac < sn ) {
  388.               facp1 = (double)(ac+(long)1);
  389.               fnmac = (double)(sn-ac);
  390.               quan = 100.0 * qbeta( omalpha, facp1, fnmac );
  391.               }
  392.           else
  393.               quan=100.0;
  394.  
  395.           printf(    "\n&  ...the UQL is %7.3lf%%", quan );
  396.           if( outfile != NULL )
  397.               fprintf( outfile, "\n&  ...the UQL is %7.3lf%%", quan );
  398.  
  399.           yesno = 'y';
  400.           printf("\n\nA/I/U-QLV calculations for the current single-stage plan are complete.");
  401.           printf("\n\nDo you want to specify additional plans [Y|n] ? ");
  402.           fgets( buf, 80, inpfile );
  403.           if( buf[0] != '\0' && buf[0] != '\r' && buf[0] != '\n' )
  404.               sscanf( buf, "%c", &yesno );
  405.           if( yesno == 'n' || yesno == 'N' )
  406.               break;
  407.           }
  408.  
  409.       printf( "\n\nExiting A/I/U-QLV...\n\n" );
  410.       }
  411.  
  412.