python类erfc()的实例源码

spikesorting.py 文件源码 项目:NeoAnalysis 作者: neoanalysis 项目源码 文件源码 阅读 33 收藏 0 点赞 0 评论 0
def __test_ks(self,x):
        x = x[~np.isnan(x)]
        n = x.size
        x.sort()
        yCDF = np.arange(1,n+1)/float(n)
        notdup = np.hstack([np.diff(x,1),[1]])
        notdup = notdup>0
        x_expcdf = x[notdup]
        y_expcdf = np.hstack([[0],yCDF[notdup]])
        zScores = (x_expcdf-np.mean(x))/np.std(x,ddof=1);
        mu = 0
        sigma = 1
        theocdf = 0.5*erfc(-(zScores-mu)/(np.sqrt(2)*sigma))

        delta1 = y_expcdf[:-1]-theocdf
        delta2 = y_expcdf[1:]-theocdf
        deltacdf = np.abs(np.hstack([delta1,delta2]))
        KSmax = deltacdf.max()
        return KSmax
spikesorting.py 文件源码 项目:NeoAnalysis 作者: neoanalysis 项目源码 文件源码 阅读 29 收藏 0 点赞 0 评论 0
def __test_ks(self,x):
        x = x[~np.isnan(x)]
        n = x.size
        x.sort()
        yCDF = np.arange(1,n+1)/float(n)
        notdup = np.hstack([np.diff(x,1),[1]])
        notdup = notdup>0
        x_expcdf = x[notdup]
        y_expcdf = np.hstack([[0],yCDF[notdup]])
        zScores = (x_expcdf-np.mean(x))/np.std(x,ddof=1);
        mu = 0
        sigma = 1
        theocdf = 0.5*erfc(-(zScores-mu)/(np.sqrt(2)*sigma))

        delta1 = y_expcdf[:-1]-theocdf
        delta2 = y_expcdf[1:]-theocdf
        deltacdf = np.abs(np.hstack([delta1,delta2]))
        KSmax = deltacdf.max()
        return KSmax
nist_sp_800_22_tests.py 文件源码 项目:py-prng 作者: czechnology 项目源码 文件源码 阅读 22 收藏 0 点赞 0 评论 0
def frequency(generator, n_bits, misc=None):
    """Frequency (Monobit) Test.

    Test purpose as described in [NIST10, section 2.1]:
    "The focus of the test is the proportion of zeroes and ones for the entire sequence. The purpose
    of this test is to determine whether the number of ones and zeros in a sequence are
    approximately the same as would be expected for a truly random sequence. The test assesses the
    closeness of the fraction of ones to 1/2, that is, the number of ones and zeroes in a sequence
    should be about the same. All subsequent tests depend on the passing of this test."
    """

    s_n = 0
    for _ in range(n_bits):
        s_n += 2 * generator.random_bit() - 1  # 1 if generator.random_bit() else -1

    s_obs = abs(s_n) / sqrt(n_bits)

    p_value = erfc(s_obs / sqrt(2))

    if type(misc) is dict:
        misc.update(s_n=s_n, s_obs=s_obs)

    return p_value
basic_tests.py 文件源码 项目:py-prng 作者: czechnology 项目源码 文件源码 阅读 21 收藏 0 点赞 0 评论 0
def frequency_test(generator, n_bits, sig_level=None, misc=None, n1=None):
    if n1 is None:
        n0, n1 = _calculate_n0_n1(generator, n_bits)
    else:
        n0 = n_bits - n1

    x1 = ((n0 - n1) ** 2) / n_bits
    p_value = erfc(sqrt(x1 / 2))

    if type(misc) is dict:
        misc.update(n0=n0, n1=n1, p_value=p_value)

    if sig_level is None:
        return x1
    else:
        limit = chi2.ppf(1 - sig_level, 1)
        if type(misc) is dict:
            misc.update(x=x1, limit=limit)
        return x1 <= limit
randtest.py 文件源码 项目:prbg 作者: Lakate 项目源码 文件源码 阅读 29 收藏 0 点赞 0 评论 0
def monobitfrequencytest(binin):
    ''' The focus of the test is the proportion of zeroes and ones for the entire sequence. The purpose of this test is to determine whether that number of ones and zeros in a sequence are approximately the same as would be expected for a truly random sequence. The test assesses the closeness of the fraction of ones to 1/2, that is, the number of ones and zeroes in a sequence should be about the same.'''

    ss = [int(el) for el in binin]
    sc = map(sumi, ss)
    sn = reduce(su, sc)
    sobs = np.abs(sn) / np.sqrt(len(binin))
    pval = spc.erfc(sobs / np.sqrt(2))
    return pval > 0.01
randtest.py 文件源码 项目:prbg 作者: Lakate 项目源码 文件源码 阅读 29 收藏 0 点赞 0 评论 0
def runstest(binin):
    ''' The focus of this test is the total number of zero and one runs in the entire sequence, where a run is an uninterrupted sequence of identical bits. A run of length k means that a run consists of exactly k identical bits and is bounded before and after with a bit of the opposite value. The purpose of the runs test is to determine whether the number of runs of ones and zeros of various lengths is as expected for a random sequence. In particular, this test determines whether the oscillation between such substrings is too fast or too slow.'''
    binin2 = [str(el) for el in binin]
    binin = ''.join(binin2)
    ss = [int(el) for el in binin]
    n = len(binin)
    pi = 1.0 * reduce(su, ss) / n
    vobs = len(binin.replace('0', ' ').split()) + len(binin.replace('1' , ' ').split())
    pval = spc.erfc(abs(vobs-2*n*pi*(1-pi)) / (2 * pi * (1 - pi) * np.sqrt(2*n)))
    return pval > 0.01
randtest.py 文件源码 项目:prbg 作者: Lakate 项目源码 文件源码 阅读 29 收藏 0 点赞 0 评论 0
def spectraltest(binin):
    '''The focus of this test is the peak heights in the discrete Fast Fourier Transform. The purpose of this test is to detect periodic features (i.e., repetitive patterns that are near each other) in the tested sequence that would indicate a deviation from the assumption of randomness. '''

    n = len(binin)
    ss = [int(el) for el in binin]
    sc = map(sumi, ss)
    ft = sff.fft(sc)
    af = abs(ft)[1:n/2+1:]
    t = np.sqrt(np.log(1/0.05)*n)
    n0 = 0.95*n/2
    n1 = len(np.where(af<t)[0])
    d = (n1 - n0)/np.sqrt(n*0.95*0.05/4)
    pval = spc.erfc(abs(d)/np.sqrt(2))
    return pval > 0.01
randtest.py 文件源码 项目:prbg 作者: Lakate 项目源码 文件源码 阅读 28 收藏 0 点赞 0 评论 0
def getfreq(linn, nu):
    val = 0
    for (x, y) in linn:
        if x == nu:
            val = y
    return val


    ''' The focus of this test is the number of times that a particular state occurs in a cumulative sum random walk. The purpose of this test is to detect deviations from the expected number of occurrences of various states in the random walk.'''
    ss = [int(el) for el in binin]
    sc = map(sumi, ss)
    cs = np.cumsum(sc)
    li = []
    for xs in sorted(set(cs)):
        if np.abs(xs) <= 9:
            li.append([xs, len(np.where(cs == xs)[0])])
    j = getfreq(li, 0) + 1
    pval = []
    for xs in xrange(-9, 9 + 1):
        if not xs == 0:
            # pval.append([xs, spc.erfc(np.abs(getfreq(li, xs) - j) / np.sqrt(2 * j * (4 * np.abs(xs) - 2)))])
            pval.append(spc.erfc(np.abs(getfreq(li, xs) - j) / np.sqrt(2 * j * (4 * np.abs(xs) - 2))))
    return pval


    ''' The focus of this test is the frequency of each and every overlapping m-bit pattern. The purpose of the test is to compare the frequency of overlapping blocks of two consecutive/adjacent lengths (m and m+1) against the expected result for a random sequence.'''
    n = len(binin)
    f1a = [(binin + binin[0:m - 1:])[xs:m + xs:] for xs in xrange(n)]
    f1 = [[xs, f1a.count(xs)] for xs in sorted(set(f1a))]
    f2a = [(binin + binin[0:m:])[xs:m + 1 + xs:] for xs in xrange(n)]
    f2 = [[xs, f2a.count(xs)] for xs in sorted(set(f2a))]
    c1 = [1.0 * f1[xs][1] / n for xs in xrange(len(f1))]
    c2 = [1.0 * f2[xs][1] / n for xs in xrange(len(f2))]
    phi1 = reduce(su, map(logo, c1))
    phi2 = reduce(su, map(logo, c2))
    apen = phi1 - phi2
    chisqr = 2.0 * n * (np.log(2) - apen)
    pval = spc.gammaincc(2 ** (m - 1), chisqr / 2.0)
    return pval
lilacsat1.py 文件源码 项目:gr-satellites 作者: daniestevez 项目源码 文件源码 阅读 22 收藏 0 点赞 0 评论 0
def berawgn(EbN0):
    """ Calculates theoretical bit error rate in AWGN (for BPSK and given Eb/N0) """
    return 0.5 * erfc(math.sqrt(10**(float(EbN0)/10)))
ewaldmatrix.py 文件源码 项目:describe 作者: SINGROUP 项目源码 文件源码 阅读 21 收藏 0 点赞 0 评论 0
def _calc_real_matrix(self):
        """
        Determines the Ewald real-space sum.
        """
        fcoords = self.system.relative_pos
        coords = self.system.cartesian_pos
        n_atoms = len(self.system)
        prefactor = 0.5
        ereal = np.zeros((n_atoms, n_atoms), dtype=np.float)

        for i in range(n_atoms):

            # Get points that are within the real space cutoff
            nfcoords, rij, js = self.system.lattice.get_points_in_sphere(
                fcoords,
                coords[i],
                self.rmax,
                zip_results=False
            )
            # Remove the rii term, because a charge does not interact with
            # itself (but does interact with copies of itself).
            mask = rij > 1e-8
            js = js[mask]
            rij = rij[mask]
            nfcoords = nfcoords[mask]

            qi = self.q[i]
            qj = self.q[js]

            erfcval = erfc(self.a * rij)
            new_ereals = erfcval * qi * qj / rij

            # Insert new_ereals
            for k in range(n_atoms):
                ereal[k, i] = np.sum(new_ereals[js == k])

        ereal *= prefactor
        return ereal
general_tests.py 文件源码 项目:PorousMediaLab 作者: biogeochemistry 项目源码 文件源码 阅读 22 收藏 0 点赞 0 评论 0
def transport_equation_test(self):
        '''Check the transport equation integrator'''
        lab = create_lab()
        D = 40
        lab.add_species(True, 'O2', D, 0, bc_top=1, bc_top_type='dirichlet', bc_bot=0, bc_bot_type='neumann')
        lab.dcdt.O2 = '0'
        lab.solve()
        x = np.linspace(0, lab.length, lab.length / lab.dx + 1)
        sol = 1 / 2 * (special.erfc((x - lab.w * lab.tend) / 2 / np.sqrt(D * lab.tend)) + np.exp(
            lab.w * x / D) * special.erfc((x + lab.w * lab.tend) / 2 / np.sqrt(D * lab.tend)))

        assert max(lab.O2.concentration[:, -1] - sol) < 1e-4
analytical_solutions.py 文件源码 项目:PorousMediaLab 作者: biogeochemistry 项目源码 文件源码 阅读 27 收藏 0 点赞 0 评论 0
def transport_equation_boundary_effect():
    '''Check the transport equation integrator'''

    w = 5
    tend = 5
    dx = 0.1
    length = 30
    phi = 1
    dt = 0.001
    lab = Column(length, dx, tend, dt, w)
    D = 5
    lab.add_species(phi, 'O2', D, 0, bc_top=1,
                    bc_top_type='dirichlet', bc_bot=0, bc_bot_type='flux')
    lab.solve()
    x = np.linspace(0, lab.length, lab.length / lab.dx + 1)
    sol = 1 / 2 * (special.erfc((
        x - lab.w * lab.tend) / 2 / np.sqrt(D * lab.tend)) +
        np.exp(lab.w * x / D) * special.erfc((
            x + lab.w * lab.tend) / 2 / np.sqrt(D * lab.tend)))

    plt.figure()
    plt.plot(x, sol, 'k', label='Analytical solution')
    plt.scatter(lab.x[::10], lab.species['O2'].concentration[
                :, -1][::10], marker='x', label='Numerical')
    plt.xlim([x[0], x[-1]])
    ax = plt.gca()
    ax.ticklabel_format(useOffset=False)
    ax.grid(linestyle='-', linewidth=0.2)
    plt.legend()
    plt.tight_layout()
    plt.show()
analytical_solutions.py 文件源码 项目:PorousMediaLab 作者: biogeochemistry 项目源码 文件源码 阅读 22 收藏 0 点赞 0 评论 0
def transport_equation_plot():
    '''Check the transport equation integrator'''

    w = 5
    tend = 5
    dx = 0.1
    length = 100
    phi = 1
    dt = 0.001
    lab = Column(length, dx, tend, dt, w)
    D = 5
    lab.add_species(phi, 'O2', D, 0, bc_top=1,
                    bc_top_type='dirichlet', bc_bot=0, bc_bot_type='dirichlet')
    lab.solve()
    x = np.linspace(0, lab.length, lab.length / lab.dx + 1)
    sol = 1 / 2 * (special.erfc((
        x - lab.w * lab.tend) / 2 / np.sqrt(D * lab.tend)) +
        np.exp(lab.w * x / D) * special.erfc((
            x + lab.w * lab.tend) / 2 / np.sqrt(D * lab.tend)))

    plt.figure()
    plt.plot(x, sol, 'k', label='Analytical solution')
    plt.scatter(lab.x[::10], lab.species['O2'].concentration[
                :, -1][::10], marker='x', label='Numerical')
    plt.xlim([x[0], x[-1]])
    ax = plt.gca()
    ax.ticklabel_format(useOffset=False)
    ax.grid(linestyle='-', linewidth=0.2)
    plt.legend()
    plt.tight_layout()
    plt.show()
heat_transfer_test.py 文件源码 项目:pygfunction 作者: MassimoCimmino 项目源码 文件源码 阅读 23 收藏 0 点赞 0 评论 0
def fls_double(z2, z1, t, dis, alpha, reaSource=True, imgSource=True):
    """ FLS expression for double integral solution.
    """
    r_pos = np.sqrt(dis**2 + (z2 - z1)**2)
    r_neg = np.sqrt(dis**2 + (z2 + z1)**2)
    fls = 0.
    if reaSource:
        fls += 0.5*erfc(r_pos/np.sqrt(4*alpha*t))/r_pos
    if imgSource:
        fls += -0.5*erfc(r_neg/np.sqrt(4*alpha*t))/r_neg
    return fls
utils.py 文件源码 项目:python-psignifit 作者: wichmann-lab 项目源码 文件源码 阅读 23 收藏 0 点赞 0 评论 0
def my_normcdf(x,mu,sigma):
    z = (x-mu) /sigma
    p = .5*erfc(-z/sqrt(2))

    return p
nist_sp_800_22_tests.py 文件源码 项目:py-prng 作者: czechnology 项目源码 文件源码 阅读 26 收藏 0 点赞 0 评论 0
def runs(generator, n_bits, misc=None):
    """Runs Test.

    Test purpose as described in [NIST10, section 2.3]:
    "The focus of this test is the total number of runs in the sequence, where a run is an
    uninterrupted sequence of identical bits. A run of length k consists of exactly k identical bits
    and is bounded before and after with a bit of the opposite value. The purpose of the runs test
    is to determine whether the number of runs of ones and zeros of various lengths is as expected
    for a random sequence. In particular, this test determines whether the oscillation between such
    zeros and ones is too fast or too slow."
    """
    pi = 0
    v_obs = 0
    b0 = None
    for _ in range(n_bits):
        b1 = generator.random_bit()
        pi += b1
        v_obs += b0 != b1
        b0 = b1
    pi /= n_bits

    if type(misc) is dict:
        misc.update(pi=pi, v_obs=v_obs)

    if abs(pi - 1 / 2) >= 2 / sqrt(n_bits):
        return 0

    p_value = erfc(abs(v_obs - 2 * n_bits * pi * (1 - pi)) / (2 * sqrt(2 * n_bits) * pi * (1 - pi)))

    return p_value
nist_sp_800_22_tests.py 文件源码 项目:py-prng 作者: czechnology 项目源码 文件源码 阅读 28 收藏 0 点赞 0 评论 0
def discrete_fourier_transform(generator, n_bits, misc=None):
    """Discrete Fourier Transform (Spectral) Test.

    Test purpose as described in [NIST10, section 2.6]:
    "The focus of this test is the peak heights in the Discrete Fourier Transform of the sequence.
    The purpose of this test is to detect periodic features (i.e., repetitive patterns that are near
    each other) in the tested sequence that would indicate a deviation from the assumption of
    randomness. The intention is to detect whether the number of peaks exceeding the 95 % threshold
    is significantly different than 5 %."
    """
    if n_bits < 1000:
        print("Warning: Sequence should be at least 1000 bits long", file=stderr)

    x = [1 if generator.random_bit() else -1 for _ in range(n_bits)]
    s = fft(x)  # S = DFT(X)
    # print(s)
    m = abs(s)[0:n_bits // 2]  # modulus(S')
    # print("m =", m)
    t = sqrt(log(1 / 0.05) * n_bits)  # the 95% peak height threshold value
    n_0 = 0.95 * n_bits / 2  # expected theoretical number of peaks under t
    n_1 = len([1 for p in m if p < t])
    d = (n_1 - n_0) / sqrt(n_bits * 0.95 * 0.05 / 4)
    p_value = erfc(abs(d) / sqrt(2))

    if type(misc) == dict:
        misc.update(m=m, t=t, n_0=n_0, n_1=n_1, d=d)

    return p_value
nist_sp_800_22_tests.py 文件源码 项目:py-prng 作者: czechnology 项目源码 文件源码 阅读 24 收藏 0 点赞 0 评论 0
def random_excursions_variant(generator, n_bits, misc=None):
    """Random Excursions Variant Test.

    Test purpose as described in [NIST10, section 2.15]:
    "The focus of this test is the total number of times that a particular state is visited (i.e.,
    occurs) in a cumulative sum random walk. The purpose of this test is to detect deviations from
    the expected number of visits to various states in the random walk. This test is actually a
    series of eighteen tests (and conclusions), one test and conclusion for each of the states:
    -9, -8, ..., -1 and +1, +2, ..., +9."
    """

    if n_bits < 10 ** 6:
        print("Warning: Sequence should be at least 10^6 bits long", file=stderr)

    s = [0] * n_bits
    for i in range(n_bits):
        x = 1 if generator.random_bit() else - 1
        s[i] = s[i - 1] + x
    s.append(0)  # leading zero not needed for our implementation

    ksi = [0] * 18
    j = 0
    for x in s:
        if x == 0:
            j += 1
        elif -9 <= x <= 9:
            ksi[x + 9 - (x > 0)] += 1

    if j < 500:
        print("Warning: The number of cycles (zero crossings) should be at least 500 (is %d)" % j,
              file=stderr)

    p_value = list(map(lambda ksi_i, x: erfc(abs(ksi_i - j) / sqrt(2 * j * (4 * abs(x) - 2))),
                       ksi, list(range(-9, 0)) + list(range(1, 10))))

    if type(misc) == dict:
        misc.update(cumsum=s, j=j, ksi=ksi)

    return p_value
rfunc.py 文件源码 项目:pastas 作者: pastas 项目源码 文件源码 阅读 26 收藏 0 点赞 0 评论 0
def polder_function(self, x, y):
        s = .5 * np.exp(2 * x) * erfc(x / y + y) + \
            .5 * np.exp(-2 * x) * erfc(x / y - y)
        return s
analytical.py 文件源码 项目:SimplePetro 作者: ishovkun 项目源码 文件源码 阅读 21 收藏 0 点赞 0 评论 0
def analytical_rate(td):
    '''
    Returns the dimensionless production rate
    Input:
        td: np.array(n)
            Array of dimensionless times
    Returns:
        qd: np.array(n)
            Array of dimensionless production rates
    '''
    pi = np.pi
    term1 = 2*np.exp(-pi**2*td/4.)
    term2 = erfc(1.5*pi*np.sqrt(td))/np.sqrt(pi*td)
    return term1 + term2
PsiMarginal.py 文件源码 项目:Psi-staircase 作者: NNiehof 项目源码 文件源码 阅读 27 收藏 0 点赞 0 评论 0
def pf(parameters, psyfun='cGauss'):
    """Generate conditional probabilities from psychometric function.

    Arguments
    ---------
        parameters: ndarray (float64) containing parameters as columns
            mu   : threshold

            sigma    : slope

            gamma   : guessing rate (optional), default is 0.2

            lambda  : lapse rate (optional), default is 0.04

            x       : stimulus intensity

        psyfun  : type of psychometric function.
                'cGauss' cumulative Gaussian

                'Gumbel' Gumbel, aka log Weibull

    Returns
    -------
    1D-array of conditional probabilities p(response | mu,sigma,gamma,lambda,x)
    """

    # Unpack parameters
    if np.size(parameters, 1) == 5:
        [mu, sigma, gamma, llambda, x] = np.transpose(parameters)
    elif np.size(parameters, 1) == 4:
        [mu, sigma, llambda, x] = np.transpose(parameters)
        gamma = llambda
    elif np.size(parameters, 1) == 3:
        [mu, sigma, x] = np.transpose(parameters)
        gamma = 0.2
        llambda = 0.04
    else:  # insufficient number of parameters will give a flat line
        psyfun = None
        gamma = 0.2
        llambda = 0.04
    # Psychometric function
    ones = np.ones(np.shape(mu))
    if psyfun == 'cGauss':
        # F(x; mu, sigma) = Normcdf(mu, sigma) = 1/2 * erfc(-sigma * (x-mu) /sqrt(2))
        z = np.divide(np.subtract(x, mu), sigma)
        p = 0.5 * erfc(-z / np.sqrt(2))
    elif psyfun == 'Gumbel':
        # F(x; mu, sigma) = 1 - exp(-10^(sigma(x-mu)))
        p = ones - np.exp(-np.power((np.multiply(ones, 10.0)), (np.multiply(sigma, (np.subtract(x, mu))))))
    elif psyfun == 'Weibull':
        # F(x; mu, sigma)
        p = 1 - np.exp(-(np.divide(x, mu)) ** sigma)
    else:
        # flat line if no psychometric function is specified
        p = np.ones(np.shape(mu))
    y = gamma + np.multiply((ones - gamma - llambda), p)
    return y


问题


面经


文章

微信
公众号

扫码关注公众号