Python linalg.pascal函数代码示例

 def test_threshold(self):
     # Regression test.  An early version of `pascal` returned an
     # array of type np.uint64 for n=35, but that data type is too small
     # to hold p[-1, -1].  The second assert_equal below would fail
     # because p[-1, -1] overflowed.
     p = pascal(34)
     assert_equal(2*p.item(-1, -2), p.item(-1, -1), err_msg="n = 34")
     p = pascal(35)
     assert_equal(2*p.item(-1, -2), p.item(-1, -1), err_msg="n = 35")

 def check_invpascal(n, kind, exact):
     ip = invpascal(n, kind=kind, exact=exact)
     p = pascal(n, kind=kind, exact=exact)
     # Matrix-multiply ip and p, and check that we get the identity matrix.
     # We can't use the simple expression e = ip.dot(p), because when
     # n < 35 and exact is True, p.dtype is np.uint64 and ip.dtype is
     # np.int64. The product of those dtypes is np.float64, which loses
     # precision when n is greater than 18.  Instead we'll cast both to
     # object arrays, and then multiply.
     e = ip.astype(object).dot(p.astype(object))
     assert_array_equal(e, eye(n), err_msg="n=%d  kind=%r exact=%r" % (n, kind, exact))

 def check_case(self, n, sym, low):
     assert_array_equal(pascal(n), sym)
     assert_array_equal(pascal(n, kind='lower'), low)
     assert_array_equal(pascal(n, kind='upper'), low.T)
     assert_array_almost_equal(pascal(n, exact=False), sym)
     assert_array_almost_equal(pascal(n, exact=False, kind='lower'), low)
     assert_array_almost_equal(pascal(n, exact=False, kind='upper'), low.T)

 def test_big(self):
     p = pascal(50)
     assert_equal(p[-1, -1], comb(98, 49, exact=True))

    def __init__(self, r, r_min, r_max, c, r_0=0.0, s=1.0, reduced=False):
        # remove zero high-order terms
        c = np.array(np.trim_zeros(c, 'b'), float)
        # if all coefficients are zero
        if len(c) == 0:
            # then both func and abel are also zero everywhere
            self.func = np.zeros_like(r)
            self.abel = self.func
            return
        # polynomial degree
        K = len(c) - 1

        if reduced:
            # rescale r to [0, 1] (to avoid FP overflow)
            r = r / r_max
            r_0 /= r_max
            s /= r_max
            abel_scale = r_max
            r_min /= r_max
            r_max = 1.0

        if s != 1.0:
            # apply stretch
            S = np.cumprod([1.0] + [1.0 / s] * K)  # powers of 1/s
            c *= S
        if r_0 != 0.0:
            # apply shift
            P = pascal(1 + K, 'upper', False)  # binomial coefficients
            rk = np.cumprod([1.0] + [-float(r_0)] * K)  # powers of -r_0
            T = toeplitz([1.0] + [0.0] * K, rk)  # upper-diag. (-r_0)^{l - k}
            c = (P * T).dot(c)

        # whether even and odd powers are present
        even = np.any(c[::2])
        odd = np.any(c[1::2])

        # index limits
        n = r.shape[0]
        i_min = np.searchsorted(r, r_min)
        i_max = np.searchsorted(r, r_max)

        # Calculate all necessary variables within [0, r_max]

        # x, x^2
        x = r[:i_max]
        x2 = x * x

        # integration limits y = sqrt(r^2 - x^2) or 0
        def sqrt0(x): return np.sqrt(x, np.zeros_like(x), where=x > 0)
        y_up = sqrt0(r_max * r_max - x2)
        y_lo = sqrt0(r_min * r_min - x2)

        # y r^k |_lo^up
        # (actually only even are neded for "even", and only odd for "odd")
        Dyr = np.outer(np.cumprod([1.0] + [r_max] * K), y_up) - \
              np.outer(np.cumprod([1.0] + [r_min] * K), y_lo)

        # ln(r + y) |_lo^up, only for odd k
        if odd:
            # ln x for x > 0, otherwise 0
            def ln0(x): return np.log(x, np.zeros_like(x), where=x > 0)
            Dlnry = ln0(r_max + y_up) - \
                    ln0(np.maximum(r_min, x) + y_lo)

        # One-sided Abel integral \int_lo^up r^k dy.
        def a(k):
            odd_k = k % 2
            # max. x power
            K = k - odd_k  # (k - 1 for odd k)
            # generate coefficients for all x^m r^{k - m} terms
            # (only even indices are actually used;
            #  for odd k, C[K] is also used for x^{k+1} ln(r + y))
            C = [0] * (K + 1)
            C[0] = 1 / (k + 1)
            for m in range(k, 1, -2):
                C[k - m + 2] = C[k - m] * m / (m - 1)
            # sum all terms using Horner's method in x
            a = C[K] * Dyr[k - K]
            if odd_k:
                a += C[K] * x2 * Dlnry
            for m in range(K - 2, -1, -2):
                a = a * x2 + C[m] * Dyr[k - m]
            return a

        # Generate the polynomial function
        func = np.zeros(n)
        span = slice(i_min, i_max)
        # (using Horner's method)
        func[span] = c[K]
        for k in range(K - 1, -1, -1):
            func[span] = func[span] * x[span] + c[k]
        self.func = func

        # Generate its Abel transform
        abel = np.zeros(n)
        span = slice(0, i_max)
        if reduced:
            c *= abel_scale
        for k in range(K + 1):
            if c[k]:
#.........这里部分代码省略.........

 def time_pascal(self, size):
     sl.pascal(size)

def locmle(z, xlim = None, Jmle = 35, d = 0., s = 1., ep = 1/100000., sw = 0, Cov_in = None):
	"""Uses z-values in [-xlim,xlim] to find mles for p0, del0, sig0 .

	Jmle is the number of iterations, beginning at (del0, sig0) = (d, s).
	sw = 1 returns the correlation matrix.
	z can be a numpy/scipy array or an ordinary Python array.
	Note that this function returns pandas Series."""
	N = len(z)
	if xlim is None:
		if N > 500000:
			b = 1
		else:
			b = 4.3 * np.exp(-0.26*np.log10(N))
		xlim = np.array([np.median(z), b*(np.percentile(z, 75)-np.percentile(z, 25))/(2*stats.norm.ppf(.75))])
	aorig = xlim[0] - xlim[1]
	borig = xlim[0] + xlim[1]
	z0 = np.array([el for el in z if el >= aorig and el <= borig])
	N0 = len(z0)
	Y = np.array([np.mean(z0), np.mean(np.power(z0, 2))])
	that = float(N0) / N
	# find MLE estimates
	for j in xrange(Jmle):
		bet = np.array([d/(s*s), -1/(2*s*s)])
		aa = (aorig - float(d)) / s
		bb = (borig - float(d)) / s
		H0 = stats.norm.cdf(bb) - stats.norm.cdf(aa)
		fa = stats.norm.pdf(aa)
		fb = stats.norm.pdf(bb)
		H1 = fa - fb
		H2 = H0 + aa * fa - bb * fb
		H3 = (2 + aa*aa) * fa - (2 + bb*bb) * fb
		H4 = 3 * H0 + (3 * aa + np.power(aa, 3)) * fa - (3 * bb + np.power(bb, 3)) * fb
		H = np.array([H0, H1, H2, H3, H4])
		r = float(d) / s
		I = pascal(5, kind = 'lower', exact = False)
		u1hold = np.power(s, range(5))
		u1 = np.matrix([u1hold for k in range(5)])
		II = np.power(r, np.matrix([[max(k-i, 0) for i in range(5)] for k in range(5)]))
		I = np.multiply(np.multiply(I, II), u1.transpose())
		E = np.array(I * np.matrix(H).transpose()).transpose()[0]/H0
		mu = np.array([E[1], E[2]])
		V = np.matrix([[E[2] - E[1]*E[1], E[3] - E[1] * E[2]],[E[3] - E[1] * E[2], E[4] - E[2]*E[2]]])
		addbet = np.linalg.solve(V, (Y - mu).transpose()).transpose()/(1+1./((j+1)*(j+1)))
		bett = bet + addbet
		if bett[1] > 0:
			bett = bet + .1 * addbet
		if pd.isnull(bett[1]) or bett[1] >= 0:
			break
		d = -bett[0]/(2 * bett[1])
		s = 1 / np.sqrt(-2. * bett[1])
		if np.sqrt(sum(np.array(np.power(bett - bet, 2)))) < ep:
			break
	if pd.isnull(bett[1]) or bett[1] >= 0:
		mle = np.array([np.nan for k in xrange(6)])
		Cov_lfdr = np.nan
		if pd.isnull(bett[1]):
			Cov_out = np.nan
		Cor = np.matrix([[np.nan]*3]*3)
	else:
		aa = (aorig - d) / s
		bb = (borig - d) / s
		H0 = stats.norm.cdf(bb) - stats.norm.cdf(aa)
		p0 = that / H0
		# sd calcs
		J = s*s * np.matrix([[1, 2 * d],[0, s]])
		JV = J * np.linalg.inv(V)
		JVJ = JV * J.transpose()
		mat = np.zeros((3,3))
		mat[1:,1:] = JVJ/N0
		mat[0,0] = (p0 * H0 * (1 - p0 * H0)) / N
		h = np.array([H1/H0, (H2 - H0)/H0])
		matt = np.eye(3)
		matt[0,:] = np.array([1/H0] + (-(p0/s) * h).tolist())
		matt = np.matrix(matt)
		C = matt * (mat * matt.transpose())
		mle = np.array([p0, d, s] + np.sqrt(np.diagonal(C)).tolist())
		if sw == 1:
			sd = mle[3:]
			Co = C/np.outer(sd, sd)
			Cor = Co[:,[1,2,0]][[1,2,0]]
			# switch to pandas dataframe for labeling
			Cor = pd.DataFrame(Cor, index=['d', 's','p0'], columns=['d','s','p0'])
		if Cov_in is not None:
			i0 = [i for i,x in enumerate(Cov_in['x']) if x > aa and x < bb]
			Cov_out = loccov(N, N0, p0, d, s, Cov_in['x'], Cov_in['X'], Cov_in['f'], JV, Y, i0, H, h, Cov_in['sw'])
	#label with pandas Series
	mle = pd.Series(mle[[1,2,0,4,5,3]], index=['del0', 'sig0', 'p0', 'sd_del0', 'sd_sig0', 'sd_p0'])
	out = {}
	out['mle'] = mle
	if sw == 1:
		out['Cor'] = Cor
	if Cov_in is not None:
		if Cov_in['sw'] == 2:
			out['pds_'] = Cov_out
		elif Cov_in['sw'] == 3:
			out['Ilfdr'] = Cov_out
		else:
			out['Cov_lfdr'] = Cov_out
	if sw == 1 or Cov_in is not None:
		return pd.Series(out)
#.........这里部分代码省略.........