Returns the regularized incomplete beta function I_x(a, b) Includes the
continued fraction way of computing, based on the book "Numerical Recipes".

@param alpha Parameter a
@param beta Parameter b
@param x Parameter x
@return result
349443
protected static double regularizedIncBetaCF(double alpha, double beta, double x) {
    final double FPMIN = Double.MIN_VALUE / NUM_PRECISION;
    double qab = alpha + beta;
    double qap = alpha + 1.0;
    double qam = alpha - 1.0;
    double c = 1.0;
    double d = 1.0 - qab * x / qap;
    if(Math.abs(d) < FPMIN) {
      d = FPMIN;
    }
    d = 1.0 / d;
    double h = d;
    for(int m = 1; m < 10000; m++) {
      int m2 = 2 * m;
      double aa = m * (beta - m) * x / ((qam + m2) * (alpha + m2));
      d = 1.0 + aa * d;
      if(Math.abs(d) < FPMIN) {
        d = FPMIN;
      }
      c = 1.0 + aa / c;
      if(Math.abs(c) < FPMIN) {
        c = FPMIN;
      }
      d = 1.0 / d;
      h *= d * c;
      aa = -(alpha + m) * (qab + m) * x / ((alpha + m2) * (qap + m2));
      d = 1.0 + aa * d;
      if(Math.abs(d) < FPMIN) {
        d = FPMIN;
      }
      c = 1.0 + aa / c;
      if(Math.abs(c) < FPMIN) {
        c = FPMIN;
      }
      d = 1.0 / d;
      double del = d * c;
      h *= del;
      if(Math.abs(del - 1.0) <= NUM_PRECISION) {
        break;
      }
    }
    return h;
  }
673946
def up_alpha_beta(returns, factor_returns, **kwargs):
    
    return up(returns, factor_returns, function=alpha_beta_aligned, **kwargs)
487231
public static double regularizedUpperIncompleteGamma(double s, double x) {
        if (s < 0.0) {
            throw new IllegalArgumentException("Invalid s: " + s);
        }

        if (x < 0.0) {
            throw new IllegalArgumentException("Invalid x: " + x);
        }

        double igf = 0.0;

        if (x != 0.0) {
            if (x == 1.0 / 0.0) {
                igf = 1.0;
            } else {
                if (x < s + 1.0) {
                    // Series representation
                    igf = 1.0 - regularizedIncompleteGammaSeries(s, x);
                } else {
                    // Continued fraction representation
                    igf = 1.0 - regularizedIncompleteGammaFraction(s, x);
                }
            }
        }
        return igf;
    }
218968
public static double differenceBetweenAngles(double alpha, double beta) {
    double phi = Math.abs(beta - alpha) % 360;
    return phi > 180 ? 360 - phi : phi;
  }
523198
public static void jacobian_Control4( DMatrixRMaj L_full ,
										  double beta[] , DMatrixRMaj A )
	{
		int indexA = 0;

		double b0 = beta[0]; double b1 = beta[1]; double b2 = beta[2]; double b3 = beta[3];

		final double ld[] = L_full.data;

		for( int i = 0; i < 6; i++ ) {
			int li = L_full.numCols*i;
			A.data[indexA++] = 2*ld[li+0]*b0 +   ld[li+1]*b1 +   ld[li+2]*b2 +   ld[li+3]*b3;
			A.data[indexA++] =   ld[li+1]*b0 + 2*ld[li+4]*b1 +   ld[li+5]*b2 +   ld[li+6]*b3;
			A.data[indexA++] =   ld[li+2]*b0 +   ld[li+5]*b1 + 2*ld[li+7]*b2 +   ld[li+8]*b3;
			A.data[indexA++] =   ld[li+3]*b0 +   ld[li+6]*b1 +   ld[li+8]*b2 + 2*ld[li+9]*b3;
		}
	}
352832
@Override
    public void scal(long N, double alpha, INDArray X) {
        if (Nd4j.getExecutioner().getProfilingMode() == OpExecutioner.ProfilingMode.ALL)
            OpProfiler.getInstance().processBlasCall(false, X);

        if (X.isSparse()) {
            Nd4j.getSparseBlasWrapper().level1().scal(N, alpha, X);
        } else if (X.data().dataType() == DataType.DOUBLE)
            dscal(N, alpha, X, BlasBufferUtil.getBlasStride(X));
        else if (X.data().dataType() == DataType.FLOAT)
            sscal(N, (float) alpha, X, BlasBufferUtil.getBlasStride(X));
        else if (X.data().dataType() == DataType.HALF)
            Nd4j.getExecutioner().exec(new ScalarMultiplication(X, alpha));
    }
833865
public static DoubleMatrix gemm(double alpha, DoubleMatrix a,
			DoubleMatrix b, double beta, DoubleMatrix c) {
		NativeBlas.dgemm('N', 'N', c.rows, c.columns, a.columns, alpha, a.data, 0,
				a.rows, b.data, 0, b.rows, beta, c.data, 0, c.rows);
		return c;
	}
