Fix CodeRabbit review findings

- Guard gradient computations against theta=0 division by zero
  in geometric_lpmf, geometric_cdf, and geometric_lcdf
- Fix include ordering in prob.hpp (gaussian before geometric)
- Fix signed/unsigned comparison in geometric_rng loop index

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>

Author: yasumorishima

stan/math/prim/prob.hpp

@@ -111,13 +111,13 @@
 #include <stan/math/prim/prob/gamma_lcdf.hpp>
 #include <stan/math/prim/prob/gamma_lpdf.hpp>
 #include <stan/math/prim/prob/gamma_rng.hpp>
+#include <stan/math/prim/prob/gaussian_dlm_obs_lpdf.hpp>
+#include <stan/math/prim/prob/gaussian_dlm_obs_rng.hpp>
 #include <stan/math/prim/prob/geometric_cdf.hpp>
 #include <stan/math/prim/prob/geometric_lccdf.hpp>
 #include <stan/math/prim/prob/geometric_lcdf.hpp>
 #include <stan/math/prim/prob/geometric_lpmf.hpp>
 #include <stan/math/prim/prob/geometric_rng.hpp>
-#include <stan/math/prim/prob/gaussian_dlm_obs_lpdf.hpp>
-#include <stan/math/prim/prob/gaussian_dlm_obs_rng.hpp>
 #include <stan/math/prim/prob/gumbel_ccdf_log.hpp>
 #include <stan/math/prim/prob/gumbel_cdf.hpp>
 #include <stan/math/prim/prob/gumbel_cdf_log.hpp>

stan/math/prim/prob/geometric_cdf.hpp

@@ -75,8 +75,10 @@ inline return_type_t<T_prob> geometric_cdf(const T_n& n,
     cdf *= cdf_i;
 
     if constexpr (is_autodiff_v<T_prob>) {
-      partials<0>(ops_partials)[i]
-          += np1 * ccdf_i / ((1.0 - theta_val) * cdf_i);
+      if (cdf_i > 0.0 && theta_val > 0.0) {
+        partials<0>(ops_partials)[i]
+            += np1 * ccdf_i / ((1.0 - theta_val) * cdf_i);
+      }
     }
   }
 

stan/math/prim/prob/geometric_lcdf.hpp

@@ -75,8 +75,10 @@ inline return_type_t<T_prob> geometric_lcdf(const T_n& n,
 
     if constexpr (is_autodiff_v<T_prob>) {
       const auto ccdf = stan::math::exp(log_ccdf);
-      partials<0>(ops_partials)[i]
-          += np1 * ccdf / ((1.0 - theta_val) * (1.0 - ccdf));
+      if (theta_val > 0.0 && ccdf < 1.0) {
+        partials<0>(ops_partials)[i]
+            += np1 * ccdf / ((1.0 - theta_val) * (1.0 - ccdf));
+      }
     }
   }
 

stan/math/prim/prob/geometric_lpmf.hpp

@@ -78,6 +78,11 @@ inline return_type_t<T_prob> geometric_lpmf(const T_n& n,
     const auto theta_val = theta_vec.val(i);
     const auto n_val = n_vec.val(i);
 
+    // When theta == 0.0, P(n) = 0 for all n
+    if (theta_val == 0.0) {
+      return negative_infinity();
+    }
+
     // When theta == 1.0, P(n=0) = 1, P(n>0) = 0
     if (theta_val == 1.0) {
       if (n_val > 0) {

stan/math/prim/prob/geometric_rng.hpp

@@ -58,7 +58,7 @@ inline auto geometric_rng(T_prob&& theta, RNG& rng) {
   } else {
     auto theta_arr = as_array_or_scalar(theta_ref);
     std::vector<int> result(theta_arr.size());
-    for (int i = 0; i < theta_arr.size(); i++) {
+    for (size_t i = 0; i < theta_arr.size(); i++) {
       if (theta_arr[i] == 1.0) {
         result[i] = 0;
       } else {