MCMC for twin diagnosis

twinsmodel.py

from pymc import Model, MCMC, deterministic, stochastic, Normal, Bernoulli, Beta, distributions
import numpy as np

def twinModelVars( p_identical = 0.5, p_down = 0.01, measurement = 2.5, measurementPrecision = 2.0 ):
    """ Stylized model for twin Down Syndrome likelihood after a blood sample is taken from mother
        containing a mixture of indicators from each foetus

    :param p_identical: (Prior) probability of identical twins versus fraternal (0.1 for population)
    :param p_down:      (Prior) probability of Down Syndrome
    :param measurement: Measured marker value
    :param measurementPrecision: Precision of marker measurement
    :return: List of variables than can be used in a pymc.Model
    """
    identical     = Bernoulli('identical', p = p_identical )   # Are twins identical or fraternal?
    down1         = Bernoulli('down1', p = p_down )            # Will first have Down Syndrome?
    down2_ind     = Bernoulli('down2_ind', p = p_down )        # Would the second, if independent?
    downMarker1   = Normal('downMarker1', value = 3.0, mu = 3.0, tau = 1.0 )           # Stylized lab results, Down Syndrome hypothesis
    normalMarker1 = Normal('normalMarker1', value = 0.0, mu = 0.0, tau = 1.0)          # Stylized lab results, Null
    downMarker2   = Normal('downMarker2', value = 3.0, mu = 3.0, tau = 1.0 )           # Independent lab results, DS
    normalMarker2 = Normal('normalMarker2', value = 0.0, mu = 0.0, tau = 1.0 )         # Independent lab results, Null
    mixingCoef    = Beta('mixingCoef',value = 0.5, alpha = 50,beta = 100 )

    @deterministic( dtype=bool, plot = False )
    def down2( identical = identical, down1 = down1, down2_ind = down2_ind ):
        # Indicator variable: is the second twin Down Syndrome?
        if identical:
            return down1
        else:
            return down2_ind

    @deterministic( dtype=float, plot = False )
    def marker1( down1 = down1, downMarker1 = downMarker1, normalMarker1 = normalMarker1 ):
        # Lab result for first (which we can't know, as results are commingled with that of its sibling)
        if down1:
            return downMarker1
        else:
            return normalMarker1

    @deterministic( dtype=float, plot = False )
    def marker2( down2 = down2, downMarker2 = downMarker2, normalMarker2 = normalMarker2 ):
        # Lab result for second (which we can't observe either)
        if down2:
            return downMarker2
        else:
            return normalMarker2

    @stochastic( dtype = float, observed = True )                    # Noisy lab result we do observe (avg. of twin's readings)
    def measurement( value = 1.0, marker1 = marker1, marker2 = marker2, mixingCoef = mixingCoef ):
        mixedMarker =  mixingCoef*marker1 + (1-mixingCoef)*marker2
        return distributions.normal_like( x = value, mu = mixedMarker, tau = measurementPrecision )

    @deterministic( dtype = bool)
    def bothDown( down1 = down1, down2=down2 ):
        return down1 and down2

    return locals()

twinsModel      = Model( twinModelVars() )
mc              = MCMC( twinsModel )
mc.sample( iter = 500 )

def main():
    # Show booleans in basis points
    for vn in [ "down1","down2","bothDown"]:
        samples = mc.trace(vn).gettrace()
        print "Mean for " + vn + " is " + str( 10000.*np.mean( samples ) ) + " bp "

    # Show markers
    for vn in ["downMarker1","downMarker1","marker1", "marker2","downMarker1","downMarker2","mixingCoef"]:
        markers = mc.trace(vn).gettrace()
        print "Mean for " + vn + " is " + str( np.mean( markers ) )

    betas = mc.trace("mixingCoef").gettrace()
    print "Here are the betas..."
    print betas
    print "WTF?"