Embedded Files
strangerland®gmail.com
<iframe src="https://docs.google.com/file/d/0B2lzqB92j17jbDVfdHpRVnNzSEE/preview" style="width:1400px; height:500px;" frameborder="0"></iframe>
Imagine we had a set of N distinguishable coins any of which may be fair or unfair, thus each has a two state probability schema. We identify heads with the boolean value True and tails with the boolean value False. Thus we have N independent boolean random variables. Assume we have an assistant toss the coins, and compute a set of Boolean polynomial functions of the outcomes of the toss. We are allowed to see the values of the Boolean polynomials and are given the form of each polynomial, but we are not allowed to see the outcomes of the individual coins (unless of course any of the functions are the identity function for a given coin). If we had to make an engineering estimate of the probability schemas for each coin we might collect statistics on the Boolean functions for a large number of coin tosses and estimate each functions two state probability schema once the probability polynomials for each boolean function is computed. These can be computed recursively according to:
P(A ∧ B) = P(A|B)P(B) = P(B|A)P(A)
P(A ∨ B) = P(A) + P(B) - P(A ∧ B)
P(¬A) = P(1 ⊻ A) = 1 - P(A)
P(A → B) = 1 - P(A) + P(B) - P(¬A ∧ B)
P(A ↔ B) = P(A ∧ B) + P(¬A ∧¬ B)
P(A ⊻ B) = 1 - P(A ∧ B) - P(¬A ∧¬ B) = P(A) + P(B) - 2P(A ∧ B)
if A and B are independent then P(A ∧ B) = P(A)P(B) which we might suppose for generators A, B
if A and B are mutually exclusive then P(A ∨ B) = P(A) + P(B) which we suppose for atoms A, B
P(A) = P(A|B)P(B) + P(A|¬B)P(¬B)
Bayes Theorem:
P(B|A) = P(A|B)P(B)/P(A)
If statements X, Y, Z... are independent use x = P(X) etc
P(X ∧ Y) = xy
P(X ∨ Y) = x + y - xy
P(X → Y) = 1 - x + xy
P(X ⊻ Y) = -2xy + x + y
P(X ↔ Y) = 2xy - x - y + 1
P(¬X) = 1 - x
P(X ∧ ¬X) = 0
P(X ∨ ¬X) = 1
P(X ⊼ Y) = -xy + 1
P(X ⊽ Y) = xy - x - y + 1
We could then make a maximum entropy "estimate" using Lagrange multipliers. Or we could use a least squares "estimate" possibly. Neither estimate need be close however to the true probabilities for each coin. The least squares method isn't subject to inconsistent or zero dimensional systems of constraint polynomials (I don't think it is unlike with max entropy)
new discovery - I can insert equation images into google sites
nifty
Anyway the problem I have to think about here is that if certain boolean polynomials are assumed true or if the corresponding probability polynomials are set to specified probabilities as constraints, then doesn't that violate the assumption that the generators were independent random variables which was the assumption used to compute the form of said polynomials? Well we assume thereexists certain two stae prob schemes for each independent generator (ind. random variable) then any given prob poly will evaluate to some probability. So if we have a consistent set of prob polynomials, there should exist solns for the generator probs (which need to be on [0,1] to be acceptable) the generators can still be assumed independent with these 2-state prob schemas. Then which solns maximize the entropy. If over specified system of polynomials, use least squares est of generator probs, still can assume independence. I think we're ok here.
SAGE and python code.
#******************************************************************************************
def ProbM(myBP): # myBP is a boolean polynomial being converted to probabilty expression Method 2
if type(myBP) == sage.rings.integer.Integer:
if myBP%2 == 0:
return QQ(0)
else:
return QQ(1)
mystr = str(myBP.variables()) # converting all uppers to lower and lowers to uppercase
newstr = ""
for l in mystr:
if l.isupper():
newstr = newstr + "1" # gotta mark them as upper or lower
else:
newstr = newstr + "0"
NewGens = ""
for ii in range(len(mystr)):
if newstr[ii] is "1":
NewGens = NewGens + mystr[ii].lower()
else:
NewGens = NewGens + mystr[ii].upper() # if boolean variable XcGGt1u then, prob var is xCggT1U etc.
NewGens = NewGens[1:len(NewGens) - 1] # generators for probabilty polynomial assuming boolean generators are independent
if len(myBP.variables()) == 1:
NewGens = NewGens.replace(',','')
QPR = PolynomialRing(QQ,len(myBP.variables()),NewGens.replace(', ',','),order='lex') # needs order to be lex to work
QPR.inject_variables()
nn = myBP.nvariables() # number of generators in boolean polynomial
BPvars = myBP.variables()
BPmons = myBP.monomials()
atoms = 2**nn # number of atoms in corresponding boolean algebra
# should abort this with error, if too many variables in polynomial
Mvect = [0]*atoms
for jj in range(len(BPmons)):
mm = 0
for kk in range(len(BPmons[jj].variables())):
mm = mm + 2**BPvars.index(BPmons[jj].variables()[kk])
Mvect[mm] = 1
Mvect.reverse()
# Mvect order for instance is (X111,X110,X101,X100,X011,X010,X001,X000) has two nonequiv intepretations
# as monomials or atoms example was for three variables
Mvect = vector(Mvect)
MGinv = matrix(QQ,[1])
for ii in range(1,nn+1):
MGinv = block_matrix(QQ,2,2,[MGinv,-MGinv,0,MGinv],subdivide=False) # creating inverse of multigrade operator 'AND' matrix MG
# modulo 2 it's an involution however...
# MG and MGinv are upper triangular, MG is like Sierpinski triangle see
# https://commons.wikimedia.org/wiki/File:Multigrade_operator_AND.svg
Avect = (MGinv*Mvect)%2
Pvect = MGinv*Avect # convert from boolean ring (mons) to boolean algebra (atoms) then to probability polynomial
# ((MGinv*Mvect)%2) is the atoms vector equivalent to the monnomials vector Mvect
# Pvect are the coefficients of prob poly in CPR
newpoly = QPR(1)
for ii in range(0,len(QPR.gens())):
newpoly = newpoly*(QPR(1) + QPR.gen(ii)) # to create all possible monomials in CPR
return Pvect*vector(newpoly.monomials()) # inner product of the two vectors to create prob poly
#******************************************************************************************
#Prob.py caculate the probabilty polynomials Method 1
def Prob(myBP):
if type(myBP) == sage.rings.integer.Integer:
if myBP%2 == 0:
return QQ(0)
else:
return QQ(1)
BPR = myBP.parent()
mystr=str(myBP.variables()) #converting all uppers to lower and lowers to uppercase
newstr=""
for l in mystr:
if l.isupper():
newstr=newstr+"1" #gotta mark them as upper or lower
else:
newstr=newstr+"0"
NewGens=""
for ii in range(len(mystr)):
if newstr[ii] is "1":
NewGens=NewGens+mystr[ii].lower()
else:
NewGens=NewGens+mystr[ii].upper()
NewGens = NewGens[1:len(NewGens) - 1] #generators for probabilty polynomial assuming boolean generators are independent
# should abort this with error, if too many variables in polynomial
if len(myBP.variables()) == 1:
NewGens = NewGens.replace(',','')
NewGens.replace(', ',',')
ProbRing = PolynomialRing(QQ,len(myBP.variables()),NewGens,order='lex')#order=BPR.term_order())
ProbRing.inject_variables()
return ProbRecurse(myBP,BPR,ProbRing)
#*******************************************************************************************
#ProbRecurse.py this is the recursive part
def ProbRecurse(myBP,BPR,ProbRing):
MyString = str(myBP)
if MyString == '1':
return 1
elif MyString == '0':
return 0
elif MyString.find(' 1') != -1:
return 1 - ProbRecurse(myBP + 1,BPR,ProbRing)
elif MyString.find('+') == -1:
return ProbRing(MyString.lower())
else:
kk = MyString.find(' + ')
myPoly1 = BPR(MyString[:kk])
myPoly2 = BPR(MyString[kk+3:])
return ProbRecurse(myPoly1,BPR,ProbRing) + ProbRecurse(myPoly2,BPR,ProbRing) - 2*ProbRecurse(myPoly1*myPoly2,BPR,ProbRing)
load('ProbRecurse.py')
load('Prob.py')
import numpy
import scipy.optimize
KDELTA = lambda A, B: A.parent(A == B)
NOT = lambda A: A.parent(1) + A
XOR = lambda A, B: A + B
OR = lambda A, B: A + B + A*B
IF = lambda A, B: A.parent(1) + A + A*B
IFF = lambda A, B: A.parent(1) + A + B
AND = lambda A, B: A*B
NAND = lambda A, B: A.parent(1) + A*B
NOR = lambda A, B: A.parent(1) + A + B + A*B
FORALL = lambda mylist: reduce(AND, mylist)
EXISTS = lambda mylist: reduce(OR, mylist)
UNIQUE = lambda mylist: reduce(XOR,[FORALL([A + NOT(KDELTA(A,B)) for A in mylist]) for B in mylist])
GIVEN = lambda A, B: Prob(AND(A,B))/Prob(B)
CONVERT = lambda f: reduce(lambda a,b:a+b,[reduce(lambda a,b:a*b,[(m.degree(g)+1+g) for g in BPR.gens()]) for m in f.monomials()])
BPR = BooleanPolynomialRing(6,'X', order='degneglex')
X = BPR.gens() #generators X[n] are independent random variables that can take on values 0 or 1
NewGens = str(X).lower()
NewGens = NewGens[1:len(NewGens) - 1]
x = var(NewGens) #x = Pr(X)
C1=AND(OR(X[0],X[1]),OR(X[0],X[2]))
C2=IF(AND(X[1],X[2]),IF(X[4],X[5]))
C3=NOT(IF(IF(X[0],X[3]),NOT(IF(X[3],X[5]))))
C4=NOT(X[5])
prems=[C1,C2,C3,C4]
if BPR(0) in prems:
print "explicit contradiction in premises"
prems = list(set(prems)) # remove dup
if BPR(1) in prems: prems.remove(BPR(1)) # remove tautologies
prb = 1.0
constraint=[Prob(FORALL(prems)) - prb]
HH = lambda s: sum([s[ii]*ln(s[ii])+(1-s[ii])*ln(1-s[ii]) for ii in range(len(s))])
dHH = lambda s: [ln(s[ii])-ln(1-s[ii]) for ii in range(len(s))]
str1='{\'type\': \'eq\', \'fun\': lambda s: ((constraint['
str2='])(*s))},'
mystr='['
for ii in range(len(constraint)):
mystr=mystr+str1+str(ii)+str2
mystr=mystr+']'
cons = eval(mystr)
numpy.set_printoptions(precision=15)
y0 = [.5 for xx in x]
bnds = tuple([(0.0, 1.0) for xx in x])
ssoln = scipy.optimize.minimize(HH, y0, jac=dHH, method='SLSQP', bounds=bnds, constraints=cons)
print('ssoln.x');print(ssoln.x);print(' ')
print('entropy ssoln.x');print(-HH(ssoln.x)/RR(ln(2)));print(' ')
pvar = lambda s: sum([(2*s[i]-1)^2 for i in range(len(s))])
dpvar = lambda s: [4*(2*s[i]-1) for i in range(len(s))]
pvarsoln = scipy.optimize.minimize(pvar, y0, jac=dpvar, method='SLSQP', bounds=bnds, constraints=cons)
print('pvarsoln.x');print(pvarsoln.x);print(' ')
print('entropy pvarsoln.x');print( -HH(pvarsoln.x)/RR(ln(2)));print(' ')
#output
ssoln.x [ 5.693223670277803e-13 9.999999999996521e-01 9.999999999993416e-01 7.759903830617532e-13 1.550981565401344e-13 1.461053500406706e-13] entropy ssoln.x 1.12073139297e-10 pvarsoln.x [ 5.693223670277803e-13 9.999999999996521e-01 9.999999999993416e-01 7.759903830617532e-13 1.550981565401344e-13 1.461053500406706e-13] entropy pvarsoln.x 1.12073139297e-10
#************************************************************************************************************
load('ProbRecurse.py')
load('Prob.py')
import numpy
import scipy.optimize
from operator import add
from operator import mul
KDELTA = lambda A, B: A.parent(A == B)
NOT = lambda A: A.parent(1) + A
XOR = lambda A, B: A + B
OR = lambda A, B: A + B + A*B
IF = lambda A, B: A.parent(1) + A + A*B
IFF = lambda A, B: A.parent(1) + A + B
AND = lambda A, B: A*B
NAND = lambda A, B: A.parent(1) + A*B
NOR = lambda A, B: A.parent(1) + A + B + A*B
FORALL = lambda mylist: reduce(AND, mylist)
EXISTS = lambda mylist: reduce(OR, mylist)
UNIQUE = lambda mylist: reduce(XOR,[FORALL([A + NOT(KDELTA(A,B)) for A in mylist]) for B in mylist])
GIVEN = lambda A, B: Prob(AND(A,B))/Prob(B)
convert = lambda f: reduce(lambda a,b:a+b,[reduce(lambda a,b:a*b,[(m.degree(g)+1+g) for g in BPR.gens()]) for m in f.monomials()]) #self inverse convert poly in gens to poly in atoms (missing var x implies, interpret as ~x)
# convert monomial terms A^i*B^j*C^k --> (A+1+i)*(B+1+j)*(C+1+k) in Boolean Ring i,j,k=0,1
objects = 6
BPR = BooleanPolynomialRing(objects^2,'X',order='degneglex')
BPR.inject_variables()
print(str(objects^2) + ' Free Boolean Algebra generators: X0 = R[0,0], X1 = R[0,1], X2 = R[0,2],..., X' + str(objects^2-1) + ' = R[' + str(objects-1) + ',' + str(objects-1) + ']')
X = BPR.gens()
NewGens = str(X).lower()
NewGens = NewGens[1:len(NewGens) - 1]
x = var(NewGens)
R = matrix(BPR, objects, objects, BPR.gens())
prems = [];reflex=[];symmtrc=[];trans=[];direct=[]
for i in range(0, objects): #reflexive
reflex.append(R[i,i])
for i in range(0,objects): #symmetric
for j in range(0, objects):
symmtrc.append(IF(R[i,j],R[j,i]))
direct.append(IF(R[i,j],NOT(R[j,i])))
symmtrc=list(set(symmtrc))
if BPR(1) in symmtrc: symmtrc.remove(BPR(1)) # get rid on tautologies
for i in range(0, objects): #transitive
for j in range(0, objects):
for k in range(0, objects):
trans.append(IF(AND(R[i,j],R[j,k]),R[i,k]))
#print(len(trans))
trans = list(set(trans))
M=matrix([[125,25,5,1],[64,16,4,1],[27,9,3,1],[8,4,2,1]]);M
A=vector([81,37,13,3]);A
M.inverse()
B=M.inverse()*A # B=(1, -2, 1, 1)
var('t')
mypoly=B*vector([t^3,t^2,t,1]);mypoly # gives number of unique transitive polynomials
prems=reflex+symmtrc+trans
print(len(reflex));print(len(symmtrc));print(len(trans));print(len(prems))
if BPR(0) in prems:
print "explicit contradiction in premises"
prems = list(set(prems)) # remove dup
prems = [p+1 for p in prems] #make true
#if BPR(1) in prems: prems.remove(BPR(1))
REI = ideal(prems); print('REI.gens() = ' + str(len(REI.gens())))
#GB = REI.groebner_basis()
QR = BPR.quotient_ring(REI)
print(prems);print(len(prems))
matrix(objects,objects,QR.gens())
#output
Defining X0, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14, X15, X16, X17, X18, X19, X20, X21, X22, X23, X24, X25, X26, X27, X28, X29, X30, X31, X32, X33, X34, X35 36 Free Boolean Algebra generators: X0 = R[0,0], X1 = R[0,1], X2 = R[0,2],..., X35 = R[5,5] [125 25 5 1] [ 64 16 4 1] [ 27 9 3 1] [ 8 4 2 1] (81, 37, 13, 3) [ 1/6 -1/2 1/2 -1/6] [ -3/2 5 -11/2 2] [ 13/3 -31/2 19 -47/6] [ -4 15 -20 10] t t^3 - 2*t^2 + t + 16 30 151 187 REI.gens() = 187 [X10*X13*X16 + X10*X13, X15*X26*X27 + X15*X26, X0*X2*X12 + X2*X12, X19*X25*X27 + X19*X27, X23*X33 + X33, X10*X11*X34 + X11*X34, X13*X25*X26 + X13*X26, X10*X25 + X10, X5*X12*X17 + X5*X12, X10*X11*X29 + X10*X29, X3*X4*X22 + X3*X22, X16*X26 + X26, X4*X12*X16 + X4*X12, X0*X4*X24 + X4*X24, X14 + 1, X1*X18*X19 + X1*X18, X1*X3*X19 + X3*X19, X13*X19*X20 + X13*X20, X9*X19 + X9, X11*X31*X35 + X11*X31, X5*X30 + X30, X6*X11*X30 + X11*X30, X13*X16*X25 + X16*X25, X13*X17*X31 + X17*X31, X3*X12*X15 + X3*X12, X16*X26*X28 + X16*X26, X24*X29*X30 + X29*X30, X12*X24*X26 + X12*X26, X18*X22*X24 + X22*X24, X8*X9*X15 + X8*X15, X5*X30 + X5, X20*X22*X26 + X22*X26, X22*X23*X34 + X23*X34, X1*X5*X31 + X5*X31, X26*X29*X32 + X29*X32, X26*X32*X34 + X26*X34, X21 + 1, X9*X19*X21 + X9*X19, X12*X15*X18 + X15*X18, X8*X10*X16 + X8*X16, X27*X33*X34 + X27*X34, X29*X34*X35 + X29*X34, X7*X9*X19 + X9*X19, X17*X20*X23 + X17*X20, X14*X16*X26 + X16*X26, X11*X13*X17 + X11*X13, X5*X18*X23 + X5*X18, X6*X9*X18 + X9*X18, X2*X6*X8 + X2*X6, X10*X25 + X25, X9*X11*X33 + X11*X33, X4*X24 + X24, X2*X24*X26 + X2*X24, X11*X25*X29 + X11*X25, X17*X26*X29 + X17*X26, X1*X12*X13 + X1*X12, X15*X20*X21 + X15*X20, X6*X24*X25 + X6*X25, X8*X10*X26 + X10*X26, X9*X19 + X19, X1*X30*X31 + X1*X30, X9*X11*X23 + X9*X23, X2*X18*X20 + X2*X18, X28 + 1, X3*X24*X27 + X3*X24, X2*X5*X32 + X5*X32, X3*X18 + X3, X16*X17*X34 + X17*X34, X2*X12 + X12, X28*X29*X34 + X29*X34, X3*X18 + X18, X8*X13 + X8, X17*X32*X35 + X17*X32, X0*X1*X6 + X1*X6, X9*X13*X15 + X9*X13, X5*X6*X11 + X5*X6, X10*X31*X34 + X10*X31, X3*X5*X23 + X3*X23, X16*X32*X34 + X16*X32, X8*X25*X26 + X8*X25, X5*X30*X35 + X5*X30, X21*X22*X27 + X22*X27, X6*X12*X13 + X6*X13, X15*X32*X33 + X15*X32, X7*X8*X13 + X8*X13, X13*X15*X19 + X15*X19, X8*X11*X32 + X11*X32, X16*X17*X29 + X16*X29, X35 + 1, X20*X32*X33 + X20*X33, X8*X19*X20 + X8*X19, X1*X2*X13 + X2*X13, X2*X4*X26 + X4*X26, X18*X30*X33 + X18*X33, X12*X30*X32 + X12*X32, X13*X31*X32 + X13*X32, X14*X15*X20 + X15*X20, X22*X27 + X22, X21*X23*X33 + X23*X33, X20*X23*X32 + X23*X32, X16*X20*X22 + X16*X20, X4*X5*X34 + X5*X34, X8*X13 + X13, X19*X31*X33 + X19*X33, X25*X31*X34 + X25*X34, X17*X32 + X32, X27*X29*X33 + X29*X33, X15*X20 + X20, X1*X5*X11 + X1*X11, X12*X18*X20 + X12*X20, X2*X12*X14 + X2*X12, X15*X17*X23 + X15*X23, X2*X5*X17 + X2*X17, X24*X30*X34 + X24*X34, X10*X19*X22 + X10*X19, X3*X6*X9 + X3*X6, X3*X5*X33 + X5*X33, X12*X17*X30 + X17*X30, X1*X6 + X6, X1*X3*X9 + X1*X9, X6*X18*X19 + X6*X19, X3*X18*X21 + X3*X18, X2*X4*X16 + X2*X16, X9*X10*X27 + X10*X27, X2*X3*X15 + X2*X15, X8*X31*X32 + X8*X31, X15*X16*X27 + X16*X27, X23*X33 + X23, X10*X25*X28 + X10*X25, X8*X11*X17 + X8*X17, X9*X10*X22 + X9*X22, X11*X19*X23 + X11*X19, X4*X5*X29 + X4*X29, X17*X32 + X17, X11*X31 + X11, X4*X30*X34 + X4*X30, X15*X17*X33 + X17*X33, X3*X30*X33 + X3*X30, X2*X3*X20 + X3*X20, X15*X20 + X15, X22*X33*X34 + X22*X33, X1*X6*X7 + X1*X6, X2*X12 + X2, X7*X11*X31 + X11*X31, X1*X2*X8 + X1*X8, X2*X30*X32 + X2*X30, X1*X4*X10 + X1*X10, X6*X30*X31 + X6*X31, X1*X4*X25 + X4*X25, X19*X22*X25 + X22*X25, X0*X5*X30 + X5*X30, X6*X10*X24 + X10*X24, X9*X31*X33 + X9*X31, X23*X27*X29 + X23*X27, X16*X26 + X16, X12*X16*X24 + X16*X24, X3*X4*X27 + X4*X27, X6*X8*X12 + X8*X12, 0, X1*X6 + X1, X29*X34 + X29, X18*X23*X30 + X23*X30, X22*X27*X28 + X22*X27, X7*X10*X25 + X10*X25, X22*X27 + X27, X11*X31 + X31, X18*X24*X27 + X18*X27, X1*X24*X25 + X1*X24, X29*X34 + X34, X4*X24 + X4, X22*X23*X29 + X22*X29, X8*X9*X20 + X9*X20, X5*X24*X29 + X5*X24, X8*X13*X14 + X8*X13, X9*X25*X27 + X9*X25, X19*X23*X31 + X23*X31, X0 + 1, X0*X3*X18 + X3*X18, X23*X33*X35 + X23*X33, X7 + 1, X25*X29*X31 + X29*X31, X15*X16*X22 + X15*X22, X14*X17*X32 + X17*X32, X4*X24*X28 + X4*X24, X4*X6*X10 + X4*X6, X20*X26*X27 + X20*X27, X4*X18*X22 + X4*X18] 187 [ 1 X1 X2 X3 X4 X5] [ X1 1 X8 X9 X10 X11] [ X2 X8 1 X15 X16 X17] [ X3 X9 X15 1 X22 X23] [ X4 X10 X16 X22 1 X29] [ X5 X11 X17 X23 X29 1]
Page updated
Google Sites
Report abuse