R help - Sep 2012 - choose() function returning anomalous results (zero instead of one)

Hello,

(Apologies for length, wanted to get all the relevant detail in that I know
of).

I've been having a lot of trouble with some code for an inventory analysis
problem I was doing, and finally came to the conclusion that it appears
that choose() is returning incorrect values. Specifically:

-------------

Browse[1]> nn
[1] 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3
Browse[1]> Q1
[1] 3
Browse[1]> choose(Q1,3)
[1] 0
Browse[1]> choose(Q1,nn)
[1] 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 0 0 0 0 0 0 0


------------------

Browse[1]> nn
[1] 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3
Browse[1]> Q1
[1] 3
Browse[1]> choose(Q1,nn)
[1] 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 0 0 0 0 0 0 0

--------------

The first was from RStudio, the second from RGui. (Windows x64, version
2.14.2). (That's executed from the point of view of browser() in
ugly_function(), see code below). It's returning choose(3,3) as 0, instead
of 1. It does not do this consistently - I've put the simplest code I've
been able to use to consistently generate the problem beneath my sig.

Something about the loop seems to generate the problem. If I remove the
loop over Q, and just set Q = 3, then choose(Q1,testn) (where testn is
equal to the nn shown above) works without a problem.

I'm not sure if it's related, but I'm also having problems with
browser()
not firing the first time through and/or having code continue to execute
before it interrupts, i.e.:

------
Browse[2]> pprime <- sum(testma)
Error: object 'testma' not found
Browse[2]>
Browse[2]> loopvar <- loopvar + 1
Browse[2]>
debug at #3: partial1 <- (((1 - f11) * (1 - f12))^(Qone - nn)) * ((f11 + (1
-
    f11) * f12)^nn)
Browse[2]> print(proc.time() - ptm)
   user  system elapsed
   0.33    0.11    0.59
-----

Everything was generated by the program, I didn't type any of it in;
browser() is not in the code below, but I generally put it somewhere in
ugly_function() to check on choose(), but it looks like ugly finishes
executing and goes back to body of the program before letting me do
anything.

(Commented out is my problematic solution, using R.basic's nChooseK on the
values most likely to return anomalous results. and bits of debug code).

Shayne Hodge
schodge@ieee.org

# VeridianDynamics Inventory Optimization Script
# Copyright 2012 Shayne Hodge
# Two-vendor Sensitivity Analysis

rm(list=ls()) # Clear previous session data
library("R.basic")
ptm <- proc.time()

ugly_function <- function(nn,k,Qone,Qtwo,f11,f12,f21,f22)
{
partial1 <- (((1-f11)*(1-f12))^(Qone-nn))*((f11+(1-f11)*f12)^nn)
partial2 <- (((1-f21)*(1-f22))^(Qtwo-k))* ((f21+(1-f21)*f22)^(k))
#ifelse ( ((Qone < 125)|((Qone-nn)<75)|(nn<75)), test1 <-
nChooseK(Qone,nn), test1 <- choose(Qone,nn))
#ifelse ( ((Qtwo < 125)|((Qtwo-k)<75)|(k<75)), test2 <-
nChooseK(Qtwo,k),
test2 <- choose(Qtwo,k))
test1 <- choose(Qone,nn)
test2 <- choose(Qtwo,k)
#print(test1[test1==0])
ptest <- ifelse( ((partial1 == 0) | (partial2 == 0)),0,
test1*partial1*partial2*test2)
return(ptest)
}

gen_nk_old <- function (k2,Q1,Q2)
{
nkmatrix <- matrix(nrow = sum(1:k2),ncol=2)
index <- 1
for (n in 0:min(Q1, k2-1)) # This corresponds to the outer sigma for P'
{
for (k in 0:min(Q2, k2-1-n)) # This is the inner sigma for P'
{
nkmatrix[index,1] <- n
nkmatrix[index,2] <- k
index <- index + 1
}
}
return(nkmatrix[(1:(index-1)),])
}

