Pacific-salmon-assess / samSim

3 stars 0 forks source link

LRP Branch: calcTAC_fixedER.R #12

Closed carrieholt closed 3 years ago

carrieholt commented 3 years ago

This function calls rbeta. I noticed in the other function that calls rbeta, "calcRealCatch()", there is some code to reset the seed:

for whatever reason rbeta appears to draw a different number of randoms

depending on the shape number; until I can figure out why reset unless

running random chains intentionally

if (random != TRUE) { set.seed(setSeedInput) }

However, in my testing, it looks like rbeta always calls 2 random #s for each random-beta draw (same as rnorm). set.seed(1) runif(10) set.seed(1) canER->0.1 cvER->0.42 sigCanER<-cvERcanER shape1<- canER^2 (((1-canER)/sigCanER^2)-(1/canER)) shape2<-shape1 * (1/canER-1) rebeta(1,shape1,shape2) runif(10) You can see that the first 2 random numbers were "used up" in the one rbeta draw. I get the same answer when I change canER and cvER (shape parameters). So, I think as we create new functions with beta calls, we don't need to worry about this. Correct?

krHolt commented 3 years ago

I agree - it seems that we should be ok to use rbeta without worrying about the seed.

seananderson commented 3 years ago

The code for rbeta illustrates where the issue is:

rbeta <- function (n, shape1, shape2, ncp = 0) {
    if (missing(ncp)) 
        .Call(C_rbeta, n, shape1, shape2)
    else {
        X <- rchisq(n, 2 * shape1, ncp = ncp)
        X/(X + rchisq(n, 2 * shape2))
    }
}

It depends on if ncp is set or left at the default ('missing'), but shouldn't depend on the shape1 or 2 values.

The C rbeta calls a uniform random sample twice: https://github.com/wch/r-source/blob/79298c499218846d14500255efd622b5021c10ec/src/nmath/rbeta.c

image

carrieholt commented 3 years ago

Thanks Sean.

CamFreshwater commented 3 years ago

Sorry for not replying earlier! I remember this being a real annoyance for me so glad you two have resolved it. I'll be sure to revisit this the next time I need to use rbeta ;)