t.nex

Begin data;
	Dimensions NTAX=4 NCHAR=1;
	Format MISSING=? GAP=- DATATYPE=DNA;
	Matrix
S1 G
S2 C
S3 A
S4 T
	;
END;

begin trees;
    [lnL=-10.264] tree 1  = [&U] (S1:0.05,S2:0.05,(S3:0.05,S4:0.05):0.05);
end;

begin paup;
    set warnreset=no warntree=no;
    set criterion=likelihood;
    [JC] lset nst=1
            basefreq=equal
            rates=equal
            pinvar=0;
    lset userbrlens=yes;
    lset logiter=yes;
    lscores 1;
end;

t.py

#! /usr/bin/env python
# -*- coding: utf-8 -*-

import math

def jc(start, end, tlen):
    if start != end:
        return (0.25 - (0.25 * math.exp((-4/3) * tlen)))
    else:
        return (0.25 + (0.75 * math.exp((-4/3) * tlen)))

def jc_check(start, end, tlen):
    from dendropy.model import discrete
    hky = discrete.Hky85(kappa=1)
    return hky.pij(start, end, tlen=tlen)

# S1 = G, S2 = C, S3 = A, S4 = T
# [lnL=-10.264] tree 1  = [&U] (S1:0.05,S2:0.05,(S3:0.05,S4:0.05):0.05);
def calc(
        external1=2,
        external2=1,
        external3=0,
        external4=3,
        nu=0.05,
        ):
    sum_p = 0.0
    for iidx1, internal12 in enumerate((0,1,2,3)):
        for iidx2, internal34 in enumerate((0,1,2,3)):
            p0 = 0.25
            p1 = jc(external1, internal12, nu)
            p2 = jc(external2, internal12, nu)
            p3 = jc(external3, internal34, nu)
            p4 = jc(external4, internal34, nu)
            p5 = jc(internal12, internal34, nu)
            p = p0 * p1 * p2 * p3 * p4 * p5
            sum_p += p
    print(math.log(sum_p))

calc()