I accidentally posted the wrong version of my code, omitting the Q1 loop, and, sure enough that loop generates non-integers:
qfracseq <- seq(from = 1, to = 0, by = -.1) Q <- 10 Qseq <- Q*qfracseq > Qseq [1] 10.00000000000000000 9.00000000000000000 8.00000000000000000 > 7.00000000000000000 [5] 6.00000000000000000 5.00000000000000000 3.99999999999999911 2.99999999999999911 [9] 1.99999999999999956 0.99999999999999978 0.00000000000000000 ----- I'm a tad puzzled by why my code would do this; is this a result of some base-10 fractions not representing as cleanly in base-2 as in -10? At any rate, I would have never thought to check this, thanks much, Shayne On Thu, Sep 6, 2012 at 11:57 AM, William Dunlap <wdun...@tibco.com> wrote: > 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. > [[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.
