this is the code: %let userN = 40; %let p01 = 0.05; %let p02= 0.05; %let p11 = 0.25; %let p12 = 0.25; data s1_A; do n1 = 2 to &userN-2; do m1 = 2 to &userN-2 while ((m1 + n1) < &userN-2); do y1 = 0 to n1-1 ; do x1 = 0 to m1-1 ; term1_p1 = pdf('BINOMIAL', x1, &p11, m1)* pdf('BINOMIAL', y1, &p12, n1); term1_p0 = pdf('BINOMIAL', x1, &p01, m1)* pdf('BINOMIAL', y1, &p02, n1); term1_p0p1 = pdf('BINOMIAL', x1, &p01, m1)* pdf('BINOMIAL', y1, &p12, n1); output; end; end; end; end; run; data s2_B; set s1_A; do n = n1+1 to &userN-(m1+1); do m = m1+1 to &userN-(n1+1) while (n+m <= &userN); do y = y1 to n while (y1<=y<n); do x = x1 to m while (x1<=x<m); term2_p0 = pdf('BINOMIAL', x1, &p01, m1)* pdf('BINOMIAL', y1, &p02, n1) *pdf('BINOMIAL', (x-x1), &p01, (m-m1))* pdf('BINOMIAL', (y-y1), &p02, (n-n1)); term2_p1 = pdf('BINOMIAL', x1, &p11, m1)* pdf('BINOMIAL', y1, &p12, n1) *pdf('BINOMIAL', (x-x1), &p11, (m-m1))* pdf('BINOMIAL', (y-y1), &p12, (n-n1)); term2_p0p1 = pdf('BINOMIAL', x1, &p01, m1)* pdf('BINOMIAL', y1, &p12, n1) *pdf('BINOMIAL', (x-x1), &p01, (m-m1))* pdf('BINOMIAL', (y-y1), &p12, (n-n1)); output; end; end; end; end; run;
... View more