code: def intervalimpl(
      n,
      sumw,
      sumwsq,
      sumwr,
      sumwsqr,
      sumwsqrsq,
      wmin,
      wmax,
      alpha=0.05,
      rmin=0,
      rmax=1,
      raiseonerr=False,
  ):
      from math import inf, isclose, sqrt
      from scipy.stats import f
  
      assert wmin < 1
      assert wmax > 1
      assert rmin <= rmax
  
      uncwfake = wmax if sumw < n else wmin
      if uncwfake == inf:
          uncgstar = 1 + 1 / n
      else:
          unca = (uncwfake + sumw) / (1 + n)
          uncb = (uncwfake**2 + sumwsq) / (1 + n)
          uncgstar = (n + 1) * (unca - 1) ** 2 / (uncb - unca * unca)
      Delta = f.isf(q=alpha, dfn=1, dfd=n)
      phi = (-uncgstar - Delta) / (2 * (n + 1))
  
      bounds = []
      for r, sign in ((rmin, 1), (rmax, -1)):
          candidates = []
          for wfake in (wmin, wmax):
              if wfake == inf:
                  x = sign * (r + (sumwr - sumw * r) / n)
                  y = (r * sumw - sumwr) ** 2 / (n * (1 + n)) - (
                      r**2 * sumwsq - 2 * r * sumwsqr + sumwsqrsq
                  ) / (1 + n)
                  z = phi + 1 / (2 * n)
                  if isclose(y * z, 0, abs_tol=1e-9):
                      y = 0
  
                  if z <= 0 and y * z >= 0:
                      kappa = sqrt(y / (2 * z))
                      if isclose(kappa, 0):
                          candidates.append((sign * r, None))
                      else:
                          gstar = x - sqrt(2 * y * z)
                          gamma = (
                              -kappa * (1 + n) / n
                              + sign * (r * sumw - sumwr) / n
                          )
                          beta = -sign * r
                          candidates.append(
                              (
                                  gstar,
                                  {
                                      "kappastar": kappa,
                                      "betastar": beta,
                                      "gammastar": gamma,
                                      "wfake": wfake,
                                      # Q_{w,r} &= -\frac{\gamma + \beta w + w r}{(N+1) \kappa} \\
                                      "qfunc": lambda c, w, r, k=kappa, g=gamma, b=beta, s=sign, num=n: -c
                                      * (g + (b + s * r) * w)
                                      / ((num + 1) * k),
                                  },
                              )
                          )
              else:
                  barw = (wfake + sumw) / (1 + n)
                  barwsq = (wfake * wfake + sumwsq) / (1 + n)
                  barwr = sign * (wfake * r + sumwr) / (1 + n)
                  barwsqr = (
                      sign * (wfake * wfake * r + sumwsqr) / (1 + n)
                  )
                  barwsqrsq = (wfake * wfake * r * r + sumwsqrsq) / (
                      1 + n
                  )
  
                  if barwsq > barw**2:
                      x = barwr + (
                          (1 - barw)
                          * (barwsqr - barw * barwr)
                          / (barwsq - barw**2)
                      )
                      y = (barwsqr - barw * barwr) ** 2 / (
                          barwsq - barw**2
                      ) - (barwsqrsq - barwr**2)
                      z = phi + (1 / 2) * (1 - barw) ** 2 / (
                          barwsq - barw**2
                      )
                      if isclose(y * z, 0, abs_tol=1e-9):
                          y = 0
  
                      if z <= 0 and y * z >= 0:
                          kappa = sqrt(y / (2 * z)) if y * z > 0 else 0
                          if isclose(kappa, 0):
                              candidates.append((sign * r, None))
                          else:
                              gstar = x - sqrt(2 * y * z)
                              beta = (
                                  -kappa * (1 - barw)
                                  - (barwsqr - barw * barwr)
                              ) / (barwsq - barw * barw)
                              gamma = -kappa - beta * barw - barwr
                              candidates.append(
                                  (
                                      gstar,
                                      {
                                          "kappastar": kappa,
                                          "betastar": beta,
                                          "gammastar": gamma,
                                          "wfake": wfake,
                                          # Q_{w,r} &= -\frac{\gamma + \beta w + w r}{(N+1) \kappa} \\
                                          "qfunc": lambda c, w, r, k=kappa, g=gamma, b=beta, s=sign, num=n: -c
                                          * (g + (b + s * r) * w)
                                          / ((num + 1) * k),
                                      },
                                  )
                              )
  
          best = min(candidates, key=lambda x: x[0])
          vbound = min(rmax, max(rmin, sign * best[0]))
          bounds.append((vbound, best[1]))
  
      return (bounds[0][0], bounds[1][0]), (bounds[0][1], bounds[1][1])
  
label_desc: Use image node_id 5 for calling the CrMinusTwo obj's underlying member method code with example usage: obj.intervalimpl(n, sumw, sumwsq, sumwr, sumwsqr, sumwsqrsq, wmin, wmax, alpha, rmin, rmax, raiseonerr) and returns: bounds,  bounds,  bounds,  bounds