[R] choose() function returning anomalous results (zero instead of one)

[William Dunlap]

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.