cgrandin / csas-latex

Latex style file and template for CSAS Research Documents
0 stars 6 forks source link

Plot: Beverton-Holt stock-recruitment #24

Closed grinnellm closed 6 years ago

grinnellm commented 6 years ago

Here's what I have so far, plotting MPDs for AM2 only, spawning biomass vs recruitment. Still needs the lines. bevertonholt

grinnellm commented 6 years ago

On that note, are the Bev-Holt parameters saved somewhere?

Shawkshaw commented 6 years ago

These should all be in the iscam.rep file so kappa beta ro

jaclyncleary commented 6 years ago

The function I used in my old R code:

plot.curve <- function(sbo, clr) { st=seq(0, max(sbt, sbo), length=100) #correct, sbo if(rectype==1) {

Beverton-Holt

            rrt=kappa*ro*st/(sbo+(kappa-1)*st)*exp(-0.5*tau^2)  
        }
        if(rectype==2)
        {
            #Ricker
            rrt=kappa*ro*st*exp(-log(kappa)*st/sbo)/sbo *exp(-0.5*tau^2) 
        }
        lines(st, rrt, col=clr)
        ro=ro*exp(-0.5*tau^2)
        points(sbo, ro, pch="O", col=clr, cex=1.25)
        points(sbo, ro, pch="+", col=clr, cex=1.25)
    }
    plot.curve(sbo, clr=4);     
    print("blue curve and blue circle+= sbo per SSB")
grinnellm commented 6 years ago

Hmm, I must be missing something, sorry. Can you spell it out for me? The BH equation above: rrt=kapparost/(sbo+(kappa-1)st)exp(-0.5tau^2) looks like it's missing some punctuation? Can you write it with so, kappa, beta, and ro?

grinnellm commented 6 years ago

Nope, I think I have it.

grinnellm commented 6 years ago

bevertonholt

jaclyncleary commented 6 years ago

great. we'll make an assessment biologist out of you yet :)

jaclyncleary commented 6 years ago

Matt- Just double checking with you on eqns, that you plotted: xx = sbt[1:(length(yr)-min(age))] yy = rt

grinnellm commented 6 years ago

You mean that I omitted the first two years of Recruits?

jaclyncleary commented 6 years ago

that you omit the first two years of sbt. the rt vector is already 2 years shorter than sbt. basically we don't have rt values to match with sbt_1951 and sbt_1952 because we would have had to observe them in 1959 and 1950.

grinnellm commented 6 years ago

You got it - I omitted the first two years of rt and sbt, so I use years 1953 to 2017. I just want to confirm my equation is correct (assessment biologist in training). I extract the following variables from the iscam.rep file:

And then I calculate Recruitment as

where Abundance is sbt from the iscam.rep file. Is that correct?

jaclyncleary commented 6 years ago

I think we're saying the same thing. For the points: plot the vector of rt values (vector length is 65 representing 1953-2017) on the y-axis against sbt [1953-2017] from iscam.rep file, removing the first two values of the sbt vector.

Then for the BH curve, use the equation as you've written but use tau in the last term. Ie, exp(-0.5 tau^2) yes abundance == sbt from iscam.rep, and tau is also in the iscam rep file.

Can you also please add a symbol to show sbo vs. ro (both from iscam.rep)?

In looking at my old code, I see the following equation used in the plotting of ro: ro=roexp(-0.5tau^2)

@Shawkshaw - do you know why/ what this adjustment to ro represents?

Add to caption: symbol x corresponds to the MPD estimate of unfished spawning biomass (SB0) and unfished age-2 recruitment R0

grinnellm commented 6 years ago

bevertonholt

grinnellm commented 6 years ago

Getting closer:) Triangles represent MPD of SB_0 vs R_0. Is that the bias correction for the log-transformed data?

jaclyncleary commented 6 years ago

perhaps eh? In the old figs, the "triangles represent MPD of SB_0 vs R_0" fall directly on the line.

grinnellm commented 6 years ago

Yes, I thought it would be on the line too -- might be bug somewhere? This is the current equation:

jaclyncleary commented 6 years ago

the equation is correct. let's see if @Shawkshaw has something to add :)

jaclyncleary commented 6 years ago

@Shawkshaw - clarification: the ro=ro*exp(-0.5tau^2) changes the triangle in the figs so that it's on the line, and I'm not sure why (i.e., what this equation represents, except that Steve did it...)

Shawkshaw commented 6 years ago

Jaclyn, I can't seem to find any reference to this particular equations but if you look at this parameterization of the B-H model Rt=(soBt-k/(1+betaBt-k))exp(residuals-0.5tau^2) and the parameterization of the equation we are using Rt = k ro Bt / (sbo + (k-1) Bt) exp(-0.5 * tau^2) he set the equation up for estimating unfished recruits (ro) by setting everything but the process error (tau) at 0 ...

grinnellm commented 6 years ago

Can we close this? The MPD points are not on the lines in the current version -- is that how you want it @jaclyncleary ? It's easy to apply the correction (above) to put them on the lines if you prefer.

jaclyncleary commented 6 years ago

let's leave them as-is. I don't actually know whether they are supposed to be on the line or not...