Skip to content

Instantly share code, notes, and snippets.

@notbanker
Created August 29, 2014 20:39
Show Gist options
  • Save notbanker/1375bca2e5283c5a4d1f to your computer and use it in GitHub Desktop.
Save notbanker/1375bca2e5283c5a4d1f to your computer and use it in GitHub Desktop.
MCMC for twin diagnosis
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?"
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment