diff --git a/pytensor/scalar/math.py b/pytensor/scalar/math.py
index 8dacc61f6f..22ec00defa 100644
--- a/pytensor/scalar/math.py
+++ b/pytensor/scalar/math.py
@@ -1357,13 +1357,13 @@ def c_code(self, *args, **kwargs):
 
 
 def betainc_grad(p, q, x, wrtp: bool):
-    """Gradient of the regularized lower gamma function (P) wrt to the first
-    argument (k, a.k.a. alpha).
-
-    Adapted from STAN `grad_reg_lower_inc_gamma.hpp`
+    """
+    Gradient of the regularized incomplete beta function wrt to the first
+    argument `p` (aka alpha) or the second argument `q` (aka beta),
+    depending on whether `wrtp` is true.
 
-    Reference: Gautschi, W. (1979). A computational procedure for incomplete gamma functions.
-    ACM Transactions on Mathematical Software (TOMS), 5(4), 466-481.
+    Reference: Boik, R. J., & Robison-Cox, J. F. (1998). Derivatives of the incomplete beta function.
+    Journal of Statistical Software, 3(1), 1-20.
     """
 
     def _betainc_der(p, q, x, wrtp, skip_loop):