1module igamma; 2 3% Author : Daniel Hobbs , University of Bath, 1995 - 1996 4% 5%-------------------------------------------------------------------------- 6% 7% The incomplete gamma function. 8% 9% igamma!:iter!:series(a,x,iter,sum,last_term) - iteratively computes the 10% value of an approximation to an infinite series used in 11% igamma (for x<=1 or x<a). 12% 13% igamma!:cont!:frac(a,x,iter,iter_max) - iteratively computes the value of 14% the continuous fraction used in igamma (for other values of x). 15% 16% igamma!:eval(a,x) - returns the value at point x of the 17% incomplete gamma function of order ord. 18% 19% The incomplete beta function. 20% 21% ibeta!:cont!:frac(iter,iter_max,a,b,x) - recursively computes 22% the value of the continuous fraction used to 23% approximate to the incomplete beta function. 24% 25% ibeta!:eval(a,b,x) - returns the value of the incomplete beta 26% function with 27% parameters a and b at point x, by approximating to the 28% incomplete beta function using a continued fraction. 29% 30%-------------------------------------------------------------------------- 31 32% Redistribution and use in source and binary forms, with or without 33% modification, are permitted provided that the following conditions are met: 34% 35% * Redistributions of source code must retain the relevant copyright 36% notice, this list of conditions and the following disclaimer. 37% * Redistributions in binary form must reproduce the above copyright 38% notice, this list of conditions and the following disclaimer in the 39% documentation and/or other materials provided with the distribution. 40% 41% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 42% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 43% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 44% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNERS OR 45% CONTRIBUTORS 46% BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 47% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 48% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS 49% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 50% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 51% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 52% POSSIBILITY OF SUCH DAMAGE. 53% 54 55 56% --------------------------- global variables ---------------------------- 57 58fluid '(!:sfiterations); 59 60global '(ibeta!:max!:iter); 61 62algebraic << 63 64%% operator igamma , igamma!:eval, ibeta, ibeta!:eval; 65 66% Set the maximum number of iterations for the continued fraction 67% used in ibeta to be 200. 68 69ibeta!:max!:iter := 200; 70 71%%% Set up rule definitions for igamma and ibeta functions. 72%% 73%%let 74%%{ 75%% igamma(~a,~x) => igamma!:eval(a,x) 76%% when numberp(a) and numberp(x) and a>0 and x>=0 and lisp !*rounded 77%%}; 78%% 79%%let 80%%{ 81%% ibeta(~a,~b,~x) => ibeta!:eval(a,b,x) 82%% when numberp(a) and numberp(b) and numberp(x) and lisp !*rounded 83%% and repart(a)>0 and repart(b)>0 and x>=0 and x<=1 84%%}; 85 86 87% Function igamma!:iter!:series: -- cum_gamma_iter x^i 88% \ ------------- 89% / (a+1)...(a+i) 90% -- i=1 91% Uses Battacharjee's method (1970) (computed recursively). 92 93expr procedure igamma!:iter!:series(a,x,iter,sum,last_term); 94begin 95 scalar value,this_term; 96 97 if (last_term < 10^-(precision(0)+3)) then 98 value := sum 99 else 100 << 101 this_term := (last_term * x / (a+iter)); 102 value := igamma!:iter!:series(a,x,iter+1,sum+this_term,this_term) 103 >>; 104 105 return value; 106end; 107 108 109% Function igamma!:cont!:frac 1 1-a 1 2-a 2 110% --- --- --- --- --- ... 111% x + 1 + x + 1 + x + 112% Recursively computes fraction using 113% Abramowitz and Stegun's method (1964), formula 6.5.31 114% see also DLMF, http://dlmf.nist.gov/8.9.E2 115 116expr procedure igamma!:cont!:frac(a,x,iter,iter_max); 117begin 118 scalar value; 119 120 if (iter>iter_max) then 121 value := 0 122 else 123 value := (iter - a)/ 124 (1 + (iter/ 125 (x + igamma!:cont!:frac(a,x,iter + 1,iter_max)))); 126 127 return value; 128end; 129 130 131% Function igamma!:eval returns the value at point x of the 132% incomplete gamma function with order ord. 133 134expr procedure igamma!:eval(a,x); 135begin 136 scalar arg,frac,last!:frac,acc,value; 137 138 % Decide whether to use a series expansion or a continued fraction. 139 if (x<=1 or x<a+2) then 140 << 141 value := (exp(-x) * x^a) * (1 + igamma!:iter!:series(a,x,1,0,1)) / 142 gamma(a + 1) 143 >> 144 else 145 << 146 % Set required accuracy to be 3 decimal places more than 147 % current precision. 148 acc := 10 ^ -(precision(0)+3); 149 % Obtain a starting value. 150 frac := igamma!:cont!:frac(a,x,1,1); 151 !:sfiterations := 1; 152 % Repeat loop until successive results of continued fraction converge. 153 repeat 154 << 155 !:sfiterations := !:sfiterations + 1; 156 last!:frac := frac; 157 frac := igamma!:cont!:frac(a,x,1,!:sfiterations) 158 >> 159 until (last!:frac - frac) < acc; 160 161 arg := exp(-x) * x^a / gamma(a); 162 value := 1 - arg / (x + frac) 163 >>; 164 165 return value; 166end; 167 168 169% Function ibeta!:cont!:frac: calculates 1 c(2) c(3) 170% --- ---- ---- ... 171% 1 + 1 + 1 + 172% where 173% c(2i) = - (a + i - 1) (b - i) * x 174% --------------------------------- 175% (a + 2i - 2) (a + 2i - 1) (1 - x) 176% and 177% c(2i+1) = i (a + b + i - 1) * x 178% ----------------------------- 179% (a + 2i - 1) (a + 2i) (1 - x) 180 181%expr procedure ibeta!:cont!:frac(iter,iter_max,a,b,x); 182%begin 183% scalar value,c_odd,c_even; 184% 185% if not (fixp(iter) and fixp(iter_max) and numberp(x)) then 186% rederr("ibeta!:cont!:frac called illegally"); 187% 188% if (iter>iter_max) then 189% value := 0 190% else 191% << 192% c_even := -(a+iter-1)*(b-iter)*x / ((a+2*iter-2)*(a+2*iter-1)*(1-x)); 193% c_odd := iter*(a+b+iter-1)*x / ((a+2*iter-1)*(a+2*iter)*(1-x)); 194% value := c_even / 195% (1 + (c_odd / 196% (1 + ibeta!:cont!:frac(iter+1,iter_max,a,b,x)))) 197% >>; 198% 199% return value; 200%end; 201 202% The above version recurses to a depth of iter_max which may be reasonably 203% large. I now provide an alternative version that does the calculation 204% from the inside out and hence avoids that nesting. 205 206expr procedure ibeta!:cont!:frac(iter_first,iter_max,a,b,x); 207begin 208 scalar iter,value,c_odd,c_even; 209 210 if not (fixp(iter_first) and fixp(iter_max) and numberp(x)) then 211 rederr("ibeta!:cont!:frac called illegally"); 212 213 value := 0; 214 iter := iter_max; 215 while iter >= iter_first do << 216 c_even := -(a+iter-1)*(b-iter)*x / ((a+2*iter-2)*(a+2*iter-1)*(1-x)); 217 c_odd := iter*(a+b+iter-1)*x / ((a+2*iter-1)*(a+2*iter)*(1-x)); 218 value := c_even / 219 (1 + (c_odd / 220 (1 + value))); 221 iter := iter - 1 >>; 222 223 return value; 224end; 225 226 227% Function ibeta!:eval: returns the value of the incomplete beta% 228% function with parameters a and b at point x. Method due to Muller (1931). 229 230expr procedure ibeta!:eval(a,b,x); 231begin 232 scalar last_value,value,arg,!:sfiterations; 233 234 if (x=0 or x=1) then 235 value := x 236 else 237 << 238 % 239 if (repart(a+b)-2)*x > (repart(a)-1) then 240 value := 1 - ibeta(b,a,1-x) 241 else 242 << 243 arg := gamma(a+b) * x^a * (1-x)^(b-1) / (a * gamma(a) * gamma(b)); 244 % A starting point of 30 levels of continued fraction. 245 !:sfiterations := 30; 246 % Starting value that will force calculation a second time at least. 247 value := -1; 248 repeat 249 << 250 last_value := value; 251 value := arg * (1/(1 + ibeta!:cont!:frac(1,!:sfiterations,a,b,x))); 252 !:sfiterations := !:sfiterations + 10 253 >> 254 until (abs(value - last_value) < 10^-(precision(0)+3)) 255 or !:sfiterations > ibeta!:max!:iter; 256 >> 257 >>; 258 259 % Error condition should not occur, but in case it does... 260 if !:sfiterations > ibeta!:max!:iter then 261 write 262 "*** Warning: max iteration limit exceeded; result may not be accurate"; 263 264 return value; 265end; 266 267>>; 268 269endmodule; 270 271end; 272 273