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

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

1 ACCEPTED SOLUTION

Accepted Solutions
Rick_SAS
SAS Super FREQ

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.

View solution in original post

14 REPLIES 14
ballardw
Super User

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;

art297
Opal | Level 21

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.

FriedEgg
SAS Employee

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

FriedEgg
SAS Employee

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

PGStats
Opal | Level 21

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
data_null__
Jade | Level 19

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

PGStats
Opal | Level 21

For the sake of completeness, that would be:

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


PG

PG
Rick_SAS
SAS Super FREQ

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.

FriedEgg
SAS Employee

Great post Rick.

art297
Opal | Level 21

I was impressed by all 4 alternative solutions!

FriedEgg
SAS Employee

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

Rick_SAS
SAS Super FREQ

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.

MagCarAddMc
Calcite | Level 5

This solved the issue I was having. Thank you!

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!

What is Bayesian Analysis?

Learn the difference between classical and Bayesian statistical approaches and see a few PROC examples to perform Bayesian analysis in this video.

Find more tutorials on the SAS Users YouTube channel.

Click image to register for webinarClick image to register for webinar

Classroom Training Available!

Select SAS Training centers are offering in-person courses. View upcoming courses for:

View all other training opportunities.

Discussion stats
  • 14 replies
  • 1641 views
  • 7 likes
  • 7 in conversation