Convergence Troubles and Pathologies in an already Reparameterized Model

Modeling
, , , ,
6

That code looks a bit off. Before you had

  1 - exp(-(max_dist[i,k]^2 / exp(log_tau[species[i]])^2))

but now you have

  1 - exp(-(max_dist[i,k]^2 / exp(log_tau[species[i]]^2)))

and squaring log_tau does not really make sense–it’s tau that’s squared in the model. Maybe you were thinking of

  1 - exp(-(max_dist[i,k]^2 / exp(log_tau[species[i]]*2)))

Simplifying the expressions should speed things up a bit. This should be equivalent to the first function but with fewer operations.

  real partial_sum_lpmf(array[, ] int slice_abund_per_band,
                        int start, int end,
                        int max_intervals,
                        array[] int bands_per_sample,
                        array[, ] real max_dist,
                        array[] int species,
                        vector log_tau)
  {
    int Pi_size = end - start + 1;
    int Pi_index = 1;
    array[Pi_size] real lps;
    for (i in start:end)
    {
      int n_max = bands_per_sample[i];
      real neg_inv_tau = -exp(-2*log_tau[species[i]]);
      vector[n_max] Ri = exp(neg_inv_tau * square(to_vector(max_dist[i])));
      real C = inv(1 - Ri[n_max]);
      vector[n_max] Pi;
      Pi[1] = C*(1 - Ri[1]);
      for (k in 2 : n_max) {
        Pi[k] = C*(Ri[k-1] - Ri[k]);
      }
      lps[Pi_index] = multinomial_lupmf(slice_abund_per_band[Pi_index, 1:n_max] | Pi);
      Pi_index = Pi_index + 1;
    }
    
    return sum(lps);
  }

I changed multinomial_lpmf to multinomial_lupmf because the normalizing constant is unnecessary but could be expensive to compute. To take advantage of that _lupmf form must be used with reduce_sum() as well.

  target += reduce_sum(partial_sum_lupmf,
                       abund_per_band,
                       grainsize,
                       max_intervals,
                       bands_per_sample,
                       max_dist,
                       species,
                       log_tau);
1 Like