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-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On > Behalf > Of Shayne Hodge > Sent: Thursday, September 06, 2012 11:10 AM > To: r-help@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 > scho...@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@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. ______________________________________________ R-help@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.
