# Keep getting infinite gradient although likelihood is differentiable

**URL:** <https://discourse.mc-stan.org/t/keep-getting-infinite-gradient-although-likelihood-is-differentiable/41541>\
**Category:** Modeling\
**Tags:** initialization, rstan, fitting-issues\
**Created:** [August 21, 2026, 8:10pm UTC](https://discourse.mc-stan.org/t/keep-getting-infinite-gradient-although-likelihood-is-differentiable/41541 "2026-08-21T20:10:27Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![SeHellmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sehellmann/32/21103_2.png) [@SeHellmann](https://discourse.mc-stan.org/u/SeHellmann)\
**Post date:** [August 21, 2026, 8:10pm UTC](https://discourse.mc-stan.org/t/keep-getting-infinite-gradient-although-likelihood-is-differentiable/41541/1 "2026-08-21T20:10:27Z")

</div>

I am trying to fit a model to a couple (19) datasets. For 16 datasets, the model runs smoothly (no divergences, reasonable Rhat, and very good ESSs). But it seems, three datasets are a bit hard to estimate. I determined the best possible initial values to prevent infinite log-likelihoods, however, I cannot get rid of the “Gradient evaluated at the initial value is not finite.” error.

This is my model (which is cumulative prospect theory):

```stan
functions {
  real prob_weighting ( real p, real gamma) {
    return exp(-(-log(p))^gamma);
  }
}

data {
  ////////////////////////
  // Context variables:
  int<lower=0> N;
  int<lower=0> Ntrials;
  int<lower=0> S1; // number of states of first option
  int<lower=0> S2; // number of states of second option
  int<lower=1, upper=N> sbj[Ntrials]; // subject identifier
  vector[4] prior_mu;
  vector[4] prior_sd;

  
  ////////////////////////
  // Predictor variables:
  
  // IMPORTANT: Note, that outcomes and corresponding probabilities are expected 
  // to be sorted according to the outcomes within each option and trial!
  // Positive outcomes should come first and be sorted decreasingly, losses 
  // should come second and be sorted increasingly!
  // In addition, probabilities should already be included as CUMULATIVE Probabilities!
  vector[S1] outcomes1[Ntrials]; //values of the outcomes of the first option in each trial
  vector[S1] cumprobs1[Ntrials]; // cumulative outcome probabilities of first option in each trial
  vector[S2] outcomes2[Ntrials]; //values of the outcomes of the second option in each trial
  vector[S2] cumprobs2[Ntrials]; // cumulative outcome probabilities of second option in each trial
  
  ////////////////////////
  // Outcome variables: 
  int<lower=0, upper=1> choice[Ntrials]; // Final choice: 1=chooseX; 2 = chooseY

}

parameters {

  // population means
  real mu_log_dscale;
  real mu_log_alpha;
  real mu_log_gamma;
  real mu_log_lambda;

  // population SDs
  real<lower=0> sd_log_dscale;
  real<lower=0> sd_log_alpha;
  real<lower=0> sd_log_gamma;
  real<lower=0> sd_log_lambda;

  // non-centered subject effects
  vector[N] z_log_dscale;
  vector[N] z_log_alpha;
  vector[N] z_log_gamma;
  vector[N] z_log_lambda;
}

transformed parameters {

    // subject parameters
  vector<lower=0>[N] dscale;
  vector<lower=0>[N] alphas;
  vector<lower=0>[N] gammas;
  vector<lower=0>[N] lambdas;
  dscale = log1p_exp(mu_log_dscale + sd_log_dscale * z_log_dscale);
  alphas = log1p_exp(mu_log_alpha + sd_log_alpha * z_log_alpha);
  gammas = log1p_exp(mu_log_gamma + sd_log_gamma * z_log_gamma);
  lambdas = log1p_exp(mu_log_lambda + sd_log_lambda * z_log_lambda);
}

model {
  // population means
  mu_log_dscale ~ normal(prior_mu[1] , prior_sd[1]);
  mu_log_alpha ~ normal(prior_mu[2] , prior_sd[2]);
  mu_log_gamma ~ normal(prior_mu[3] , prior_sd[3]);
  mu_log_lambda ~ normal(prior_mu[4] , prior_sd[4]);
  
  // population SDs
  sd_log_dscale ~ normal(0, 1);
  sd_log_alpha ~ normal(0, 1);
  sd_log_gamma ~ normal(0, 1);
  sd_log_lambda ~ normal(0, 1);

  // non-centered latent effects
  z_log_dscale ~ std_normal();
  z_log_alpha ~ std_normal();
  z_log_gamma ~ std_normal();
  z_log_lambda ~ std_normal();

  
  // Likelihood
  real drift;
  int curM;
  vector[S1] cur_outcomes1;
  vector[S2] cur_outcomes2;
  vector[S1] weights1;
  vector[S2] weights2;
  real prev_weighted_cumprob, weighted_cumprob; 
  int switchsign;
  
  for (j in 1:Ntrials){
    cur_outcomes1 = outcomes1[j];
    cur_outcomes2 = outcomes2[j];
    switchsign = 0;
    prev_weighted_cumprob = 0.0; // reset prop_weight
    for (out in 1:S1) {
      if (cur_outcomes1[out] > 0) {
       weighted_cumprob = prob_weighting(cumprobs1[j][out], gammas[sbj[j]]);
        weights1[out] = weighted_cumprob - prev_weighted_cumprob;
        prev_weighted_cumprob = weighted_cumprob;
        
        cur_outcomes1[out] = cur_outcomes1[out]^alphas[sbj[j]];
      } else {
        if (switchsign==0) {
          prev_weighted_cumprob = 0.0;
          switchsign=1;
        }
        weighted_cumprob = prob_weighting(cumprobs1[j][out], gammas[sbj[j]]);
        weights1[out] = weighted_cumprob - prev_weighted_cumprob;
        prev_weighted_cumprob = weighted_cumprob;
        
        cur_outcomes1[out] = - lambdas[sbj[j]] * (-cur_outcomes1[out])^alphas[sbj[j]];
      }
    }
    
    switchsign = 0;
    prev_weighted_cumprob = 0.0;  
    for (out in 1:S2) {
      if (cur_outcomes2[out] > 0) {
        weighted_cumprob = prob_weighting(cumprobs2[j][out], gammas[sbj[j]]);
        weights2[out] = weighted_cumprob - prev_weighted_cumprob;
        prev_weighted_cumprob = weighted_cumprob;
        
        cur_outcomes2[out] = cur_outcomes2[out]^alphas[sbj[j]];
      } else {
        if (switchsign==0) {
          prev_weighted_cumprob = 0.0;
          switchsign=1;
        }
         weighted_cumprob = prob_weighting(cumprobs2[j][out], gammas[sbj[j]]);
        weights2[out] = weighted_cumprob - prev_weighted_cumprob;
        prev_weighted_cumprob = weighted_cumprob;
        
        cur_outcomes2[out] = - lambdas[sbj[j]] * (-cur_outcomes2[out])^alphas[sbj[j]];
      }
    }
    drift = dscale[sbj[j]] * 
      (sum( cur_outcomes1 .* weights1)-
      sum( cur_outcomes2 .* weights2));
    
  // Choice probability
    choice[j] ~ bernoulli_logit(drift);
  }
}

```

There you see that I even started tuning my priors to be suitable for the data. With this, I could actually fit two of the three problematic datasets, but the last one still resists all my attempts.  
It contains about 270 choices from 64 participants each (in total 17320 observations).  
I used the mean of participant-level maximum likelihood estimates for the parameters for the initialization of the group-level priors for the means (mu\_…) their standard deviation for the sd\_…'s, and the scaled individual offsets for the random effects (z\_…), so the initials should be fine actually.

Also, pooling all participants by ignoring the hierarchical structure and investigating the likelihood surface for the population means produces reasonably smooth curves.

I attach the problematic sub-data and a minimal R-script to reproduce the error.

Is there something problematic with my model or any strategy that I could use to find out what exactly is going wrong?

[Fit\_choice\_data.R](https://discourse.mc-stan.org/uploads/short-url/1HdIKJBdu6ZOQ6V6iZm6HSZrrX.R) (9.0 KB)

[choice\_data.csv](https://discourse.mc-stan.org/uploads/short-url/nXaiW5KmqJSH36icP6XadNWDmsR.csv) (780.6 KB)

---

<div class="post-metadata">

**Author:** ![kaskogsholm](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/kaskogsholm/32/19209_2.png) [@kaskogsholm](https://discourse.mc-stan.org/u/kaskogsholm)\
**Post date:** [August 24, 2026, 2:56pm UTC](https://discourse.mc-stan.org/t/keep-getting-infinite-gradient-although-likelihood-is-differentiable/41541/2 "2026-08-24T14:56:48Z")

</div>

Welcome to the forums SeHellmann! Thank you for including your code and data.

I think the root cause is your prob\_weighting function. It is defined for p=0, but its gradient is not well behaved. When I altered the function to include a guard, the initialization issue was resolved (in my testing):

```Stan
functions {
  real prob_weighting(real p, real gamma) {
    if (p == 0) {
      return 0;
    } else if (p == 1) {
      return 1;
    }

    return exp(-(-log(p))^gamma);
  }
}

```

The issue is data-dependent because the dataset includes many exact 0s for the argument p, and I presume the others were not as extreme.

---

<div class="post-metadata">

**Author:** ![SeHellmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sehellmann/32/21103_2.png) [@SeHellmann](https://discourse.mc-stan.org/u/SeHellmann)\
**Post date:** [August 24, 2026, 8:26pm UTC](https://discourse.mc-stan.org/t/keep-getting-infinite-gradient-although-likelihood-is-differentiable/41541/3 "2026-08-24T20:26:49Z")

</div>

> [@kaskogsholm](#):
>
> k the root cause is your prob\_weighting function. It is defined for p=0, but its gradient is not well behaved

I wished, I had posted earlier. Thanks a lot! It seems to work just fine. Such a simple solution, that would’ve saved me quite some time.  
I checked and visualized the gradients just today, but I expected that defining the border cases “by hand” would not solve the differentiability issue, since the function’s behaviour does not really improve. Also, there are other datasets with much more 0’s. Anyways, now the model runs without problems.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [August 30, 2026, 5:32pm UTC](https://discourse.mc-stan.org/t/keep-getting-infinite-gradient-although-likelihood-is-differentiable/41541/4 "2026-08-30T17:32:41Z")

</div>

I thought this would be an instructive example to rewrite. I think the result is much more readable. I removed most of the documentation, which was redundant with the good naming conventions already in use. If I were documenting this, I’d include a link to a paper or somewhere that describes the model if such a place exists. I did remove the redundant `_log` markings on all the variables and also compressed the outcomes and cumulative probabilities into a 2D array. I also switched to the new array syntax.

In the end, this won’t be any more efficient, but I hope it’s easier to read for intent.

```stan
functions {
  real prob_weighting(real p, real gamma) {
    return (p == 0 || p == 1) ? p : exp(-(-log(p))^gamma);
  }

  vector center_neg_prob(real mu, real sd, vector z) {
    return log1p_exp(mu + sd * z); // equiv., -log_inv_logit(mu + sd * z);
  }    

  real unscaled_drift(vector outcomes, vector cumprobs, real lambda, real alpha, real gamma) {
    int S = rows(outcomes);
    int switchsign = 0;
    real prev_weighted_cumprob = 0.0;
    real lp = 0;
    for (s in 1:S) {
      real weighted_cumprob = prob_weighting(cumprobs[s], gamma);
      if (outcomes[s] > 0) {
        real weight = weighted_cumprob - prev_weighted_cumprob;
        prev_weighted_cumprob = weighted_cumprob;
        real cur_outcome = outcomes[s]^alpha;
        lp += weight * cur_outcome;
      } else {
        if (!switchsign) {
          prev_weighted_cumprob = 0;
          switchsign = 1;
        }
        real weight = weighted_cumprob - prev_weighted_cumprob;
        prev_weighted_cumprob = weighted_cumprob;
        real cur_outcome = -lambda * (-outcomes[s])^alpha;
        lp += weight * cur_outcome;
      }
    }
    return lp;
  }
  void validate_decreasing(array[] vector xs) {
   for (x in xs) {
     for (n in 2:rows(x)) {
        if (x[n] > x[n - 1]) {
          reject("require outcomes in decreasing order");
        }
      }      
    }
  }    
  void validate_increasing(array[] vector xs) {
    for (x in xs) {
      for (n in 2:rows(x)) {
        if (x[n] < x[n - 1]) {
          reject("require cumulative probabilities in increasing order");
        } 
      }
    }
  }    
}
data {
  int<lower=0> N, Ntrials;
  int<lower=0> S1, S2;
  array[Ntrials] int<lower=1, upper=N> sbj;
  vector[4] prior_mu;
  vector<lower=0>[4] prior_sd;
  array[2, Ntrials] vector[S1] outcomes;
  array[2, Ntrials] vector<lower=0, upper=1>[S1] cumprobs;
  array[Ntrials] int<lower=0, upper=1> choice; // 1 = chooseX, 2 = chooseY
}
transformed data {
  validate_decreasing(outcomes[1]);
  validate_decreasing(outcomes[2]);
  validate_increasing(cumprobs[1]);
  validate_increasing(cumprobs[2]);
}
parameters {
  real mu_dscale, mu_alpha, mu_gamma, mu_lambda;
  real<lower=0> sd_dscale, sd_alpha, sd_gamma, sd_lambda;
  vector[N] z_dscale, z_alpha, z_gamma, z_lambda;
}
transformed parameters {
  vector<lower=0>[N] dscale = center_neg_prob(mu_dscale, sd_dscale, z_dscale);
  vector<lower=0>[N] alphas = center_neg_prob(mu_alpha, sd_alpha, z_alpha);
  vector<lower=0>[N] gammas = center_neg_prob(mu_gamma, sd_gamma, z_gamma);
  vector<lower=0>[N] lambdas = center_neg_prob(mu_lambda, sd_lambda, z_lambda);
}
model {
  {mu_dscale, mu_alpha, mu_gamma, mu_lambda} ~ normal(prior_mu, prior_sd);
  {sd_dscale, sd_alpha, sd_gamma, sd_lambda} ~ std_normal();
  z_dscale ~ std_normal();
  z_alpha ~ std_normal();
  z_gamma ~ std_normal();
  z_lambda ~ std_normal();

  vector[Ntrials] drift;
  for (j in 1:Ntrials) {
    int s = sbj[j];
    real d1 = unscaled_drift(outcomes[1, j], cumprobs[1, j], lambdas[s], alphas[s], gammas[s]);
    real d2 = unscaled_drift(outcomes[2, j], cumprobs[1, j], lambdas[s], alphas[s], gammas[s]);
    drift[j] = d1 - d2;
  }
  choice ~ bernoulli_logit(dscale .* drift);
}

```

I refactored it a piece at a time and checked that it compiled, but I did not check that it gets the same answer. The efficiency could probably be improved by vectorizing the operations within the `unscaled_drift` function, but I didn’t want to think too hard about this. It would be easier if the `cur_outcome` definition wasn’t different in the two branches, but we can’t change the model to make it easier to compute :-).

---

<div class="post-metadata">

**Author:** ![SeHellmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sehellmann/32/21103_2.png) [@SeHellmann](https://discourse.mc-stan.org/u/SeHellmann)\
**Post date:** [August 31, 2026, 4:01pm UTC](https://discourse.mc-stan.org/t/keep-getting-infinite-gradient-although-likelihood-is-differentiable/41541/5 "2026-08-31T16:01:43Z")

</div>

Thank you so much for putting the effort in this! This reads indeed much simpler. I had to undo some of the changes, because in general S1 is not S2, that means the two lotteries potentially differ in the number of outcomes. But I still kept many of your re-arrangements and the model looks much tiidier with it.

Also, I used the rather complicated way to define the log-likelihood as a vector because I compared two models using loo, so I originally had it in the output as transformed parameter. I just simplified it for the purpose of sharing a reproducible example. (Apparently, I should have it in the generated quantities block, but I never understood why this should be faster, because the computations need to be done anyways, but this would be a thing for a different question.)

---

<div class="post-metadata">

**Author:** ![SeHellmann](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/sehellmann/32/21103_2.png) [@SeHellmann](https://discourse.mc-stan.org/u/SeHellmann)\
**Post date:** [September 1, 2026, 9:22am UTC](https://discourse.mc-stan.org/t/keep-getting-infinite-gradient-although-likelihood-is-differentiable/41541/6 "2026-09-01T09:22:55Z")

</div>

I probably should open a new question for this, but I came across a different issue, now. For most datasets, the sampling runs smoothly, but for one dataset the chains use up all my memory (which is 1.5TB so it should be enough).  
I have the same model (used most of your refactoring @Bob_Carpenter), but instead of only modelling choices, I use the wiener\_lpdf to model choices and response times. I think, vectorization is a bit difficult because the wiener\_lpdf does not support choice-directed response times.

```stan
functions {
  real prob_weighting(real p, real gamma) {
    return (p == 0 || p == 1) ? p : exp(-(-log(p))^gamma);
  }

  vector center_neg_prob(real mu, real sd, vector z) {
    return log1p_exp(mu + sd * z); // equiv., -log_inv_logit(mu + sd * z);
  }    

  real unscaled_drift(vector outcomes, vector cumprobs, real lambda, real alpha, real gamma) {
    int S = rows(outcomes);
    int switchsign = 0;
    real prev_weighted_cumprob = 0.0;
    real lp = 0;
    for (s in 1:S) {
      real weighted_cumprob = prob_weighting(cumprobs[s], gamma);
      if (outcomes[s] > 0) {
        real weight = weighted_cumprob - prev_weighted_cumprob;
        prev_weighted_cumprob = weighted_cumprob;
        real cur_outcome = outcomes[s]^alpha;
        lp += weight * cur_outcome;
      } else {
        if (!switchsign) {
          prev_weighted_cumprob = 0;
          switchsign = 1;
        }
        real weight = weighted_cumprob - prev_weighted_cumprob;
        prev_weighted_cumprob = weighted_cumprob;
        real cur_outcome = -lambda * (-outcomes[s])^alpha;
        lp += weight * cur_outcome;
      }
    }
    return lp;
  }
}

data {
  int<lower=0> N, Ntrials;
  int<lower=0> S1, S2;
  array[Ntrials] int<lower=1, upper=N> sbj;
  vector[S1] outcomes1[Ntrials]; //values of the outcomes of the first option in each trial
  vector[S1] cumprobs1[Ntrials]; // cumulative outcome probabilities of first option in each trial
  vector[S2] outcomes2[Ntrials]; //values of the outcomes of the second option in each trial
  vector[S2] cumprobs2[Ntrials]; // cumulative outcome probabilities of second option in each trial
    // Outcome variables: 
  real<lower=0> RTs[Ntrials];
  real<lower=0> minRTs[N];
  int<lower=1, upper=2> choice[Ntrials]; // Final choice: 1=chooseX; 2 = chooseY
}

parameters {
  // population means
  real mu_log_A, mu_log_dscale, mu_logit_z, mu_logit_tau,mu_log_alpha,mu_log_gamma,mu_log_lambda;
  // population SDs
  real<lower=0> sd_log_A, sd_log_dscale, sd_logit_z, sd_logit_tau, sd_log_alpha, sd_log_gamma, sd_log_lambda;
  // non-centered subject effects
  vector[N] z_log_A, z_log_dscale, z_logit_z, z_logit_tau, z_log_alpha, z_log_gamma, z_log_lambda;
}

transformed parameters {
    // subject parameters
  vector<lower=0>[N] As = log1p_exp(mu_log_A + sd_log_A * z_log_A);
  vector<lower=0>[N] dscale = log1p_exp(mu_log_dscale + sd_log_dscale * z_log_dscale);
  vector<lower=0,upper=1>[N] Zs = inv_logit(mu_logit_z + sd_logit_z * z_logit_z);
  vector<lower=0,upper=1>[N] tau = inv_logit(mu_logit_tau + sd_logit_tau * z_logit_tau);
  vector<lower=0>[N] alphas = log1p_exp(mu_log_alpha + sd_log_alpha * z_log_alpha);
  vector<lower=0>[N] gammas = log1p_exp(mu_log_gamma + sd_log_gamma * z_log_gamma);
  vector<lower=0>[N] lambdas = log1p_exp(mu_log_lambda + sd_log_lambda * z_log_lambda);
}

model {
  // population means
  mu_log_A ~ normal(1 , 10);
  mu_log_dscale ~ normal(0 , 10);
  mu_logit_z ~ normal(0 , 1);
  mu_logit_tau ~ normal(0 , 2);
  mu_log_alpha ~ normal(0.3 , 1);
  mu_log_gamma ~ normal(0.3 , 1);
  mu_log_lambda ~ normal(0.5 , 2);
  
  // population SDs
  {sd_log_A, sd_log_dscale, sd_logit_z, sd_logit_tau, sd_log_alpha, sd_log_gamma, sd_log_lambda} ~ std_normal(); 

  // non-centered latent effects
  z_log_A ~ std_normal(); 
  z_log_dscale ~ std_normal(); 
  z_logit_z ~ std_normal();
  z_logit_tau ~ std_normal();
  z_log_alpha ~ std_normal();
  z_log_gamma ~ std_normal();
  z_log_lambda ~ std_normal();

  for (j in 1:Ntrials) {
    int s = sbj[j];
    real d1 = unscaled_drift(outcomes1[j], cumprobs1[j], lambdas[s], alphas[s], gammas[s]);
    real d2 = unscaled_drift(outcomes2[j], cumprobs2[j], lambdas[s], alphas[s], gammas[s]);
    real drift = dscale[s] * (d1 - d2);
    if (choice[j]==1) { // upper boundary
      target += wiener_lpdf(RTs[j] | As[s], tau[s]*minRTs[s], Zs[s], drift);
    } else {
      target += wiener_lpdf(RTs[j] | As[s], tau[s]*minRTs[s],1-Zs[s], -drift);
    }
  }
}

```

I am using rstan for this with following options:  
This script fits the CPT-DDM model to the different data sets using the stan

# 

```r
## This script fits the CPT-DDM model to the different data sets using the stan
# models included in the folder "stan_models"

###### Compile the models ######
rstan::rstan_options(auto_write = TRUE)
fit_model = rstan::stan_model("stan_models/stan_riskyDDM_prelec_withlosses_transformation_refactored.stan")

###### General MCMC simulation settings ######
niter = 3000; nwarmup = 400; nchain = 10; ncore = nchain 
nthin = 4; adapt_delta = 0.93
stepsize = 1; max_treedepth = 12

fitted_RiskyDDM = rstan::sampling(
   fit_model, data = dataList, 
   init=make_inits, 
   warmup = nwarmup, iter = niter, chains = nchain, cores = ncore, thin = nthin,
   control = list(adapt_delta = adapt_delta, max_treedepth = max_treedepth, 
                  stepsize = stepsize),
   pars = c("z_log_A", "z_logit_z", "z_logit_tau", "z_log_dscale", "z_log_alpha", "z_log_gamma", rep("z_log_lambda", has_losses)),
   include = FALSE,
   seed=1234)
  )
  

```

As far as I understand, when people recommending CmdStan instead of rstan mostly solves problems, when the chains are finished and the results are parsed into R. But for me, the chains already use all my memory before reaching the first iteration. I ran a single chain and the memory usage converges towards roughly 200GB, which still seems too large for me. The problematic dataset has 103 subjects and a total of 21,740 trials (which is much smaller than the largest dataset with 537 subjects and 113,270 trials, which takes a lot of time of course but uses reasonable memory space; about 800MB per chain when it is running). Would it still make sense to try CmdStan?

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [September 8, 2026, 8:24pm UTC](https://discourse.mc-stan.org/t/keep-getting-infinite-gradient-although-likelihood-is-differentiable/41541/7 "2026-09-08T20:24:35Z")

</div>

> [@SeHellmann](#):
>
> (Apparently, I should have it in the generated quantities block, but I never understood why this should be faster, because the computations need to be done anyways, but this would be a thing for a different question.)

The speedup of the vectorization in the model block for automatic differentiation more than makes up for the cost of having to compute it again using `double` values in the generated quantities block. The deeper reason is that the reverse pass of reverse-mode autodiff is essentially interpreted and anything we can do to vectorize allows us to use more precompiled code, which is going to be faster for two reasons: (a) less pointer chasing, and (b) shared computations.

For example, if I write this:

```stan
y ~ normal(alpha + beta * x, sigma);

```

and `y` is a long vector and `x` is a long vector, then the result only has to compute `log(sigma)` once. All the intermediate calculations for each of the `x[n], y[n]` pairs then reuse that value and can be calculated with efficient vectorized (in the CPU sense, not the Stan sense) operations.

Then in generated quantities, it’s super fast to execute the loop because there’s no automatic differentiation.

```stan
generated quantities {

  vector[N] log_lik;
  for (n in 1:N) {
    log_lik[n] = normal_lpdf(y[n] | alpha + beta * x[n], sigma);
  }

```

It’s just painful to have to duplicate code like this. We’re working on actually adding efficient vector-return pdfs so that the generated quantities part will be a one-liner.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [September 8, 2026, 8:26pm UTC](https://discourse.mc-stan.org/t/keep-getting-infinite-gradient-although-likelihood-is-differentiable/41541/8 "2026-09-08T20:26:12Z")

</div>

> [@SeHellmann](#):
>
> I am using rstan

If you don’t need something RStan-specific, I’d suggest switching to `cmdstanr` as it should be faster and more up to date with Stan.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [September 8, 2026, 8:28pm UTC](https://discourse.mc-stan.org/t/keep-getting-infinite-gradient-although-likelihood-is-differentiable/41541/9 "2026-09-08T20:28:57Z")

</div>

Sorry, should have put this all in one response.

> [@SeHellmann](#):
>
> I ran a single chain and the memory usage converges towards roughly 200GB

If this is memory during a single iteration, it usually means you have billions of subexpressions introducing their own local autodiff trees. Vectorizing can cut this down dramatically. So can using things like our map functions that do partial evaluations and can reduce dynamic memory overhead.

CmdStan is only going to help if the bottleneck is the draws. That is, something high dimensional that you get a lot of. CmdStan (and presumably cmdstanr) let you just stream the results directly to disk. Processing them afterward can be a problem, though.