Wednesday, July 31, 2013

Computational Molecular Evolution

OK, I trying to do the quiz for Bayesian Analysis.  I'm stuck on the first problem.  I'm doing the calculations as I understand them but instead of getting fractions for likelihood I'm getting whole numbers.  That can't be. I coded:
as:

def fac(n):
    if n == 1: return 1
    else: return n * fac(n-1)

def prob(h,n,p):
    return float((fac(n)*p*(1-p)))/float(fac(h)*fac(n-h))

I'll figure it out. I thought I did it correctly. Clearly I did not.  There is a math function factorial.  I was too lazy to import the math library.  It seems to me that the number of ways one can get heads from N flips doesn't really express the probability of it. 

Ah, I see the problem.  That superscript really does mean "Raised to the power".  Now I can move on to more misunderstandings on my part.

Interestingly, Wikipedia gives a different but similar function where r is the uniform prior distribution, somewhat analogous to p in the prior formula.  The Wikipedia example following this helped me identify my error.
While I know some about probabilities I hadn't dealt with conditional probabilities.  It looks quite powerful and useful.

On another topic.  When using Python one needs to know math.log(x) gives the natural log of x rather than the base 10 log of x.  The base 10 log of x is math.log10(x). 

It gets stranger.  math.log1p(x) is the natural log of 1+x.

No comments:

Post a Comment