penalty <- 5000000
k2 <- 10

# [ f11 f12 f21 f22 falt v2on alton parton]
f <- c(0.5667, 0.8760, 0.5667, 0.7, 1, 0, 0, 1)

ci <- 1
Q <- 10
Q1 <- 3
Q2 <- Q-Q1
loopvar <- 1

nkmatrix <- gen_nk_old(k2,Q1,Q2)
testn <- nkmatrix[,1]
testk <- nkmatrix[,2]

testma <- ugly_function(testn, testk, Q1, Q2, f[1],f[2],f[3],f[4])
pprime <- sum(testma)

loopvar <- loopvar + 1

print(proc.time() - ptm)

	[[alternative HTML version deleted]]

Your Q1 may be a tad below 3.  When printed it is, by default, rounded
to 7 significant digits.  E.g.,
  > 3 - 3*10^-16
  [1] 3
  > options(digits=17)
  > 3 - 3*10^-16
  [1] 2.9999999999999996
choose will truncate, not round, its inputs to make integers out of them
so your 3 - 3*10^-16 is treated as 2.  
  > choose(3, 0:3)
  [1] 1 3 3 1
  > choose(3-3*10^-16, 0:3)
  [1] 1 3 3 0
  > choose(2, 0:3)
  [1] 1 2 1 0

Try using Q1 <- round(Q1) before calling choose.

Bill Dunlap
Spotfire, TIBCO Software
wdunlap tibco.com

> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at
r-project.org] On Behalf
> Of Shayne Hodge
> Sent: Thursday, September 06, 2012 11:10 AM
> To: r-help at r-project.org
> Subject: [R] choose() function returning anomalous results (zero instead of
one)
> 
> Hello,
> 
> (Apologies for length, wanted to get all the relevant detail in that I know
> of).
> 
> I've been having a lot of trouble with some code for an inventory
analysis
> problem I was doing, and finally came to the conclusion that it appears
> that choose() is returning incorrect values. Specifically:
> 
> -------------
> 
> Browse[1]> nn
> [1] 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3
> Browse[1]> Q1
> [1] 3
> Browse[1]> choose(Q1,3)
> [1] 0
> Browse[1]> choose(Q1,nn)
> [1] 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 0 0 0 0 0 0 0
> 
> 
> ------------------
> 
> Browse[1]> nn
> [1] 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3
> Browse[1]> Q1
> [1] 3
> Browse[1]> choose(Q1,nn)
> [1] 1 1 1 1 1 1 1 1 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 0 0 0 0 0 0 0
> 
> --------------
> 
> The first was from RStudio, the second from RGui. (Windows x64, version
> 2.14.2). (That's executed from the point of view of browser() in
> ugly_function(), see code below). It's returning choose(3,3) as 0,
instead
> of 1. It does not do this consistently - I've put the simplest code
I've
> been able to use to consistently generate the problem beneath my sig.
> 
> Something about the loop seems to generate the problem. If I remove the
> loop over Q, and just set Q = 3, then choose(Q1,testn) (where testn is
> equal to the nn shown above) works without a problem.
> 
> I'm not sure if it's related, but I'm also having problems with
browser()
> not firing the first time through and/or having code continue to execute
> before it interrupts, i.e.:
> 
> ------
> Browse[2]> pprime <- sum(testma)
> Error: object 'testma' not found
> Browse[2]>
> Browse[2]> loopvar <- loopvar + 1
> Browse[2]>
> debug at #3: partial1 <- (((1 - f11) * (1 - f12))^(Qone - nn)) * ((f11 +
(1
> -
>     f11) * f12)^nn)
> Browse[2]> print(proc.time() - ptm)
>    user  system elapsed
>    0.33    0.11    0.59
> -----
> 
> Everything was generated by the program, I didn't type any of it in;
> browser() is not in the code below, but I generally put it somewhere in
> ugly_function() to check on choose(), but it looks like ugly finishes
> executing and goes back to body of the program before letting me do
> anything.
> 
> (Commented out is my problematic solution, using R.basic's nChooseK on
the
> values most likely to return anomalous results. and bits of debug code).
> 
> Shayne Hodge
> schodge at ieee.org
> 
> # VeridianDynamics Inventory Optimization Script
> # Copyright 2012 Shayne Hodge
> # Two-vendor Sensitivity Analysis
> 
> rm(list=ls()) # Clear previous session data
> library("R.basic")
> ptm <- proc.time()
> 
> ugly_function <- function(nn,k,Qone,Qtwo,f11,f12,f21,f22)
> {
> partial1 <- (((1-f11)*(1-f12))^(Qone-nn))*((f11+(1-f11)*f12)^nn)
> partial2 <- (((1-f21)*(1-f22))^(Qtwo-k))* ((f21+(1-f21)*f22)^(k))
> #ifelse ( ((Qone < 125)|((Qone-nn)<75)|(nn<75)), test1 <-
> nChooseK(Qone,nn), test1 <- choose(Qone,nn))
> #ifelse ( ((Qtwo < 125)|((Qtwo-k)<75)|(k<75)), test2 <-
nChooseK(Qtwo,k),
> test2 <- choose(Qtwo,k))
> test1 <- choose(Qone,nn)
> test2 <- choose(Qtwo,k)
> #print(test1[test1==0])
> ptest <- ifelse( ((partial1 == 0) | (partial2 == 0)),0,
> test1*partial1*partial2*test2)
> return(ptest)
> }
> 
> gen_nk_old <- function (k2,Q1,Q2)
> {
> nkmatrix <- matrix(nrow = sum(1:k2),ncol=2)
> index <- 1
> for (n in 0:min(Q1, k2-1)) # This corresponds to the outer sigma for P'
> {
> for (k in 0:min(Q2, k2-1-n)) # This is the inner sigma for P'
> {
> nkmatrix[index,1] <- n
> nkmatrix[index,2] <- k
> index <- index + 1
> }
> }
> return(nkmatrix[(1:(index-1)),])
> }
> 
> penalty <- 5000000
> k2 <- 10
> 
> # [ f11 f12 f21 f22 falt v2on alton parton]
> f <- c(0.5667, 0.8760, 0.5667, 0.7, 1, 0, 0, 1)
> 
> ci <- 1
> Q <- 10
> Q1 <- 3
> Q2 <- Q-Q1
> loopvar <- 1
> 
> nkmatrix <- gen_nk_old(k2,Q1,Q2)
> testn <- nkmatrix[,1]
> testk <- nkmatrix[,2]
> 
> testma <- ugly_function(testn, testk, Q1, Q2, f[1],f[2],f[3],f[4])
> pprime <- sum(testma)
> 
> loopvar <- loopvar + 1
> 
> print(proc.time() - ptm)
> 
> 	[[alternative HTML version deleted]]
> 
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide
http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

On 06/09/2012 2:09 PM, Shayne Hodge wrote:
> Hello,
>
> (Apologies for length, wanted to get all the relevant detail in that I know
> of).
>
> I've been having a lot of trouble with some code for an inventory
analysis
> problem I was doing, and finally came to the conclusion that it appears
> that choose() is returning incorrect values. Specifically:
>
> -------------
>
> Browse[1]> nn
> [1] 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3
> Browse[1]> Q1
> [1] 3
> Browse[1]> choose(Q1,3)
> [1] 0
That's easy to reproduce, with correct behaviour:

 > Q1 <- 3 - 1e-10
 > Q1
[1] 3
 > choose(Q1, 3)
[1] 0

My guess is your code (which is tl;dr) computes Q1 with some rounding 
error, and you're seeing something like this.

Duncan Murdoch