The MSS-maximum likelihood method for Fluctuation Test
<p>I am trying to implement an algorithm in the following paper (method 5) <a href="http://dx.doi.org/10.1016%2FS0076-6879(05)09012-9" rel="nofollow">http://dx.doi.org/10.1016%2FS0076-6879(05)09012-9</a> in python2.7 to improve my programming skill. Implementations can be found at these locations: Apparently I cannot post so many links. If my reputation goes up, I will post the links here.</p>

<p><a href="http://i.stack.imgur.com/RsSYZ.png" rel="nofollow"><img src="http://i.stack.imgur.com/RsSYZ.png" alt="This is the function I am trying to implement"></a></p>

<p>Essentially, the algorithm is used for biological research, and finds the mutation rate of cells under some condition. Here is my attempt, which has errors (NOTE: I updated this code to remove an error however I am still not getting the right answer):</p>

<pre><code>    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
</code></pre>

<p>This code gives the following errors:</p>

<pre><code>        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:])) &lt;= ftol):
    Warning: Maximum number of function evaluations has been exceeded.
</code></pre>

<p>Please help me make this code better if you have some time.</p>

<p>Thanks for looking, there were a couple stupid mistakes in my code. Here is the working code (as far as I can tell):</p>

<pre><code>    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
</code></pre>

python algorithm