# EXAMPLE 11: Introduction to Statistics for Astrophysicists
# SS 2012
# JOCHEN WELLER

x0=1
#fwhm = 5 = 2.35*sigma
w=2.13
n0=33.3333
A=1
B=2

#First generate random data 
#Datum is Gaussian, and measurement Poisson

xmin=-7
xmax=7
x=seq(xmin,xmax,by=1.0)

d=n0*(A*exp(-0.5*(x-x0)^2/w^2)+B)
N=d
imax=length(d)
i=1
while (i <= imax) {N[i]=rpois(1,d[i])
	i=i+1}
	
# The minimum chi2 for later	
logepmin=sum(N*log(d)-d)	
	

# posterior function for arbitrary A and B
posterior=function (A,B) 
{ summe=0.0*A
	for (i in 1:length(x)) {
	di=n0*(A*exp(-0.5*(x[i]-x0)^2/w^2)+B)
	summe=summe+N[i]*log(di)-di
	}
	return(exp(summe-logepmin))
}
pAmarg=function(Ai) {I=integrate(posterior,0.0,Inf,A=Ai)
	return(I$value)
	}

pAplot=function(Ai) {dumplot=Ai
	for (i in (1:length(Ai))) {
	     Aarg=Ai[i]*1.0
	     dumplot[i]=pAmarg(Aarg)
	            }
	return(dumplot)}		
	
pBmarg=function(Bi) {I=integrate(posterior,0.0,Inf,B=Bi)
	return(I$value)
	}

pBplot=function(Bi) {dumplot=Bi
	for (i in (1:length(Bi))) {
	     Barg=Bi[i]*1.0
	     dumplot[i]=pBmarg(Barg)
	            }
	return(dumplot)}	
	
Aplot=seq(0.1,3.0,by = 0.02)
Bplot=seq(0.1,3.0,by = 0.02)

normA=1.0/max(pAplot(Aplot))
plot(Aplot,normA*pAplot(Aplot),type='l',lwd=2,xlab='A',ylab='prob(A|{N_k},I)')
normA2=1.0/max(posterior(Aplot,2.0))
lines(Aplot,normA2*posterior(Aplot,2.0),lty=3,lwd=2)	

quartz()
normB=1.0/max(pBplot(Bplot))
plot(Bplot,normB*pBplot(Bplot),type='l',lwd=2,xlab='B',ylab='prob(B|{N_k},I)')