Hi.
This isn't specifically on O.R. question, but it has to do with a kind of numerical computing that perhaps you'll have some expertise with.
I'm looking for advice / feedback / abuse concerning a numerical computation issue. Briefly, I'm trying to solve the equation y=f(x) for the value of x, given the function f() and the value y.
I have provided some code below for a function that depends on three parameters, N, M, and r, with a target value of G. I have fixed N = 730 and M = 1, and want to find the value of r that gives G = 730. (The code is designed to look at a range of values of N, with G changing proportionally, but in this example I didn't vary N.) The function has monotonicity properties that guarantee the value of r will fall between 0 and 1; when r = 0, G = 0, and G increases as r increases until G ~ 2**N when r = 1. The solution method is a simple bisection that moves r down if G is above the target, and moves r up if G is below the target, until either the difference between the computed G and the target is within a certain tolerance level, or an upper limit of iterations is reached.
I have run the code below in SAS 9.1.3 on both Windows XP and HP-UX. I get slightly different results on the two systems, but in both cases there is a discontinuity near the solution; as the value of r converges to 0.00929, the value of G flip-flops between 753 and 540, but never seems to get any closer to the target value of 730.
This is a tricky function. It is a sum of about N terms (recall N = 730), and each term involves the product of a very large number (comb(N,k)) with a very small number.
Is there a better way to do this computation to converge on the target value, or is this just too much to ask of SAS?
************************************************************************
%let Nlo = 730 ;
%let Nhi = 730 ;
%let Ninc = 10 ;
%let N0 = 1000 ;
%let G0 = 1000 ;
%let M0 = 1 ;
%let rstart = 0.5 ;
%let gtol = 0.0001 ;
%let ntries = 100 ;
%let bouncelim = 4 ;
%let vnum = 01 ;
%*let vnum = %sysfunc(putn(%sysevalf(&vnum. + 1),z2.)) ;
%put vnum = &vnum. ;
%let dnum = 01 ;
%*let dnum = %sysfunc(putn(%sysevalf(&dnum. + 1),z2.)) ;
%put dnum = &dnum. ;
* trial and error to find r that fits N, M, G for several scaled N, G, with debug info ;
data rbyn_&vnum._d&dnum. ;
keep N M G r rfound i rhi rlo groups gdif rdif mingdif rmingdif gmingdif GBOUNCE ;
array combnk{2:&N0.} ;
M = &M0. ;
do N = &Nlo. to &Nhi. by &Ninc. ;
G = N * %sysevalf(&G0. / &N0.) ;
do k = 2 to N ;
combnk{k} = comb(N,k) ;
end ; * k ;
rfound = 0 ;
r = &rstart. ;
rhi = 1 ;
rlo = 0 ;
pgdif = . ;
mingdif = . ;
rmingdif = . ;
gmingdif = . ;
gbounce = 0 ;
do i = 1 to &ntries. while (rfound = 0) ;
groups = 0 ;
do k = 2 to N ;
qk = (1 - r**k)**M ;
pk = 1 - qk ;
groups_k = combnk{k} * pk ;
groups + groups_k ;
end ; * k ;
gdif = (groups - G) ;
absgdif = abs(gdif) ;
rdif = (rhi - rlo) ;
output ;
if (absgdif < abs(mingdif)) or missing(mingdif) then do ;
if abs(sum(gdif,-mingdif)) < %sysevalf(0.5 * >ol.) then do ;
gbounce + 1 ;
end ;
else do ;
gbounce = 0 ;
end ;
mingdif = gdif ;
rmingdif = r ;
gmingdif = groups ;
end ;
else if ((absgdif > abs(mingdif)) and (gdif * mingdif > 0)) or (gbounce ge &bouncelim.) then do ; * should be impossible ;
rfound = -1 ;
r = rmingdif ;
groups = gmingdif ;
end ;
if absgdif < >ol. then do ;
rfound = 1 ;
end ;
else if groups > G then do ;
rhi = r ;
if . < pgdif < 0 then do ;
sumdiff = (gdif - pgdif) ;
r = ((gdif * rlo) - (pgdif * rhi)) / sumdiff ;
end ;
else do ;
r = 0.5 * (rhi + rlo) ;
end ;
end ;
else if groups < G then do ;
rlo = r ;
if pgdif > 0 then do ;
sumdiff = (pgdif - gdif) ;
r = ((pgdif * rlo) - (gdif * rhi)) / sumdiff ;
end ;
else do ;
r = 0.5 * (rlo + rhi) ;
end ;
end ;
pgdif = gdif ;
* output ;
end ; * i ;
output ;
end ; * N ;
run ;
************************************************************************
Thanks!
-- TMK --
"The Macro Klutz"