home *** CD-ROM | disk | FTP | other *** search
/ Current Shareware 1994 January / SHAR194.ISO / engineer / ksprob21.zip / KSMISC.EXE / KSPRBAS.BAS < prev    next >
BASIC Source File  |  1992-12-22  |  5KB  |  127 lines
  1.  
  2.                    Probability Programs for Small Computers
  3.                                Joseph C. Hudson
  4.  
  5.         This program is written in Sharp EL-5500III BASIC. Lower case
  6.         letters and extra spaces are used for clarity. The variable L is
  7.         given as a capital letter. sqr is used instead of Sharp's √ and
  8.         pi is used instead of π. If your computer does not predefine pi,
  9.         you will have to define it if you implement line 410. BASIC
  10.         dialects vary among small computers. Expect to make  substantial
  11.         changes in i/o and initialization statements if you use this
  12.         program on another machine.
  13.  
  14.         lines       compute or approximate
  15.  
  16.         100 - 150   binomial density and cdf
  17.         200 - 290   hypergeometric density and cdf
  18.         300 - 350   Poisson density and cdf
  19.         400 - 450   Gaussian cdf
  20.         500 - 570   inverse Gaussian cdf
  21.         600 - 620   Gaussian right tail probability
  22.         630 - 695   subroutines used in lines 600 - 960
  23.         700 - 760   Student's t right tail probability
  24.         800 - 840   chi - square right tail probability
  25.         900 - 960   f distribution right tail probability
  26.  
  27.         errors
  28.  
  29.         the binomial, Poisson and hypergeometric computations are exact
  30.         to 4 decimal places for all values checked. the program does
  31.         well for in all three sections computing normal quantities. for
  32.         Student's t, upper tail probs are acceptable for 4 or more df, 8
  33.         or more df for chi-square. for the f distribution, results are
  34.         acceptable if both df are 8 or more, results are not acceptable
  35.         if both are 4 or less, and are mixed otherwise.
  36.  
  37.         the program crashes if t = 0, x² = df - 1 or f = 1 with both
  38.         degrees of freedom equal. the program will not accept df = 1 for
  39.         Student's t or df < 3 for x².
  40.  
  41.         the file ksprbas.err contains the results of extensive error
  42.         checks. use this for more detail about errors and for checking
  43.         the results of your own implementation.
  44.  
  45.         references: see ks.doc
  46.  
  47. 100 clear : restore : input "x=";x,"n=";n,"p=";p
  48. 110 t = (1-p)^n : s = t
  49. 120 if x = 0 then 140
  50. 130 for i = 1 to x : t = t *(n-i+1) * p / i / (1-p) : s = s + t : next i
  51. 140 using : print x;" ";n;" ";p
  52. 150 print using "###.####";t+0.00005;s+0.00005 : goto 100
  53.  
  54. 200 clear:restore:input "x=";x,"np=";np,"ns=";ns,"k=";k
  55. 210 mi = 0 : if k + ns - np > 0 then let mi = k + ns - np
  56. 220 ma = k : if ns < k then let ma = ns
  57. 230 if x < mi then let t = 0 : let s = 0 : goto 280
  58. 240 if x > ma then let t = 0 : let s = 1 : goto 280
  59. 250 L=np-ns : d=k :nt=ns:if mi > 0 then let L=ns : let d=np-k : let nt=np-ns
  60. 260 t=1:for i=1 to nt:L=L+1:t = t*(L-d)/L : next i : s=t : if x=mi then 280
  61. 270 for i = mi+1 to x : t = t*(k-i+1)*(ns-i+1)/i/(np-k-ns+i) : s = s+t : next i
  62. 280 using : print x;" ";np;" ";ns;" ";k
  63. 290 print using "###.####";t+0.00005;s+0.00005 : goto 200
  64.  
  65.  
  66. ksprbas.bas                                             page 2
  67.  
  68. 300 clear : restore : input "x=";x,"m=";m
  69. 310 t = exp(-m) : s = t
  70. 320 if x = 0 then 340
  71. 330 for i = 1 to x : t = t * m / i : s = s + t : next i
  72. 340 using : print x;" ";m
  73. 350 print using "###.####";t + 0.00005;s + 0.00005 : goto 300
  74.  
  75. 400 clear : restore : input "z=";z
  76. 410 y = exp(-z*z/2) / sqr(2*pi)
  77. 420 x = 1 / (1 + 0.2316419 * abs(z))
  78. 430 p=0.31938153+x*(-0.356563782+x*(1.781477937+x*(-1.821255978+x*1.33027449)))
  79. 440 p = p * y * x : if z > 0 then let p = 1 - p
  80. 450 using : print z;" p=; using "##.####";p + 0.00005 : goto 400
  81.  
  82. 500 clear : restore : input "p=";p
  83. 510 f = 0 : if p > 0.5 then let f = 1 : let p = 1 - p
  84. 520 t = sqr(-2 * ln(p))
  85. 530 n = 2.515517 + t * (0.802853 + t * 0.010328)
  86. 540 d = 1 + t * (1.432788 + t * (0.189269 + t * 0.001308))
  87. 550 z = t - (n / d) + 0.0005 : if f = 1 then let p = 1 - p
  88. 560 if p < 0.5 then let z = - z
  89. 570 using : print p;" z="; using "###.###";z : goto 500
  90.  
  91. 600 gosub 630
  92. 610 input "z=";z : gosub 670
  93. 620 using : print z;" qz="; using "###.####";q+0.00005 : goto 610
  94.  
  95. 630 clear : restore : dim a(9)
  96. 640 data -1.26551223,1.00002368,.37409196,.09678418,-.18628806
  97. 650 data .27886807,-1.13520398,1.48851587,-.82215223,.17087277
  98. 660 for i = 0 to 9 : read a(i) : next i : return
  99.  
  100. 670 y = 1 / (1 + abs(z) / 2 / sqr(2))
  101. 680 p = a(9) : for i = 8 to 0 step -1 : p = a(i) + y * p : next i
  102. 690 q = 0.5 * y * exp(p-z*z/2) : if z < 0 then let q = 1 - q
  103. 695 return
  104.  
  105. 700 gosub 630
  106. 710 input "df=";d,"t=";t : if d < 2 then 710
  107. 720 g = sqr(ln(1 + t * t / d))
  108. 730 s = 0.184 * (8 * d + 3) / d / g
  109. 740 z = sqr(d) * g * (1 - 2 * sqr(1 - exp(-s * s)) / (8 * d + 3)) : gosub 670
  110. 750 if t < 0 then let q = 1 - q
  111. 760 using : print d;t;"qt=";using "###.####";q+.00005 : goto 710
  112.  
  113. 800 gosub 630
  114. 810 input "df=";d,"x=";x:if d < 3 then 810
  115. 820 z = sqr(x + (d - 1) * (ln((d - 1) / x) - 1))
  116. 830 z = z * (x - d + 2 / 3 - 0.08 / d) / abs(x - d + 1) : gosub 670
  117. 840 using : print d;x;"qx=";using "###.####";q + 0.00005 : goto 810
  118.  
  119. 900 gosub 630
  120. 910 input "df1=";d,"df2=";e,"f=";f
  121. 920 p = e / (d * f + e) : c = 0.08 * ((1-p) / e - p/d + (.5-p)/(d+e))
  122. 930 d = d-1 : e = e-1 : z = e * ln(e/((d+e)*p)) + d * ln(d/((d+e)*(1-p)))
  123. 940 z = sqr(3 * (d + e) * z / (3 * (d + e) + 1))
  124. 950 z = z * (e + 1/3 - (d+e+2/3)*p+c) / abs(e - (d+e)*p) : gosub 670
  125. 960 using : print d+1;e+1;f;"qf=";using "###.####";q+0.00005 : goto 910
  126.  
  127.