turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

Find a Community

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

04-18-2017 10:08 PM

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, 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, 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

Accepted Solutions

Solution

04-19-2017
06:17 PM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

04-19-2017 08:18 AM

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.00114**74609374996"

In the IML program, you can use

print result[f=20.16];

**0.00114**17389904598

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.

All Replies

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

04-19-2017 03:52 AM

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 ?

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

04-19-2017 04:23 AM

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

Solution

04-19-2017
06:17 PM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

04-19-2017 08:18 AM

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.00114**74609374996"

In the IML program, you can use

print result[f=20.16];

**0.00114**17389904598

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.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Highlight
- Email to a Friend
- Report Inappropriate Content

04-19-2017 06:18 PM

Many thanks for everyone!