nissandjac / smsR

the Stochastic MultiSpecies stock assessment model (SMS)
Other
0 stars 0 forks source link

Something's confusing about logF0 #13

Open mebrooks opened 1 year ago

mebrooks commented 1 year ago

Before I started using the function getFbar, I was using this other way. Now I noticed that they don't always agree and I'm wondering why. In some years they're the same within machine precision, but in other years (e.g. 2013), they're exactly 1.0 different. I can look into it this week, but wanted to document it before I forget.

library(plyr)
FF <- exp(array(sas$reps$value[names(sas$reps$value) == 'logF0'],
                        dim = c(df.tmb$nage, df.tmb$nyears,df.tmb$nseason),
                        dimnames=list(age=ages, year=years, s=1:nseason)))
 FFm <- subset(melt(FF), age %in% 1:2)
 Fbar <- ddply(FFm, ~year, summarize, F=sum(value)/2)

> Fbar$F-getFbar(df.tmb, sas)$F
 [1]  0.000000e+00  1.110223e-16  1.110223e-16  0.000000e+00  5.551115e-17  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
[10] -1.110223e-16  0.000000e+00 -5.551115e-17  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00  0.000000e+00
[19]  2.220446e-16 -1.110223e-16  1.110223e-16  0.000000e+00  0.000000e+00  0.000000e+00  1.000000e+00 -1.110223e-16  0.000000e+00
[28]  0.000000e+00  0.000000e+00  1.000000e+00  1.000000e+00  0.000000e+00  1.000000e+00  1.000000e+00  1.110223e-16  0.000000e+00
[37] -1.110223e-16  0.000000e+00  0.000000e+00  1.000000e+00
mebrooks commented 1 year ago

I think the problem is how zeros are being reported on the log-scale.

nissandjac commented 1 year ago

Yes I think that's it. I only put the Fs that are above 1 into log scale in the c++ code. Then the R code sets them to -Inf before taking the exponent (since exp(-Inf) in R = 0).

I the above code you take the exponent of the zero values which then gives an incorrect F = 1. Probably the correct thing to do in the c++ code would be to assign the zero values to -Inf or NA.