#C11.txt; Cross-over
Help:=proc():print(`GD1(p1,p2,n,N) , GD(p1,p2,n,N,M) ,  UnNorPost(DA,p1,p2,n), Post(DA,p1,p2) `):
print(` MARGp1(DA,p1),  MARGp2(DA,p2),  MARGn(DA) `):
end:

with(Statistics):

#GD1(p1,p2,n,N): "Synthetic data" (for testing) Bernoulli(p1) for the first N tosses
#followed by Bernoulli(p2) for the last N-n tosses
GD1:=proc(p1,p2,n,N) local X1,X2,i,L1,L2:

X1:=RandomVariable(Bernoulli(p1)):
X2:=RandomVariable(Bernoulli(p2)):
L1:=Sample(X1,n):
L2:=Sample(X2,N-n):
[seq(round(L1[i]),i=1..n),seq(round(L2[i]),i=1..N-n)]:

end:

#GD(p1,p2,n,N,M): generates M samples from GD1
GD:=proc(p1,p2,n,N,M) local i:
[seq(GD1(p1,p2,n,N),i=1..M)]:

end:

#UnNorPost(DA,p1,p2,n): The unormalized posterior deduced from the data DA
#assuming uniform priors for both p1 and p2
UnNorPost:=proc(DA,p1,p2,n) local N,i,j:
N:=nops(DA[1]):
#Prob(n,p1,p2|DA)~ Prob(DA|n,p1,p2)

mul(
mul(p1^DA[j][i]*(1-p1)^(1-DA[j][i]),i=1..n)*
mul(p2^DA[j][i]*(1-p2)^(1-DA[j][i]),i=n+1..N),
j=1..nops(DA)):
end:

#Post(DA,p1,p2): The  posterior distribution from the data DA
#a list of length N, whose n-th item is the pdf for p1,p2 with cross-over n
#assuming uniform priors for both p1 and p2
Post:=proc(DA,p1,p2) local n,i,L,N,tot:
N:=nops(DA[1]):
L:=[seq(UnNorPost(DA,p1,p2,n),n=1..N)]:

tot:=add(int(int(L[i],p1=0..1),p2=0..1),i=1..N):

L/tot:
end:

#MARGp1(DA,p1): Inputs a data set DA, and a symbol p1, outputs the EXACT (symbolic)
#expression in p1 for the MARGINAL posterior of p1 (where p2 is integrated over, and n  is summed over)
#followed by the Expected (predicted) value for p1
MARGp1:=proc(DA,p1) local N,L,n,p2,f:
N:=nops(DA[1]):
L:=Post(DA,p1,p2):

f:=add( int(L[n],p2=0..1)      , n=1..N):

f,evalf(int(p1*f,p1=0..1)):


end:

#MARGp2(DA,p2): Inputs a data set DA, and a symbol p2, outputs the EXACT (symbolic)
#expression in p2 for the MARGINAL posterior of p2 (where p1 is integrated over, and n  is summed over)
#followed by the Expected (predicted) value for p2
MARGp2:=proc(DA,p2) local N,L,n,p1,f:
N:=nops(DA[1]):
L:=Post(DA,p1,p2):

f:=add( int(L[n],p1=0..1)      , n=1..N):

f,evalf(int(p2*f,p2=0..1)):


end:


#MARGn(DA): Inputs a data set DA,, outputs the probability mass function supported between 1 and N ( that nops(DA[1])), of the Marginal n
#(where p1 and p2 are "averaged out", i,e. intergated, followed by the expected (i.e. predicted) value of the cross-over point,
#judging from the data
MARGn:=proc(DA) local N,L,i,p1,p2	:
N:=nops(DA[1]):
L:=Post(DA,p1,p2):

L:=[seq( int(int(L[i],p1=0..1),p2=0..1)  , i=1..N)]:

evalf(L),evalf(add(i*L[i],i=1..N)):

end: