Hi everyone,

I am translating a fisheries stock assessment model from R to stan language for Bayesian fitting. I am trying to translate the following R code:

```
Linf = 115
nlen = 200
L0 = 1
DL = 1
L = L0:nlen
K = 0.15
cv = 0.1
no_at_age = c(0,3e+04,1e+04,5e+03,2e+03,1e+02)
maxage = length(no_at_age)
mean_pop_init <- rep(0, length(no_at_age))
sd_pop_init <- rep(0, length(no_at_age))
pop_init_distr <- rep(0, length(no_at_age))
pop_init <- rep(0, nlen)
for (m in 1:maxage) {
mean_pop_init[m] = Linf - (Linf-L0)*exp(-K*m);
}
for (y in 1:maxage) {
sd_pop_init[y] = cv * mean_pop_init[y];
}
pr_at_age = no_at_age/sum(no_at_age);
for (a in 1:maxage){
for (p in 1:nlen) {
pop_init_distr[p] = pnorm(L[p]+DL,mean_pop_init[a],sd_pop_init[a])-pnorm(L[p],mean_pop_init[a],sd_pop_init[a]);
}
pop_init = pop_init + pr_at_age[a]*pop_init_distr;
}
pop_init = pop_init/sum(pop_init)*sum(no_at_age);
```

See figure attached for the correct outcome of this code. From what I am reading online, in this case pnorm could be substituted by normal_cdf. Surprise surprise, this is not happening. This is what I am writing in stan, I donât report here the initialisation, the data and the parameters, I just report the code where I substituted pnorm with normal_cdf

:

```
// calculate pop_init
for (m in 1:maxage) {
mean_pop_init[m] = Linf - (Linf-L0)*exp(-K*m);
}
for (y in 1:maxage) {
sd_pop_init[y] = cv * mean_pop_init[y];
}
pr_at_age = no_at_age/sum(no_at_age);
for (a in 1:maxage){
for (p in 1:nlen) {
pop_init_distr[p] = normal_cdf(L[p]+DL,mean_pop_init[a],sd_pop_init[a])-normal_cdf(L[p],mean_pop_init[a],sd_pop_init[a]);
}
pop_init = pr_at_age[a]*pop_init_distr;
}
pop_init = pop_init/sum(pop_init)*sum(no_at_age);
```

See figure for the output of this code.

Thank you so much to each one who will spend a bit of time for replying, any help will be very appreciated!

It looks like youâve got a step missing in your Stan code.

Your R code has:

```
pop_init = pop_init + pr_at_age[a]*pop_init_distr;
```

But your Stan code only has:

```
pop_init = pr_at_age[a]*pop_init_distr;
```

1 Like

Hi!

Thank you so much for your reply and sorry if i am getting mack to you quite late. Actually you are right, and your answer helped me to fix part of the problem. Nevertheless, there is another issue and this time I am quite sure is especially related to the normal_cdf function. In fact, the peak of length-distribution does not show up at the mean that I have specified. Check out the plot, the grey histogram is the stan output, while the red line is the desired output. I donât understand what is happeningâŚ

The code remained the same, the problem what you pointed out, but to make it work I had to initialise pop_init with 0 before looping:

```
for (k in 1:nlen) {
pop_init[k] = 0;
}
```

Just to be sure, can you post the full R and Stan code thatâs being used now?

sure!

this is the stan code:

```
for (m in 1:maxage) {
mean_pop_init[m] = Linf - (Linf-L0)*exp(-K*m);
}
for (y in 1:maxage) {
sd_pop_init[y] = cv * mean_pop_init[y];
}
pr_at_age = no_at_age/sum(no_at_age);
for (k in 1:nlen) {
pop_init[k] = 0;
}
for (a in 1:maxage){
for (p in 1:nlen) {
pop_init_distr[p] = normal_cdf(L[p]+DL, mean_pop_init[a],sd_pop_init[a])-normal_cdf(L[p], mean_pop_init[a],sd_pop_init[a]);
}
pop_init = pop_init + pr_at_age[a]*pop_init_distr;
}
pop_init = pop_init/sum(pop_init)*sum(no_at_age);
```

this is the R code:

for (m in 1:maxage) {

mean_pop_init[m] = Linf - (Linf-L0)*exp(-K*m);

}

for (y in 1:maxage) {

sd_pop_init[y] = cv * mean_pop_init[y];

}

pr_at_age = no_at_age/sum(no_at_age);

for (a in 1:maxage){

for (p in 1:nlen) {

pop_init_distr[p] = pnorm(L[p]+DL,mean_pop_init[a],sd_pop_init[a])-pnorm(L[p],mean_pop_init[a],sd_pop_init[a]);

}

pop_init = pop_init + pr_at_age[a]*pop_init_distr;

}

pop_init = pop_init/sum(pop_init)*sum(no_at_age);

Can you post the Stan code for `mean_pop_init`

?

Ah ignore me, it was already in your code!

As far as I can see those two code snippets look the same. The next step is to debug whether it is the `normal_cdf`

function or some other aspect of the Stan model. The easiest way to do this is to export the `normal_cdf`

function to R and use it in place of `pnorm`

, and see whether you get the same results (with no other changes).

The easiest way to do this is through the `stanFunction`

framework in the `StanHeaders`

package, which allows you to export a Stan function to R. In your case you would run:

```
# Need to provide 'dummy' initial values to define the needed argument types
StanHeaders::stanFunction("normal_cdf", y = 0, mu = 0, sigma = 1)
```

Then you can just replace `pnorm`

with `normal_cdf`

in your R code. If the results of the R code with `pnorm`

differ from those with `normal_cdf`

, then that gives us a good starting point

It is the first loop, mean_pop_init is just a quantity calculated starting from some parameters (K, L0, Linf and m)

It is initialised as:

mean_pop_init â rep(0, length(no_at_age))

Also, when printed, it returns the desired values: 16.87929 30.54672 42.31039 52.43547 61.15021 68.65106