Converting R code to SAS IML

[JJLee]

I am trying to covert the following R code to SAS IML

LDA<-function(I,r,x,n) {
  total<-sum(r*x)
  options(warn=-1)
  lik<-function(phi) {
    -sum(-r*phi*x+(n-r)*log(1-exp(-phi*x)))
  }
 
  #Maximizing phi
  actual=optim(0.5,lik)$par
 
  if ((1/actual)>x[I]) {
    u=cat(as.integer(1/actual)," ",paste("Warning:Estimate exceeds maximum cells tested per well"),"\n")
  } else {
    u=(as.integer(1/actual))
  }
  print(paste("Estimated from LDA"))
 
  return(c(u,total))
}

For the following test data

I<-9
r<-c(0, 0, 0, 6, 8, 8, 8, 8, 😎
x<-c(50000, 10000, 2000, 400, 80, 16, 3.2, 0.64, 0.128)
n<-c(8, 8, 8, 8, 8, 8, 8, 8, 😎
LDA(I,r,x,n)

The result is 871 and 3199.744.

Here is my IML code.

proc iml;

I=9;
r={0, 0, 0, 6, 8, 8, 8, 8, 8};
x={50000, 10000, 2000, 400, 80, 16, 3.2, 0.64, 0.128};
n={8, 8, 8, 8, 8, 8, 8, 8, 8};

total=sum(r#x);

print I, r, x, n, total;

start lik(phi) global(r, x, n);

f=sum(-r#phi#x+(n-r)#log(1-exp(-phi#x)));

return(f);

finish lik;

phi={0.5};

optn={1, 4, 10};

con={0.0000000001, .};

parm={.,0.1};

call nlpnms(rc, result, "lik", phi, optn, con, , parm);

u=round(1/result,1);

if u>x[1] then print "Warning: Estimate exceeds maximum cells tested per well";

print u, total;

quit;

The result I got is 876 and 3199.744.

What do I need to do in the IML code in order to get the same result as R code?

Thanks,

JJ

[Ksharp]

I noticed in R , object function is 
   -sum(-r*phi*x+(n-r)*log(1-exp(-phi*x)))
in IML, there is not minus .
is it right ?

[Ksharp]

I still get 876 after running Genetic Algorithm.



proc iml;


r={0, 0, 0, 6, 8, 8, 8, 8, 8};
xx={50000, 10000, 2000, 400, 80, 16, 3.2, 0.64, 0.128};
n={8, 8, 8, 8, 8, 8, 8, 8, 8};
total=sum(r#xx);

start func(x) global(r, xx, n);
f=sum(-r#x#xx+(n-r)#log(1-exp(-x#xx)));
return(f);
finish;

id=gasetup(1,1,1234);
call gasetobj(id,1,"func");
call gasetsel(id,100,1,0.95);
call gasetcro(id,0.05,4);
call gasetmut(id,0.05);

call gainit(id,10000,{1e-6,10});
niter = 1000;
summary = j(niter,2);
mattrib summary [c = {"Min Value", "Avg Value"} l=""];
do i = 1 to niter;
 call garegen(id);
 call gagetval(value, id);
 summary[i,1] = value[1];
 summary[i,2] = value[:];
end;
call gagetmem(mem, value, id, 1);
print mem[ l="best member:"],
      "Max Value: " value[l = ""] ;
u=round(1/mem,1);
print u,total;
call gaend(id);
quit;





best member:
0.0011416
Max Value:	-6.520232
u
876
total
3199.744



Maybe there is some little different in object function between R and SAS.

Or you could check which obj has a bigger value.

proc iml;
r={0, 0, 0, 6, 8, 8, 8, 8, 8};
x={50000, 10000, 2000, 400, 80, 16, 3.2, 0.64, 0.128};
n={8, 8, 8, 8, 8, 8, 8, 8, 8};

start lik(phi) global(r, x, n);
f=sum(-r#phi#x+(n-r)#log(1-exp(-phi#x)));
return(f);
finish ;
phi={0.5};
optn={1, 2};
con={1e-6, .};

R=lik(1/871);
SAS=lik(1/876);

print R,SAS;
quit;


r
-6.520349
SAS
40.793882

[Rick_SAS · Accepted solution]

When you perform numerical optimization in statistical software, the algorithms have various "convergence criteria" which tells the software when to stop iterating and report an answer. The criteria differs from algorithm to algorithm and between different software. For this problem, both software are finding the "same" solution, but the solutions differ by about 5E-6, which is a typical range for numerical optimization. Because you are reporting 1/ result, this small difference gets magnified.

In your R program, you can see the actual solution (and/or the reciprocal value) by using

actual=optim(0.5,lik)$par
sprintf("%20.16f",actual)

      [1] " 0.0011474609374996"

In the IML program, you can use 

print result[f=20.16];

      0.0011417389904598

To resolving the issue might be tough. You need to understand the convergence criteria for both software. The SAS/IML documentation includes complete details of the termination criteria and how to control the algorithm...

Personally, I would not pursue trying to resolve the difference. The problem is not in the software. The problem is that you are reporting an ill-conditioned quantity (the reciprocal of a small value) which is highly sensitive to the numerical optimization. A better choice is to report the optimal value itself, not the reciprocal.

[JJLee]

Many thanks for everyone!