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

- Home
- /
- SAS Programming
- /
- SAS Procedures
- /
- Comb Function Limitations

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

11-05-2014 06:04 PM

I'm trying to use an equation that involves combination calculations (the COMB function in SAS). The equation results in very large numbers, but in the end the value is between 0 and 1. SAS can't handle large values within the COMB function so I keep receiving errors. For example, COMB(1045,600) is acceptable, but SAS doesn't like COMB(1046,600) since I assume it is too large. Anyone know a way to get around this? The equation I'm using is:

COMB(k,i)*COMB(N-k,n-i)/COMB(N,n)

Where k = 1044, i = 600, N = 1046, n = 600

Thanks

Accepted Solutions

Solution

11-06-2014
09:05 AM

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

Posted in reply to MagCarAddMc

11-06-2014 09:05 AM

First, SAS is a case-insensitive language, so N and n represent the same variable in your expression. In the code below I replace n by nn.

Second, one of the basic tenets of numerical analysis is that you should avoid computing factorials and combinatorial expressions in the naive way, because these expressions are prone to overflow.

Third, expressions like this occur often in probability models. Your expression looked familiar to me, and I fortunately recalled that this expresion is the probability density function for the hypergeometric distribution. Thus your problem is solved in SAS by calling the built-in PDF function, like this:

N=1046;

k=1044;

nn = 600;

i = 600;

h = pdf("Hypergeom", i, N, k, nn); /* = COMB(k,i)*COMB(N-k,nn-i)/COMB(N,nn) */

I think an important lesson to learn is that it is helpful for the OP to describe what he is doing and why he is doing it. The context is often useful in guiding those who are trying to answer the question.

All Replies

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

Posted in reply to MagCarAddMc

11-05-2014 06:29 PM

I get Comb(N,n) as 9.176E307. when N=1045. Which is right at the limit of the system. Increasing N>1045 and n=600 exceeds my system capability. (1.798E308)

You can see what your limit is with

data _null_;

x= constant('BIG');

put x;

run;

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

Posted in reply to MagCarAddMc

11-05-2014 06:39 PM

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

Posted in reply to data_null__

11-05-2014 06:50 PM

As always, data_null came up with what I think is an excellent suggestion for addressing your problem. However, in the example you provided, you have an instance where N-n equals zero. That won't make sense for a comb function, but one also can't obtain the log of 0.

Additionally, are you on a system that can differentiate between upper and lower case variable names? You may want to rename n to be n2.

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

Posted in reply to MagCarAddMc

11-05-2014 07:03 PM

This is because the answer goes outside the bounds of the double data type.

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

Posted in reply to FriedEgg

11-05-2014 07:15 PM

proc groovy;

submit;

def binomial(N, K) {

ret = BigInteger.ONE;

for (k = 0; k < K; k++) {

ret = ret.multiply(BigInteger.valueOf(N-k))

.divide(BigInteger.valueOf(k+1));

}

return ret;

}

System.out.println(binomial(1046, 600));

endsubmit;

quit;

215214135164432318969664923845300683858831944672653837461277170046182567034643018347761262897990464760845987368908944618207241161102715206382228014600474147009141057204830063147454042567661669186964749079231644777086378723140490665422937054079891989213282500050827623980250811932790782038439267537485639619200

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

Posted in reply to MagCarAddMc

11-05-2014 08:52 PM

Data_Null_ provided the most useful suggestion, but also:

depending on the relative size of your variables, you may get away with the equivalent expression

comb(N-n, k-i) / comb(N, k) * comb(n, i)

Notice that it is also a good idea to divide before you multiply.

PG

PG

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

Posted in reply to PGStats

11-06-2014 07:33 AM

How would the OPs expression or your equivalent need to be changed to use LCOMB?

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

Posted in reply to data_null__

11-06-2014 12:44 PM

For the sake of completeness, that would be:

**EXP( LCOMB(k,i) + LCOMB(N-k,n-i) - LCOMB(N,n) )**

PG

PG

Solution

11-06-2014
09:05 AM

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

Posted in reply to MagCarAddMc

11-06-2014 09:05 AM

First, SAS is a case-insensitive language, so N and n represent the same variable in your expression. In the code below I replace n by nn.

Second, one of the basic tenets of numerical analysis is that you should avoid computing factorials and combinatorial expressions in the naive way, because these expressions are prone to overflow.

Third, expressions like this occur often in probability models. Your expression looked familiar to me, and I fortunately recalled that this expresion is the probability density function for the hypergeometric distribution. Thus your problem is solved in SAS by calling the built-in PDF function, like this:

N=1046;

k=1044;

nn = 600;

i = 600;

h = pdf("Hypergeom", i, N, k, nn); /* = COMB(k,i)*COMB(N-k,nn-i)/COMB(N,nn) */

I think an important lesson to learn is that it is helpful for the OP to describe what he is doing and why he is doing it. The context is often useful in guiding those who are trying to answer the question.

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

Posted in reply to Rick_SAS

11-06-2014 09:40 AM

Great post Rick.

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

Posted in reply to FriedEgg

11-06-2014 09:49 AM

I was impressed by all 4 alternative solutions!

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

Posted in reply to art297

11-06-2014 10:11 AM

I would say my post was more for amusement than any real practicality.

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

Posted in reply to FriedEgg

11-07-2014 12:47 PM

In the same vein, here is a blog post from 2011 that uses SAS to define and manipulate "BigIntegers" for computing the factorials of numbers when the factorial value exceeds the value that can be stored in a double (1.798E308). I could write a similar program for combinations. Notice that I use the LFACT function, which is the analog to the LCOMB function that Pierre and Data_Null use here. Like Fried Egg's post, this way of solving the problem is more for amusement than for practical scientidic computing, so I think the post uses all of the features that Art likes.

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

Posted in reply to Rick_SAS

11-06-2014 12:54 PM

This solved the issue I was having. Thank you!