1148108
public final Result solve(GradientSolver gslvr, double [] beta, GradientInfo ginfo, ProgressMonitor pm) {
    if(_hist == null)
      _hist = new History(_historySz, beta.length);
    int iter = 0;
    double rel_improvement = 1;
    final double [] pk = new double[beta.length];
    double minStep = 1e-16;
    LineSearchSolver lineSearch = new MoreThuente(gslvr,beta,ginfo);
    while(!ArrayUtils.hasNaNsOrInfs(beta) && (ArrayUtils.linfnorm(ginfo._gradient,false) > _gradEps  && rel_improvement > _objEps) && iter != _maxIter) {
      ++iter;
      _hist.getSearchDirection(ginfo._gradient,pk);
      if(!lineSearch.evaluate(pk))
        break;
      lineSearch.setInitialStep(Math.max(minStep, lineSearch.step()));
      GradientInfo newGinfo = lineSearch.ginfo();
      _hist.update(pk, newGinfo._gradient, ginfo._gradient);
      rel_improvement = (ginfo._objVal - newGinfo._objVal)/ginfo._objVal;
      ginfo = newGinfo;
      if(!pm.progress(lineSearch.getX(), ginfo))break;
    }
    return new Result((ArrayUtils.linfnorm(ginfo._gradient,false) <= _gradEps  || rel_improvement <= _objEps),iter,lineSearch.getX(), lineSearch.ginfo(),rel_improvement);
  }
1066480
def Betainc(a, b, x):
    
    return sp.special.betainc(a, b, x),
652366
def _B(self, x, a, b):
        
        return special.betainc(a, b, x) * special.beta(a, b)
552652
def MarginalBeta(self, i):
        
        alpha0 = self.params.sum()
        alpha = self.params[i]
        return Beta(alpha, alpha0 - alpha)
1020201
def down_alpha_beta(returns, factor_returns, **kwargs):
    
    return down(returns, factor_returns, function=alpha_beta_aligned, **kwargs)
368104
private static double incompleteFractionSummation(double alpha, double beta, double x) {
        final int MAXITER = 500;
        final double EPS = 3.0E-7;

        double aplusb = alpha + beta;
        double aplus1 = alpha + 1.0;
        double aminus1 = alpha - 1.0;
        double c = 1.0;
        double d = 1.0 - aplusb * x / aplus1;
        if (Math.abs(d) < FPMIN) {
            d = FPMIN;
        }
        d = 1.0 / d;
        double h = d;
        double aa = 0.0;
        double del = 0.0;
        int i = 1, i2 = 0;
        boolean test = true;
        while (test) {
            i2 = 2 * i;
            aa = i * (beta - i) * x / ((aminus1 + i2) * (alpha + i2));
            d = 1.0 + aa * d;
            if (Math.abs(d) < FPMIN) {
                d = FPMIN;
            }
            c = 1.0 + aa / c;
            if (Math.abs(c) < FPMIN) {
                c = FPMIN;
            }
            d = 1.0 / d;
            h *= d * c;
            aa = -(alpha + i) * (aplusb + i) * x / ((alpha + i2) * (aplus1 + i2));
            d = 1.0 + aa * d;
            if (Math.abs(d) < FPMIN) {
                d = FPMIN;
            }
            c = 1.0 + aa / c;
            if (Math.abs(c) < FPMIN) {
                c = FPMIN;
            }
            d = 1.0 / d;
            del = d * c;
            h *= del;
            i++;
            if (Math.abs(del - 1.0) < EPS) {
                test = false;
            }
            if (i > MAXITER) {
                test = false;
                logger.error("Beta.incompleteFractionSummation: Maximum number of iterations wes exceeded");
            }
        }
        return h;
    }
944589
public static double pdf(double val, double alpha, double beta) {
    return (alpha <= 0. || beta <= 0. || Double.isNaN(alpha) || Double.isNaN(beta) || Double.isNaN(val)) ? Double.NaN : //
        (val < 0. || val > 1.) ? 0. : //
            val == 0. ? (alpha > 1. ? 0. : alpha < 1. ? Double.POSITIVE_INFINITY : beta) : //
                val == 1. ? (beta > 1. ? 0. : beta < 1. ? Double.POSITIVE_INFINITY : alpha) : //
                    FastMath.exp(-logBeta(alpha, beta) + FastMath.log(val) * (alpha - 1) + FastMath.log1p(-val) * (beta - 1));
  }
539281
function bfrac(a, b, x, y) {
      var lambda;

      lambda = a - (a + b) * x;
      return brcomp(a, b, x, y) + Math.log(contFrac(
         function(i) { // betas from DiDonato
            var n;

            n = i - 1;
            return i === 0 ? 0
                           : n + n * (b - n) * x / (a + 2 * n - 1) +
                             (a + n) / (a + 2 * n + 1) *
                                (lambda + 1 + n * (1 + y));
         },
         function(i) { // alphas from DiDonato
            var n;

            n = i - 1;
            return i === 1 ? 1
                           : (a + n - 1) * (a + b + n - 1) /
                             (a + 2 * n - 1) / (a + 2 * n - 1) *
                             n * (b - n) * x * x;
         }
      ));
   }
485434
def function(self, x, y, amp, alpha, beta, center_x, center_y):
        
        x_shift = x - center_x
        y_shift = y - center_y
        return amp * (1. + (x_shift**2+y_shift**2)/alpha**2)**(-beta)
1312826
def alpha(self, theta_x, theta_y, kwargs_lens, k=None):
        
        beta_x, beta_y = self.ray_shooting(theta_x, theta_y, kwargs_lens)
        alpha_x = theta_x - beta_x
        alpha_y = theta_y - beta_y
        return alpha_x, alpha_y