python - 波动检验的MSS最大似然法

为了提高我的编程技巧，我尝试在下面的Python2.7中实现一个算法（方法5）http://dx.doi.org/10.1016%2FS0076-6879(05)09012-9。可以在这些位置找到实现：显然我不能发布这么多链接如果我的名声提高了，我会把链接贴在这里。

该算法本质上是用于生物学研究，在一定条件下求出细胞的突变率。这是我的尝试，它有错误（注意：我更新了此代码以删除错误，但是仍然没有得到正确的答案）：

    import numpy as np
    import sympy as sp
    from scipy.optimize import minimize


    def leeCoulson(nparray):
            median=np.median(nparray)
            x=sp.Symbol('x')
            M_est=sp.solve(sp.Eq(-x*sp.log(x) - 1.24*x + median,0),x)
            return M_est

    def ctArray(nparray,max):
            list=[0] * int(max+1)
            for i in range(int(max)+1):
                    list[i]=nparray.count(i)
            return list

    values='filename1.csv'
    data=np.genfromtxt(values,delimiter=',')
    mVal=int(max(data))
    ctArray_=ctArray(np.ndarray.tolist(data),mVal)

    ef mssCalc(estM,max=mVal,count=ctArray_):
            def rec(pi,r):
                    pr=(estM/r)+sum([(pi[i]/(r-i+1)) for i in range(0,r-1)])
                    return pr
            prod=1
            pi=[0]*max
            pi[0]=np.exp(-1*estM)
            for r in range(1,max):
                    pi[r]=rec(pi,r)
                    prod=prod*(pi[r]**count[r])
            return -1*prod

    finalM=minimize (mssCalc,leeCoulson(data),method='nelder-mead',options={'xtol':1e-3,'disp':True})
            print finalM

此代码提供以下错误：

        mss-mle_calc.py:37: RuntimeWarning: overflow encountered in multiply
      prod=prod*(pi[r]**count[r])
  /opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/scipy/optimize/optimize.py:462: RuntimeWarning: invalid value encountered in subtract
      numpy.max(numpy.abs(fsim[0] - fsim[1:])) <= ftol):
    Warning: Maximum number of function evaluations has been exceeded.

如果你有时间的话，请帮我把代码改好。

最佳答案

谢谢你的查找，我的代码中有几个愚蠢的错误这是工作代码（据我所知）：

    import numpy as np
    import sympy as sp
    from scipy.optimize import minimize


    def leeCoulson(nparray):
            median=np.median(nparray)
            x=sp.Symbol('x')
            M_est=sp.solve(sp.Eq(-x*sp.log(x) - 1.24*x + median,0),x)
            return float(M_est[0])

    def ctArray(nparray,max):
            list=[0] * int(max+1)
            for i in range(int(max)+1):
                    list[i]=nparray.count(i)
            return list




    values='filename1.csv'
    data=np.genfromtxt(values,delimiter=',')
    mVal=int(max(data))
    ctArray_=ctArray(np.ndarray.tolist(data),mVal)

    def mssCalc(estM,max=mVal,count=ctArray_):
            def rec(pi,r):
                    pr=(estM/r)*sum([(pi[i]/(r-i+1)) for i in range(0,r)])
                    return pr
            prod=1
            pi=[0]*(max+1)
            pi[0]=np.exp(-1.0*estM)
            for r in range(1,max+1):
                    pi[r]=rec(pi,r)
                    prod=prod*(pi[r]**count[r])

            return -1*prod

    finalM=minimize (mssCalc,leeCoulson(data),method='nelder-mead',options={'xtol':1e-3,'disp':True})
    print finalM

关于python - 波动检验的MSS最大似然法，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/39295443/