scipy · tupui · Aug 27, 2021 · May 23, 2020 · May 23, 2020 · May 23, 2020
diff --git a/scipy/stats/_distn_infrastructure.py b/scipy/stats/_distn_infrastructure.py
@@ -1240,8 +1240,13 @@ def moment(self, n, *args, **kwds):
 
         """
         args, loc, scale = self._parse_args(*args, **kwds)
-        if not (self._argcheck(*args) and (scale > 0)):
-            return nan
+        arrs = tuple(map(np.asarray, [*args, loc, scale]))
+        arrs = np.broadcast_arrays(*arrs)
+        args = arrs[:-2]
+        loc, scale = arrs[-2], arrs[-1]
+        cond0 = (self._argcheck(*args) and scale > 0)
+        cond1 = (loc != 0)
+        output = zeros(loc.shape, dtype='d')
         if (floor(n) != n):
             raise ValueError("Moment must be an integer.")
         if (n < 0):
@@ -1257,16 +1262,18 @@ def moment(self, n, *args, **kwds):
 
         # Convert to transformed  X = L + S*Y
         # E[X^n] = E[(L+S*Y)^n] = L^n sum(comb(n, k)*(S/L)^k E[Y^k], k=0...n)
-        if loc == 0:
-            return scale**n * val
-        else:
-            result = 0
-            fac = float(scale) / float(loc)
-            for k in range(n):
-                valk = _moment_from_stats(k, mu, mu2, g1, g2, self._munp, args)
-                result += comb(n, k, exact=True)*(fac**k) * valk
-            result += fac**n * val
-            return result * loc**n
+        fac = _lazywhere(cond1, (scale, loc),
+                         lambda scale, loc: scale / loc, self.badvalue)
+        for k in range(n):
+            valk = _moment_from_stats(k, mu, mu2, g1, g2, self._munp, args)
+            output = output + comb(n, k, exact=True)*fac**k * valk
+        output = output + fac**n * val
+        output = output * loc**n
+        output = np.where(cond1, output, scale**n * val)
+        output = np.where(cond0, output, self.badvalue)
+        if output.ndim == 0:
+            return output.item()
+        return output
 
     def median(self, *args, **kwds):
         """