Help:=proc(): print(`Beta1P(DA,n), Beta2P(DA,n), PTn(DA,p1,p2) , GibbsS(DA,M1,M2), AvBay(DA) , AvGibbs(DA,M1,M2) `): end:


#From C11.txt; Cross-over
Help11:=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:


####New stuff


#Beta1P(DA,n): parameters for the Beta distribution for p1 with data set DA and cross-over n
Beta1P:=proc(DA,n) local i,j,a:
a:=add(add(DA[j][i],i=1..n),j=1..nops(DA)):
a+1,n*nops(DA)-a+1:
end:


#Beta2P(DA,n): parameters for the Beta distribution for p2 with data set DA and cross-over n
Beta2P:=proc(DA,n) local i,j,a,N:
N:=nops(DA[1]):

a:=add(add(DA[j][i],i=n+1..N),j=1..nops(DA)):

a+1,(N-n)*nops(DA)-a+1:
end:


#PTn(DA,p1,p2): The normalized proba. table for n with data set DA and prob. p1, p2
PTn:=proc(DA,p1,p2) local N,n,L,i,j:
N:=nops(DA[1]):
L:=[seq(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)),
n=1..N)]:
L/add(L):
end:


#Gibbs(DA,M1,M2): inputs a data set DA, and pos. integers M1 and M2 outputs M2 samples of (p1,p2,n) in
#the cross-over problem after a burn-out period of M1
GibbsS:=proc(DA,M1,M2) local N,i,p1,p2,n,X1,X2,S:
N:=nops(DA[1]):

n:=rand(1..N)():
for i from 1 to M1 do

p1:=Sample(RandomVariable(BetaDistribution(Beta1P(DA,n))),1)[1]:
p2:=Sample(RandomVariable(BetaDistribution(Beta2P(DA,n))),1)[1]:
n:=trunc(Sample(RandomVariable(ProbabilityTable(PTn(DA,p1,p2))),1)[1])	:

od:

S:=[];

for i from 1 to M2 do

p1:=Sample(RandomVariable(BetaDistribution(Beta1P(DA,n))),1)[1]:
p2:=Sample(RandomVariable(BetaDistribution(Beta2P(DA,n))),1)[1]:
n:=trunc(Sample(RandomVariable(ProbabilityTable(PTn(DA,p1,p2))),1)[1])	:

S:=[op(S),[n,p1,p2]]:
od:
S:

end:


#AvBay(DA): Given a data set DA outputs the expected n, p1, and p2 exactly
AvBay:=proc(DA) local p1,p2:
[MARGn(DA)[2], MARGp1(DA,p1)[2],MARGp2(DA,p2)[2]]:
end:

#AvGibbs(DA,M1,M2): Given a data set DA outputs the expected n, p1, and p2 according to Gibbs Sample with (M1,M2)
AvGibbs:=proc(DA,M1,M2) local S:
S:=GibbsS(DA,M1,M2):
evalf(add(S)/nops(S)):
end: