我正在使用python做一些贝叶斯统计。我已经在python和Fortran 95中进行了编码。Fortran代码的运行速度更快...大约是100倍。我期望Fortran能够更快，但是我真的希望通过使用numpy可以得到python。代码接近，可能在2倍以内。我已经分析了python代码，看起来大部分时间都花在了以下事情上：

scipy.stats.rvs：从分布中随机抽取。我执行了〜19000次，总共耗时3.552秒

numpy.slogdet：计算矩阵行列式的对数。我这样做〜10,000，总共需要2.48秒

numpy.solve：求解线性系统：我将此例程调用约10,000次，总时间为2.557 s

总的来说，我的代码运行时间约为11秒，而我的fortran代码则需要0.092秒。你在开玩笑吗？我真的不是想对python的期望不切实际，而且我当然不希望我的python代码与Fortran一样快。但是要慢100倍。Python必须能够做得更好。以防万一您好奇，这是我的探查器的完整输出：（我不知道为什么它将文本分成几个块）

     1290611 function calls in 11.296 CPU seconds

Ordered by: internal time, function name

ncalls  tottime  percall  cumtime  percall filename:lineno(function)

18973    0.864    0.000    3.552    0.000 /usr/lib64/python2.6/site-packages/scipy/stats/distributions.py:484(rvs)
 9976    0.819    0.000    2.480    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:1559(slogdet)
 9976    0.627    0.000    6.659    0.001 /bluehome/legoses/bce/bayes_GP_integrated_out/python/ce_funcs.py:77(evaluate_posterior)
 9384    0.591    0.000    0.753    0.000 /bluehome/legoses/bce/bayes_GP_integrated_out/python/ce_funcs.py:39(construct_R_matrix)
77852    0.533    0.000    0.533    0.000 :0(array)
37946    0.520    0.000    1.489    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:32(_wrapit)
77851    0.423    0.000    0.956    0.000 /usr/lib64/python2.6/site-packages/numpy/core/numeric.py:216(asarray)
37946    0.360    0.000    0.360    0.000 :0(all)
 9976    0.335    0.000    2.557    0.000 /usr/lib64/python2.6/sitepackages/scipy/linalg/basic.py:23(solve)
107799    0.322    0.000    0.322    0.000 :0(len)

109740    0.301    0.000    0.301    0.000 :0(issubclass)

28357    0.294    0.000    0.294    0.000 :0(prod)
 9976    0.287    0.000    0.957    0.000 /usr/lib64/python2.6/site-packages/scipy/linalg/lapack.py:45(find_best_lapack_type)
    1    0.282    0.282   11.294   11.294 /bluehome/legoses/bce/bayes_GP_integrated_out/python/ce_funcs.py:199(get_rho_lambda_draws)
 9976    0.269    0.000    1.386    0.000 /usr/lib64/python2.6/site-packages/scipy/linalg/lapack.py:60(get_lapack_funcs)
19952    0.263    0.000    0.476    0.000 /usr/lib64/python2.6/site-packages/scipy/linalg/lapack.py:23(cast_to_lapack_prefix)
19952    0.235    0.000    0.669    0.000 /usr/lib64/python2.6/site-packages/numpy/lib/function_base.py:483(asarray_chkfinite)
66833    0.212    0.000    0.212    0.000 :0(log)
18973    0.207    0.000    1.054    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:1427(product)
29931    0.205    0.000    0.205    0.000 :0(reduce)
28949    0.187    0.000    0.856    0.000 :0(map)
 9976    0.175    0.000    0.175    0.000 :0(dot)
47922    0.163    0.000    0.163    0.000 :0(getattr)
 9976    0.157    0.000    0.206    0.000 /usr/lib64/python2.6/site-packages/numpy/lib/twodim_base.py:169(eye)
19952    0.154    0.000    0.271    0.000 /bluehome/legoses/bce/bayes_GP_integrated_out/python/ce_funcs.py:32(loggbeta)
18973    0.151    0.000    0.793    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:1548(all)
19953    0.146    0.000    0.146    0.000 :0(any)
 9976    0.142    0.000    0.316    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:99(_commonType)
 9976    0.133    0.000    0.133    0.000 :0(dgetrf)
18973    0.125    0.000    0.175    0.000 /usr/lib64/python2.6/site-packages/scipy/stats/distributions.py:462(_fix_loc_scale)
39904    0.117    0.000    0.117    0.000 :0(append)
18973    0.105    0.000    0.292    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:1461(alltrue)
19952    0.102    0.000    0.102    0.000 :0(zeros)
19952    0.093    0.000    0.154    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:71(isComplexType)
19952    0.090    0.000    0.090    0.000 :0(split)
 9976    0.089    0.000    2.569    0.000 /bluehome/legoses/bce/bayes_GP_integrated_out/python/ce_funcs.py:62(get_log_determinant_of_matrix)
19952    0.087    0.000    0.134    0.000 /bluehome/legoses/bce/bayes_GP_integrated_out/python/ce_funcs.py:35(logggamma)
 9976    0.083    0.000    0.154    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:139(_fastCopyAndTranspose)
 9976    0.076    0.000    0.125    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:157(_assertSquareness)
 9976    0.074    0.000    0.097    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:151(_assertRank2)
 9976    0.072    0.000    0.119    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:127(_to_native_byte_order)
18973    0.072    0.000    0.072    0.000 /usr/lib64/python2.6/site-packages/scipy/stats/distributions.py:832(_argcheck)
 9976    0.072    0.000    0.228    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:901(diagonal)
 9976    0.070    0.000    0.070    0.000 :0(arange)
 9976    0.061    0.000    0.061    0.000 :0(diagonal)
 9976    0.055    0.000    0.055    0.000 :0(sum)
 9976    0.053    0.000    0.075    0.000 /usr/lib64/python2.6/site-packages/numpy/linalg/linalg.py:84(_realType)
11996    0.050    0.000    0.091    0.000 /usr/lib64/python2.6/site-packages/scipy/stats/distributions.py:1412(_rvs)
 9384    0.047    0.000    0.162    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:1898(prod)
 9976    0.045    0.000    0.045    0.000 :0(sort)
11996    0.041    0.000    0.041    0.000 :0(standard_normal)
 9976    0.037    0.000    0.037    0.000 :0(_fastCopyAndTranspose)
 9976    0.037    0.000    0.037    0.000 :0(hasattr)
 9976    0.037    0.000    0.037    0.000 :0(range)
 6977    0.034    0.000    0.055    0.000 /usr/lib64/python2.6/site-packages/scipy/stats/distributions.py:3731(_rvs)
 9977    0.027    0.000    0.027    0.000 :0(max)
 9976    0.023    0.000    0.023    0.000 /usr/lib64/python2.6/site-packages/numpy/core/numeric.py:498(isfortran)
 9977    0.022    0.000    0.022    0.000 :0(min)
 9976    0.022    0.000    0.022    0.000 :0(get)
 6977    0.021    0.000    0.021    0.000 :0(uniform)
    1    0.001    0.001   11.295   11.295 <string>:1(<module>)
    1    0.001    0.001   11.296   11.296 profile:0(get_rho_lambda_draws(correlations,energies,rho_priors,lambda_e_prior,lambda_z_prior,candidate_sig2_rhos,candidate_sig2_lambda_e,candidate_sig2_lambda_z,3000))
    2    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:445(__call__)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:385(__init__)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:175(_array2string)
    2    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:475(_digits)
    2    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:309(_extendLine)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:317(_formatArray)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:1477(any)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:243(array2string)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/numeric.py:1390(array_str)
    1    0.000    0.000    0.000    0.000 :0(compress)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/arrayprint.py:394(fillFormat)
    6    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/numeric.py:2166(geterr)
   12    0.000    0.000    0.000    0.000 :0(geterrobj)
    0    0.000             0.000          profile:0(profiler)
    1    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/fromnumeric.py:1043(ravel)
    1    0.000    0.000    0.000    0.000 :0(ravel)
    8    0.000    0.000    0.000    0.000 :0(rstrip)
    6    0.000    0.000    0.000    0.000 /usr/lib64/python2.6/site-packages/numpy/core/numeric.py:2070(seterr)
    6    0.000    0.000    0.000    0.000 :0(seterrobj)
    1    0.000    0.000    0.000    0.000 :0(setprofile)

编辑：

这是相关例程的副本

def get_rho_lambda_draws(correlations, energies, rho_priors, lam_e_prior, lam_z_prior,  
                         candidate_sig2_rhos, candidate_sig2_lambda_e, 
                         candidate_sig2_lambda_z, ndraws):

    nBasis = len(correlations[0])
    nStruct = len(correlations)

    rho _draws = [ [0.5 for x in xrange(nBasis)] for y in xrange(ndraws)]
    lambda_e_draws = [ 5 for x in xrange(ndraws)]
    lambda_z_draws = [ 5 for x in xrange(ndraws)]
            
    accept_rhos = array([0. for x in xrange(nBasis)])
    accept_lambda_e = 0.
    accept_lambda_z = 0.

    for i in xrange(1,ndraws):
        if i % 100 == 0:
            print i, "REP<---------------------------------------------------------------------------------"
        #do metropolis to get rho
        rho_draws[i] = [x for x in rho_draws[i-1]]
        lambda_e_draws[i] = lambda_e_draws[i-1]
        lambda_z_draws[i] = lambda_z_draws[i-1]

        rho_vec = [x for x in rho_draws[i-1]]
        R_matrix_before =construct_R_matrix(correlations,correlations,rho_vec)
        post_before = evaluate_posterior(R_matrix_before,rho_vec,energies,lambda_e_draws[i-1],lambda_z_draws[i-1],lam_e_prior,lam_z_prior,rho_priors)

        index = 0
        for j in xrange(nBasis):
            cand = norm.rvs(rho_draws[i-1][j],scale=candidate_sig2_rhos[j])
            if 0.0 < cand < 1.0:
                rho_vec[j] = cand

                R_matrix_after = construct_R_matrix(correlations,correlations,rho_vec)
                post_after = evaluate_posterior(R_matrix_after,rho_vec,energies,lambda_e_draws[i-1],lambda_z_draws[i-1],lam_e_prior,lam_z_prior,rho_priors)
                metrop_value = post_after - post_before
                unif = log(uniform.rvs(0,1))
                if metrop_value > unif:
                    rho_draws[i][j] = cand
                    post_before = post_after
                    accept_rhos[j] += 1
                else:
                    rho_vec[j] = rho_draws[i-1][j]



        R_matrix = construct_R_matrix(correlations,correlations,rho_vec)
        cand = norm.rvs(lambda_e_draws[i-1],scale=candidate_sig2_lambda_e)
        if cand > 0.0:
            post_after = evaluate_posterior(R_matrix,rho_vec,energies,cand,lambda_z_draws[i-1],lam_e_prior,lam_z_prior,rho_priors)

            metrop_value = post_after - post_before
            unif = log(uniform.rvs(0,1))
            if metrop_value > unif:
                lambda_e_draws[i] = cand
                post_before = post_after
                accept_lambda_e = accept_lambda_e + 1


        cand = norm.rvs(lambda_z_draws[i-1],scale=candidate_sig2_lambda_z)
        if cand > 0.0:
            post_after = evaluate_posterior(R_matrix,rho_vec,energies,lambda_e_draws[i],cand,lam_e_prior,lam_z_prior,rho_priors)
            metrop_value = post_after - post_before
            unif = log(uniform.rvs(0,1))
            if metrop_value > unif:
                lambda_z_draws[i] = cand
                post_before = post_after
                accept_lambda_z = accept_lambda_z + 1


    print accept_rhos/ndraws
    print accept_lambda_e/ndraws
    print accept_lambda_z/ndraws
    return [rho_draws,lambda_e_draws,lambda_z_draws]


def evaluate_posterior(R_matrix,rho_vec,energies,lambda_e,lambda_z,lam_e_prior,lam_z_prior,rho_prior_params):

    #    from scipy.linalg import solve
    #from numpy import allclose

    working_matrix = eye(len(R_matrix))/lambda_e + R_matrix/lambda_z
    logdet = get_log_determinant_of_matrix(working_matrix)

    x = solve(working_matrix,energies,sym_pos=True)
    #    if not allclose(dot(working_matrix,x),energies):
#        exit('solve routine didnt work')

    rho_priors = sum([loggbeta(rho_vec[j],rho_prior_params[j][0],rho_prior_params[j][1]) for j in xrange(len(rho_vec))])

    loggposterior = -.5 * logdet - .5*dot(energies,x) + logggamma(lambda_e,lam_e_prior[0],lam_e_prior[1]) + logggamma(lambda_z,lam_z_prior[0],lam_z_prior[1]) + rho_priors #(a_e-1)*log(lambda_e) - b_e*lambda_e + (a_z-1)*log(lambda_z) - b_z*lambda_z + rho_priors
    return loggposterior

def construct_R_matrix(listone,listtwo,rhos):

    return prod(rhos[:]**(4*(listone[:,newaxis]-listtwo)**2),axis=2)

（再次...我不知道为什么在我发帖时将输入分成几部分。希望您可以对它进行解码）

解决方案

很难准确说明代码的状态，但是我怀疑您只是拥有一些（或不可能）非常矢量化的数据。因为显然调用.rvs（）19000次将比.rvs（size = 19000）慢得多。看到：

  In [5]: %timeit x=[scipy.stats.norm().rvs() for i in range(19000)]
  1 loops, best of 3: 1.23 s per loop

  In [6]: %timeit x=scipy.stats.norm().rvs(size=19000)
  1000 loops, best of 3: 1.67 ms per loop

因此，如果您确实有一个不是非常矢量化的代码或算法，那么它应该比fortran慢。