The CDF, static version.

@param x Value
@param k Shape k
@param theta Theta = 1.0/Beta aka. "scaling" parameter
@return cdf value
512438
public static double cdf(double x, double k, double theta, double shift) {
    if(x == Double.NEGATIVE_INFINITY) {
      return 0.;
    }
    if(x == Double.POSITIVE_INFINITY) {
      return 1.;
    }
    if(x != x) {
      return Double.NaN;
    }
    final double e = FastMath.exp((x - shift) * theta);
    return e < Double.POSITIVE_INFINITY ? GammaDistribution.regularizedGammaP(k, e) : 1.;
  }
473566
def cumulative_distribution(self, X):
        
        self.check_fit()

        U, V = self.split_matrix(X)

        num = np.multiply(
            np.exp(np.multiply(-self.theta, U)) - 1,
            np.exp(np.multiply(-self.theta, V)) - 1
        )
        den = np.exp(-self.theta) - 1

        return -1.0 / self.theta * np.log(1 + num / den)
219374
public static double J0(double x) {
        double ax;

        if ((ax = Math.abs(x)) < 8.0) {
            double y = x * x;
            double ans1 = 57568490574.0 + y * (-13362590354.0 + y * (651619640.7
                    + y * (-11214424.18 + y * (77392.33017 + y * (-184.9052456)))));
            double ans2 = 57568490411.0 + y * (1029532985.0 + y * (9494680.718
                    + y * (59272.64853 + y * (267.8532712 + y * 1.0))));

            return ans1 / ans2;
        } else {
            double z = 8.0 / ax;
            double y = z * z;
            double xx = ax - 0.785398164;
            double ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4
                    + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
            double ans2 = -0.1562499995e-1 + y * (0.1430488765e-3
                    + y * (-0.6911147651e-5 + y * (0.7621095161e-6
                    - y * 0.934935152e-7)));

            return Math.sqrt(0.636619772 / ax) *
                    (Math.cos(xx) * ans1 - z * Math.sin(xx) * ans2);
        }
    }
934353
def fit(self, X):
        
        U, V = self.split_matrix(X)
        self.tau = stats.kendalltau(U, V)[0]
        self.theta = self.compute_theta()
        self.check_theta()
1374483
public static double logpdf(double x, double k, double theta) {
    if(x < 0) {
      return Double.NEGATIVE_INFINITY;
    }
    if(x == 0) {
      return (k == 1.0) ? FastMath.log(theta) : Double.NEGATIVE_INFINITY;
    }
    if(k == 1.0) {
      return FastMath.log(theta) - x * theta;
    }
    final double xt = x * theta;
    return (xt == Double.POSITIVE_INFINITY) ? Double.NEGATIVE_INFINITY : //
        FastMath.log(theta) + (k - 1.0) * FastMath.log(xt) - xt - logGamma(k);
  }
432953
def Theta(i,j,theta,omega_rescaled,omega_min,
			detunings,detuningsij,combinations,detuning_indices,
			Lij,i_d,I_nd,Nnd,
			states=None,verbose=1,other_the=None):
	""""""

	if i==j:
		return ''
	elif j>i:
		raise ValueError,'i should never be less than j.'

	Ne=len(omega_rescaled[0]); Nl=len(theta[0])
	if other_the!=None:
		the=other_the[:]
	else:
		the=[theta[j-1][l]-theta[i-1][l] for l in range(Nl)]

	###########################################################################################
	#This part is about finding detunings that reduce to -omega_ij
	if the==[0 for l in range(Nl)]:
		if omega_rescaled[i-1][j-1]==0:
			return ''
		#We need a combination of detunings that yields -omega_i,j.
		#According to equation labeled "detuning-exception1"
		#-omega_ij= delta^l_ik - delta^l_jk
		#so we seek a k that is lower than i and j
		#and an l such that l is in L_ik and in L_jk

		ii=min(i,j)
		band1=False
		band2=False
		for k in range(1,ii):
			for l in range(1,Nl+1):
				if l in Lij[i-1][k-1] and l in Lij[j-1][k-1]:
					s   ='I found that -omega_'+str(i)+','+str(j)
					s+='= delta^'+str(l)+'_'+str(i)+','+str(k)
					s+= '-delta^'+str(l)+'_'+str(j)+','+str(k)
					s+='\n and also l='+str(l)
					s+=' is in  L'+str(i)+','+str(k)+'='+str(Lij[i-1][k-1])
					s+=' and in  L'+str(j)+','+str(k)+'='+str(Lij[j-1][k-1])+'\n'

					#print s
					band1=True
					break
			if band1:
				break

		if not band1:
			#If no such a k exists then we follow te equation labeled "detuning-exception2"
			#-omega_ij= delta^l_kj - delta^l_ki
			#to look for a k greater than i and j
			#and an l such that l is in L_ik and in L_jk
			for k in range(ii+1,Ne+1):
				for l in range(1,Nl+1):
					if l in Lij[i-1][k-1] and l in Lij[j-1][k-1]:
						s   ='I found that -omega_'+str(i)+','+str(j)
						s+='= delta^'+str(l)+'_'+str(k)+','+str(j)
						s+='-delta^'+str(l)+'_'+str(k)+','+str(i)
						s+='\n and also l='+str(l)
						s+=' is in  L'+str(i)+','+str(k)+'='+str(Lij[i-1][k-1])
						s+=' and in  L'+str(j)+','+str(k)+'='+str(Lij[j-1][k-1])+'\n'
						#print s
						band2=True
						break

				if band2:
					break

		if band1:
			The=''
			#We need to find which detunings are delta^l_i,k and delta^l_j,k.
			#Since they are detunings of laser l, they must have a number greater
			#than those with smaller l:
			acum=sum([detuning_indices[ll] for ll in range(l-1)])
			#We test the indices of all detunings of laser l to find the ones we need.
			for kk in range(detuning_indices[l-1]):
				if detuningsij[l-1][kk][0]+1==I_nd(i) and detuningsij[l-1][kk][1]+1==I_nd(k):
					The+='+detuning('+str(acum+kk+1)+')'
				if detuningsij[l-1][kk][0]+1==I_nd(j) and detuningsij[l-1][kk][1]+1==I_nd(k):
					The+='-detuning('+str(acum+kk+1)+')'
			return The
		elif band2:
			The=''
			#We need to find which detunings are delta^l_k,j and delta^l_k,i.
			#Since they are detunings of laser l, they must have a number greater
			#than those with smaller l:
			acum=sum([detuning_indices[ll] for ll in range(l-1)])
			#We test the indices of all detunings of laser l to find the ones we need.
			for kk in range(detuning_indices[l-1]):
				if detuningsij[l-1][kk][0]+1==I_nd(k) and detuningsij[l-1][kk][1]+1==I_nd(j):
					The+='+detuning('+str(acum+kk+1)+')'
				if detuningsij[l-1][kk][0]+1==I_nd(k) and detuningsij[l-1][kk][1]+1==I_nd(i):
					The+='-detuning('+str(acum+kk+1)+')'
			return The
		else:
			if verbose>1: print 'WARNING: Optical frequencies will be used instead for -omega_',i,j,'=',omega_rescaled[i-1][j-1],'\n'
			return format_double(-omega_rescaled[i-1][j-1])
	###########################################################################################
	#This part is about finding detunings that reduce to (theta_j -theta_i -omega_ij)

	#We establish to which omegaij the detunings should reduce to
	omegaij=[0 for k in range(Nnd)]
	omegaij[I_nd(i)-1]=1;omegaij[I_nd(j)-1]=-1

	#We search for a combination of detunings that reduces to omegaij
	#That is a linear combination of detunings with coefficients the
	#that reduces to omegaij. Here
	detuning_failure=True
	for comb in combinations:
		suma=[sum([the[l]*detunings[l][comb[l]][ii] for l in range(Nl)])  for ii in range(Nnd)]
		if suma==omegaij:
			detuning_failure=False
			break

	#We stop if no combination reduces to the needed expression.
	if detuning_failure:

		if verbose>1: print 'We will see if it is possible to express Theta_'+str(i)+','+str(j)+'=theta_'+str(j)+'-theta_'+str(i)+'-omega_'+str(i)+','+str(j)
		if verbose>1: print 'in terms of other indices a,b such that omega_ab=omega_ij and transition i -> j is allowed by Lij.'

		band3=False
		for a in range(Ne):
			for b in range(a):
				if omega_rescaled[a][b]==omega_rescaled[i-1][j-1] and Lij[a][b]!=[]:
					band3=True
					break
			if band3: break
		if verbose>1: print band3,a,b,i,j

		a=a+1; b=b+1
		if band3 and (i!=a or j!=b):
			if verbose>1: print omega_rescaled[i-1][j-1],omega_rescaled[a-1][b-1]
			if verbose>1: print the
			if verbose>1: print 'This was possible for omega_'+str(a)+','+str(b)
			return Theta(a,b,theta,omega_rescaled,omega_min,
					detunings,detuningsij,combinations,detuning_indices
					,Lij,other_the=the,verbose=verbose,states=states)
		else:
			#verbose=2
			#print 111
			#qi=states[i-1].quantum_numbers
			#qj=states[j-1].quantum_numbers
			#print verbose
			if verbose>0: print 'WARNING: It was impossible to express laser frequencies',the
			if verbose>0: print 'and atomic transition -omega_',i,j,'=',-omega_rescaled[i-1][j-1]
			if verbose>0: print 'in terms of detunings from the transitions given by Lij'
			#for i in range(Ne):
			#	print i+1,Lij[i]
			#print
			if verbose>0: print 'I Will use optical frequencies instead.'
			The=''
			#We give the optical frequencies
			for l in range(Nl):
				a=the[l]
				if a==1:
					The+='+'+format_double(omega_min[l])+'+detuning_knob('+str(l+1)+')'
				elif a==-1:
					The+='-('+format_double(omega_min[l])+'+detuning_knob('+str(l+1)+'))'
				elif a==0:
					The+=''
				elif a>0:
					The+='+'+str(a)+'*'+format_double(omega_min[l])+'+detuning_knob('+str(l+1)+')'
				else:
					The+=    str(a)+'*('+format_double(omega_min[l])+'+detuning_knob('+str(l+1)+'))'

			#We substract omega_ij
			The+=format_double(-omega_rescaled[i-1][j-1])
			if verbose>1: print The
			if verbose>1: print
			return The

	#For each optical frequency in the, we write the corresponding detuning.
	#This way of assigining a global index ll to the detunings ammounts to
	#   ll=   number_of_previous_detunings
	#       + number_of_detuning_ordered_by_row_and_from_left_to_right_column
	The=''
	acum=0
	####################################################################
	#print 'comb',comb,'the',the,'i,j',i,j
	for l in range(Nl):
		if the[l]==1:
			The+='+detuning('+str(acum+comb[l]+1)+')'
		elif the[l]==-1:
			The+='-detuning('+str(acum+comb[l]+1)+')'
		elif the[l]==0:
			pass
		elif the[l]>0:
			The+='+'+str(the[l])+'*detuning('+str(acum+comb[l]+1)+')'
		else:
			The+= str(the[l])+'*detuning('+str(acum+comb[l]+1)+')'
		acum+=detuning_indices[l]

	if The[:1]=='+':
		The=The[1:]

	####################################################################
	#print 'The',The
	return The
1165863
def dlmk(l,m,k,theta1):
    
    
    if m >= k:
        
        factor = sqrt(factorial(l-k)*factorial(l+m)/factorial(l+k)/factorial(l-m))
        part2 = (cos(theta1/2))**(2*l+k-m)*(-sin(theta1/2))**(m-k)/factorial(m-k)
        part3 = sp.hyp2f1(m-l,-k-l,m-k+1,-(tan(theta1/2))**2)

        return factor*part2*part3
    
    else:
        
        return (-1)**(m-k) * dlmk(l,k,m,theta1)
795457
def rdirichlet(theta, size=1):
    
    gammas = np.vstack([rgamma(theta, 1) for i in xrange(size)])
    if size > 1 and np.size(theta) > 1:
        return (gammas.T / gammas.sum(1))[:-1].T
    elif np.size(theta) > 1:
        return (gammas[0] / gammas[0].sum())[:-1]
    else:
        return 1.
545223
public static double cdf(double val, double shape, double location, double scale) {
    if(val < location) {
      return 0;
    }
    val = (val - location) / scale;
    return 1. / (1. + FastMath.pow(val, -shape));
  }
357605
def probability_density(self, X):
        """"""
        self.check_fit()

        U, V = self.split_matrix(X)

        if self.theta == 1:
            return np.multiply(U, V)

        else:
            a = np.power(np.multiply(U, V), -1)
            tmp = np.power(-np.log(U), self.theta) + np.power(-np.log(V), self.theta)
            b = np.power(tmp, -2 + 2.0 / self.theta)
            c = np.power(np.multiply(np.log(U), np.log(V)), self.theta - 1)
            d = 1 + (self.theta - 1) * np.power(tmp, -1.0 / self.theta)
            return self.cumulative_distribution(X) * a * b * c * d
1003681
def v_theta(self):
        
        # compute theta
        norm = np.sqrt(np.sum(self.wxyz**2))
        theta = 2 * np.arccos(self.wxyz[0] / norm)

        # compute the unit vector
        v = np.array(self.wxyz[1:])
        v = v / np.sqrt(np.sum(v**2))

        return v, theta
1173575
public static double logcdf(double val, double shape1, double shape2) {
    if(val == Double.NEGATIVE_INFINITY) {
      return Double.NEGATIVE_INFINITY;
    }
    if(val == Double.POSITIVE_INFINITY) {
      return 0.;
    }
    if(val != val) {
      return Double.NaN;
    }
    if(shape1 == 0.) {
      val = FastMath.exp(-val);
    }
    else {
      double tmp = shape1 * val;
      if(tmp == Double.NEGATIVE_INFINITY) {
        return shape2 == 0 ? 0. : Double.NEGATIVE_INFINITY;
      }
      if(tmp >= 1.) {
        return shape2 == 0 ? Double.NEGATIVE_INFINITY : 0.;
      }
      val = FastMath.exp(FastMath.log1p(-tmp) / shape1);
    }
    if(shape2 == 0.) {
      return -val;
    }
    final double tmp = shape2 * val;
    return tmp < 1. ? FastMath.log1p(-tmp) / shape2 : Double.NEGATIVE_INFINITY;
  }
405450
def get_G(self, k, add_noise=True):
        
        G = np.ones((self.N - 1, self.N - 1))
        if (add_noise):
            noise = pi * 0.1 * np.random.rand(
                (self.N - 1) * (self.N - 1)).reshape((self.N - 1, self.N - 1))
        other_indices = np.delete(range(self.N), k)
        for idx, i in enumerate(other_indices):
            for jdx, j in enumerate(other_indices):
                if (add_noise and
                        i != j):  # do not add noise on diagonal elements.
                    thetak_ij = self.get_inner_angle(k,
                                                     (i, j)) + noise[idx, jdx]
                else:
                    thetak_ij = self.get_inner_angle(k, (i, j))
                G[idx, jdx] = cos(thetak_ij)
                G[jdx, idx] = cos(thetak_ij)
        return G
1257515
public static double cdf(double val, double loc, double scale) {
    val = (val - loc) / scale;
    return 1. / (1. + FastMath.exp(-val));
  }
1043989
def _cdf(self, xloc, dist, length, cache):
        
        output = evaluation.evaluate_forward(
            dist, xloc.reshape(1, -1)).reshape(length, -1)
        assert xloc.shape == output.shape
        return output
635035
public static double pdf(double x, double k, double theta, double shift) {
    x = (x - shift);
    return x <= 0. ? 0. : FastMath.pow(theta, k) / GammaDistribution.gamma(k) * FastMath.pow(1 + x, -(theta + 1.)) * FastMath.pow(FastMath.log1p(x), k - 1.);
  }
1246862
public static double cdf(double val, double k, double theta) {
    if(val < 0) {
      return 0.;
    }
    final double vt = val * theta;
    return (vt == Double.POSITIVE_INFINITY) ? 1. : regularizedGammaP(k, vt);
  }