# Normal\_cdf does not behave as expected (comparison with pnorm)

**URL:** <https://discourse.mc-stan.org/t/normal-cdf-does-not-behave-as-expected-comparison-with-pnorm/27127>\
**Category:** RStan\
**Tags:** techniques, rstan\
**Created:** [April 11, 2022, 9:19pm UTC](https://discourse.mc-stan.org/t/normal-cdf-does-not-behave-as-expected-comparison-with-pnorm/27127 "2022-04-11T21:19:53Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![sillymalagg](https://avatars.discourse-cdn.com/v4/letter/s/c6cbf5/32.png) [@sillymalagg](https://discourse.mc-stan.org/u/sillymalagg)\
**Post date:** [April 11, 2022, 9:19pm UTC](https://discourse.mc-stan.org/t/normal-cdf-does-not-behave-as-expected-comparison-with-pnorm/27127/1 "2022-04-11T21:19:53Z")

</div>

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:

```nohighlight
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  
 ![Rplot](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/1/1117b6881f6879ee410fb68f827eed3ed53ce5d9.png)  
 ![Rplot01](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/7/7ff92fde1633be3d70ec3acc120d3d677c65c874.png)  
:

```nohighlight
 // 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!

---

<div class="post-metadata">

**Author:** ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)\
**Post date:** [April 13, 2022, 12:46pm UTC](https://discourse.mc-stan.org/t/normal-cdf-does-not-behave-as-expected-comparison-with-pnorm/27127/2 "2022-04-13T12:46:18Z")

</div>

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

Your R code has:

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

```

But your Stan code only has:

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

```

---

<div class="post-metadata">

**Author:** ![sillymalagg](https://avatars.discourse-cdn.com/v4/letter/s/c6cbf5/32.png) [@sillymalagg](https://discourse.mc-stan.org/u/sillymalagg)\
**Post date:** [April 21, 2022, 12:44am UTC](https://discourse.mc-stan.org/t/normal-cdf-does-not-behave-as-expected-comparison-with-pnorm/27127/3 "2022-04-21T00:44:57Z")

</div>

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…

![Rplot](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/1/1014ef1bd78d8ef94b48e42aea026c34497c6dc2.png)

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;
}
```

---

<div class="post-metadata">

**Author:** ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)\
**Post date:** [April 21, 2022, 12:52am UTC](https://discourse.mc-stan.org/t/normal-cdf-does-not-behave-as-expected-comparison-with-pnorm/27127/4 "2022-04-21T00:52:41Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![sillymalagg](https://avatars.discourse-cdn.com/v4/letter/s/c6cbf5/32.png) [@sillymalagg](https://discourse.mc-stan.org/u/sillymalagg)\
**Post date:** [April 21, 2022, 12:57am UTC](https://discourse.mc-stan.org/t/normal-cdf-does-not-behave-as-expected-comparison-with-pnorm/27127/5 "2022-04-21T00:57:22Z")

</div>

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);

---

<div class="post-metadata">

**Author:** ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)\
**Post date:** [April 21, 2022, 12:59am UTC](https://discourse.mc-stan.org/t/normal-cdf-does-not-behave-as-expected-comparison-with-pnorm/27127/6 "2022-04-21T00:59:26Z")

</div>

Can you post the Stan code for `mean_pop_init`?

---

<div class="post-metadata">

**Author:** ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)\
**Post date:** [April 21, 2022, 1:36am UTC](https://discourse.mc-stan.org/t/normal-cdf-does-not-behave-as-expected-comparison-with-pnorm/27127/7 "2022-04-21T01:36:56Z")

</div>

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:

```r
# 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

---

<div class="post-metadata">

**Author:** ![sillymalagg](https://avatars.discourse-cdn.com/v4/letter/s/c6cbf5/32.png) [@sillymalagg](https://discourse.mc-stan.org/u/sillymalagg)\
**Post date:** [April 21, 2022, 2:30pm UTC](https://discourse.mc-stan.org/t/normal-cdf-does-not-behave-as-expected-comparison-with-pnorm/27127/8 "2022-04-21T14:30:11Z")

</div>

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