BookmarkSubscribeRSS Feed
🔒 This topic is solved and locked. Need further help from the community? Please sign in and ask a new question.
JJLee
Calcite | Level 5

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

1 ACCEPTED SOLUTION

Accepted Solutions
Rick_SAS
SAS Super FREQ

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.

View solution in original post

4 REPLIES 4
Ksharp
Super User
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
Super User
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
SAS Super FREQ

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
Calcite | Level 5

Many thanks for everyone!

sas-innovate-2024.png

Join us for SAS Innovate April 16-19 at the Aria in Las Vegas. Bring the team and save big with our group pricing for a limited time only.

Pre-conference courses and tutorials are filling up fast and are always a sellout. Register today to reserve your seat.

 

Register now!

Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Find more tutorials on the SAS Users YouTube channel.

From The DO Loop
Want more? Visit our blog for more articles like these.
Discussion stats
  • 4 replies
  • 1084 views
  • 4 likes
  • 3 in conversation