#implementing  "Bayesian Inference: Gibbs Sampling" by Ilker Yildrim

with(Statistics):

ez:=proc(): print(` GD1(n0,N0,L1,L2), pL1L2n(L1,L2,n,N,DA,a,b) `): 
print(`DA:=[seq(GD1(5,10,1,2,a,b),i=1..30)]: pL1L2n(L1,L2,n,10,DA,1,1);`):
print(`DA:=[seq(GD1(5,10,1,2),i=1..30)]: Post(L1,L2,10,DA,1,1);`):
print(`DA:=[seq(GD1(3,10,2,5),i=1..30)]: Post(L1,L2,10,DA,1,1): MARGn(10,DA,1,1);`):
print(`DA:=[seq(GD1(3,10,2,5),i=1..30)]: Post(L1,L2,10,DA,1,1): MARGL1(L1,10,DA,1,1);`):
print(`DA:=[seq(GD1(3,10,2,5),i=1..30)]: Post(L1,L2,10,DA,1,1): MARGL2(L2,10,DA,1,1);`):
end:


#GD1(n0,N0,L1,L2): inputs positive integers n0, N0, 1<=n0<=N0, and positive numbers L1, L2, outputs a list of integers  of length N0
#obtained by random smpling Poission with parameter L1 up to n0, and L2 after. Try:
#GD1(3,10,1,2);
GD1:=proc(n0,N0,L1,L2)  local i:
[seq(Sample(RandomVariable(Poisson(L1)),1)[1],i=1..n0),seq(Sample(RandomVariable(Poisson(L2)),1)[1],i=n0+1..N0)]:

end:


#pL1L2n(L1,L2,n,N,DA,a,b): The un-normalized trivaratie posterior-distribution of the Ilker Yildim process
#given the data, with prior Gamma(L;a,b) Try:
#DA:=[seq(GD1(5,10,1,2,a,b),i=1..30)]: pL1L2n(L1,L2,n,10,DA,1,1);
pL1L2n:=proc(L1,L2,n,N,DA,a,b) local gu,i1,j1:

gu:=mul(mul(exp(-L1)*L1^DA[i1][j1],j1=1..n)*mul(exp(-L2)*L2^DA[i1][j1],j1=n+1..N),i1=1..nops(DA)):
gu:=gu*1/GAMMA(a)*b^a*L1^(a-1)*exp(-b*L1)*1/GAMMA(a)*b^a*L2^(a-1)*exp(-b*L2)*(1/N):

end:

#Post(L1,L2,N,DA,a,b): The list of length N whose  n-th item is the posterior pdf for L1, L2. Try:
#DA:=[seq(GD1(5,10,1,2),i=1..30)]: Post(L1,L2,10,DA,1,1);
Post:=proc(L1,L2,N,DA,a,b) local n,gu,T,i,x,y:
option remember:
gu:=[seq(pL1L2n(L1,L2,n,N,DA,a,b),n=1..N)]:

T:=add(int(int(subs({L1=x,L2=y},gu[i]),x=0..infinity),y=0..infinity),i=1..N):

gu/T:

end:


#MARGn(N,DA,a,b):  The Marginal posterior distrubution for n, . Try:
#DA:=[seq(GD1(3,10,2,5),i=1..30)]: Post(L1,L2,10,DA,1,1): MARGn(10,DA,1,1);
#followed by the expectation
MARGn:=proc(N,DA,a,b) local gu, L1,L2,i,vu:
gu:=Post(L1,L2,N,DA,a,b):
vu:=[seq(int(int(gu[i],L1=0..infinity),L2=0..infinity),i=1..N)]:
vu,add(i*vu[i],i=1..N):
end:



#MARGL1(L1,N,DA,a,b):  The Marginal posterior distrubution for L1, . Try:
#DA:=[seq(GD1(3,10,2,5),i=1..30)]: Post(L1,L2,10,DA,1,1): MARGL1(L1,10,DA,1,1);
#followed by the expectation
MARGL1:=proc(L1,N,DA,a,b) local gu, L2,i,vu:
gu:=Post(L1,L2,N,DA,a,b):
vu:=add(int(gu[i],L2=0..infinity),i=1..N):
vu,int(L1*vu,L1=0..infinity):
end:


#MARGL2(L2,N,DA,a,b):  The Marginal posterior distrubution for L1, . Try:
#DA:=[seq(GD1(3,10,2,5),i=1..30)]: Post(L1,L2,10,DA,1,1): MARGL2(L2,10,DA,1,1);
#followed by the expectation
MARGL2:=proc(L1,N,DA,a,b) local gu, L2,i,vu:
gu:=Post(L1,L2,N,DA,a,b):
vu:=add(int(gu[i],L1=0..infinity),i=1..N):
vu,int(L2*vu,L2=0..infinity):
end:


SamL1:=proc(n,a,b,DA)  local X,i,j:
X:=RandomVariable(GammaDistribution(a+add(add(DA[i][j],j=1..n),i=1..nops(DA)), b+ n*nops(DA)):
Sample(X,1)[1]:

end: