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.
