SS2012-StatisticalMethods.git

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
# 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)
 
 
הההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההההה
XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX