help with some correlated binary data

[janus]

Hi All

I am observing a number of binary processes, and I can see that they
are not completely uncorrelated. And the p of the processes is not
identical

So I assume a very simple model. My system can be in a good or bad
state. If the system is in the good state each process has some (high)
probability of success, if system is in the bad state each process has
some (lower) probability of success.

So i have written the following script... that obviously dos not
works. One problem is that way that I currently use 'xGlobal' I can
see that my approach is not gonna work, but I don not know what to do
instead. I would think that 'xGlobal' only assumes one value for each
iteration of the simulation, but that there should be N=28 values for
each iteration...

Any input is much appreciated, also if I can do this in a completely
different way? :)

Thanks in advance
best Janus


import numpy as np
from pymc import *
import pylab as pl

# very correlated data
meas = np.array((
[1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,0,1,1,1,1,0,0,1,1,1,1,1,0],
[1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,1,1,0],
[1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,0,1,1,1,1,0,0,1,1,1,1,1,0],
[1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,0,1,1,1,0,0,1,1,1,1,1,1,1]
))

N = np.size(meas,1)
noErasures = np.prod(meas,0)
onlyErasures = np.sum(meas,0) == 0
print 'no erasures with prob = ' + str( np.mean(noErasures) )
print 'only erasures with prob = ' + str( np.mean(onlyErasures) )

pGlobal = pymc.Uniform('pGlobal', lower = .5, upper = 1)
xGlobal = pymc.Bernoulli('xGlobal', p = pGlobal)

p11 = list() #prob of succes when in the global good state
p01 = list() #prob of succes when in the global bad state
x = list()

for i in range(np.shape(meas)[0]):
p11.append( pymc.Uniform("p11"+str(i), lower=0.5, upper=1) )
p01.append( pymc.Uniform("p01"+str(i), lower=0.01, upper=1 ) )

@deterministic
def probNode(p11=p11[i], p01=p01[i], xGlobal=xGlobal):
if xGlobal:
return p11
else:
return p01

x.append( pymc.Bernoulli("x"+str(i), p = probNode, value =
meas[i], observed=True) )


@deterministic
def pNoErasures(pGlobal=pGlobal, p11=p11, p01=p01):
return pGlobal * np.prod(p11) + (1-pGlobal) * np.prod(p01)

@deterministic
def pOnlyErasures(pGlobal=pGlobal, p11=p11, p01=p01):
return pGlobal * (1-np.prod(p11)) + (1-pGlobal) * (1-np.prod(p01))

xNoErasure = pymc.Bernoulli('xNoErasure', p = pNoErasures, value =
noErasures, observed=True)
xOnlyErasure = pymc.Bernoulli('xOnlyErasure', p = pOnlyErasures, value
= onlyErasures, observed=True)

[janus]

After some tinkering I have the following that runs, but still do not
work. The following lines is supposed to produce a function for each
of my 4 nodes, can someone confirm that this is what happens? As I
only get one plot for probNode, this could be my problem?

"
@deterministic
def probNode(pGlobal=pGlobal, p11=p11[i], p01=p01[i]):
return pGlobal*p11 + (1.-pGlobal)*p01

"


import numpy as np
from pymc import *
import pylab as pl

# very correlated data
meas = np.array((
[1,1,1,1,1,1,1,0,0,0,1,1,1,1,1,0,1,1,1,1,0,0,1,1,1,1,1,0],
[1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,1,1,0],
[1,1,1,1,1,1,1,0,0,0,0,1,1,1,1,0,1,1,1,1,0,0,1,1,1,1,1,0],
[1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,0,1,1,1,0,0,1,1,1,1,1,1,1]
))

N = np.size(meas,1)
noErasures = np.prod(meas,0)

onlyErasures = np.zeros(np.size(noErasures), dtype=int)
onlyErasures[np.sum(meas,0) == 0 ] = 1
print 'no erasures = ' + str( noErasures ) + ', with prob = ' +
str(np.mean(noErasures))
print 'only erasures = ' + str( onlyErasures ) + ', with prob = ' +
str(np.mean(onlyErasures))

aGlobal = pymc.Exponential('aGlobal', beta = 0.01, trace=False)
bGlobal = pymc.Exponential('bGlobal', beta = 0.01, trace=False)
pGlobal = pymc.Beta('pGlobal', alpha = aGlobal, beta = bGlobal, value=.
7)

a11 = list()
b11 = list()

p11 = list() #prob of succes when in the global good state

a01 = list()
b01 = list()

p01 = list() #prob of succes when in the global bad state
x = list()

pNode = list()

for i in range(np.shape(meas)[0]):

a11.append( pymc.Exponential("a11"+str(i), beta = 0.01,
trace=False) )
b11.append( pymc.Exponential("b11"+str(i), beta = 0.01,
trace=False) )
p11.append( pymc.Beta("p11"+str(i), alpha = a11[i], beta = b11[i],
value=.7) )

a01.append( pymc.Exponential("a01"+str(i), beta = 0.01,
trace=False) )
b01.append( pymc.Exponential("b01"+str(i), beta = 0.01,
trace=False) )
p01.append( pymc.Beta("p01"+str(i), alpha = a01[i], beta = b01[i],
value=.5) )

@deterministic
def probNode(pGlobal=pGlobal, p11=p11[i], p01=p01[i]):
return pGlobal*p11 + (1.-pGlobal)*p01

x.append( pymc.Bernoulli("x"+str(i), p = probNode, value =
meas[i], observed=True) )


@deterministic
def pNoErasures(pGlobal=pGlobal, p11=p11, p01=p01):

return pGlobal * np.prod(p11) + (1.-pGlobal) * np.prod(p01)

@deterministic
def pOnlyErasures(pGlobal=pGlobal, p11=p11, p01=p01):

return pGlobal * (1.-np.prod(p11)) + (1.-pGlobal) * (1.-

np.prod(p01))

xNoErasure = pymc.Bernoulli('xNoErasure', p = pNoErasures, value =
noErasures, observed=True)
xOnlyErasure = pymc.Bernoulli('xOnlyErasure', p = pOnlyErasures, value
= onlyErasures, observed=True)

[janus]

for any future interest here follows one solution
/ Janus

import numpy as np
from pymc import *
import pylab as pl

# very correlated data
meas = np.array((

[0,1,1,1,1,1,1,1,0,0,1,1,1,1,1,0,1,1,1,1,1,0,1,1,1,1,1,1,1,1,1],
[1,0,1,1,1,1,1,0,1,1,0,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1,1,0,1,1,1],
[1,1,0,1,1,1,1,0,1,1,0,1,1,1,1,0,1,1,1,1,1,0,1,1,1,1,1,0,0,1,1],
[1,1,1,0,1,1,1,1,0,0,1,1,1,1,1,0,1,1,1,1,0,1,1,1,1,1,1,0,0,0,1]
))

noErasures = np.prod(meas,0)
onlyErasures = np.zeros(np.size(noErasures), dtype=int)
onlyErasures[np.sum(meas,0) == 0 ] = 1
print 'no erasures = ' + str( noErasures ) + ', with prob = ' +
str(np.mean(noErasures))
print 'only erasures = ' + str( onlyErasures ) + ', with prob = ' +
str(np.mean(onlyErasures))

pGlobal = pymc.Uniform('pGlobal', lower = .5, upper = 1)

p11 = list() #prob of succes when in the global good state
p01 = list() #prob of succes when in the global bad state
x = list()

pNode = list()

for i in range(np.shape(meas)[0]):

p11.append( pymc.Uniform("p11"+str(i), lower=0.5, upper=1) )
p01.append( pymc.Uniform("p01"+str(i), lower=0.01, upper=1 ) )

@deterministic(name = 'probNode'+str(i))

def probNode(pGlobal=pGlobal, p11=p11[i], p01=p01[i]):
return pGlobal*p11 + (1.-pGlobal)*p01

pNode.append(probNode)
x.append( pymc.Bernoulli("x"+str(i), p = pNode[i], value =

meas[i], observed=True) )


@deterministic
def pNoErasures(pGlobal=pGlobal, p11=p11, p01=p01):
return pGlobal * np.prod(p11) + (1.-pGlobal) * np.prod(p01)

@deterministic
def pOnlyErasures(pGlobal=pGlobal, p11=p11, p01=p01):

return pGlobal * np.prod(1.-np.array(p11)) + (1.-pGlobal) *
np.prod(1.-np.array(p01))

xNoErasure = pymc.Bernoulli('xNoErasure', p = pNoErasures, value =
noErasures, observed=True)
xOnlyErasure = pymc.Bernoulli('xOnlyErasure', p = pOnlyErasures, value
= onlyErasures, observed=True)