9/solution9.R (0d1a76b) - SS2012-StatisticalMethods.git

solution9.R

#Sheet 9 Michael Klauser
library(gplots)
library(fields)
Filepath="ProjectedDensity.dat"
data=scan(Filepath, what=list(x=0., D=0., sigma=0.))
data=data.frame(x=data$x, D=data$D, sigma=data$sigma)

plotCI(data$x, data$D, uiw=data$sigma, xlab="Angular Position", ylab="Projected Density")

M=10
N=length(data$x)
theta=rep(0,M)
delta.x=(max(data$x)-min(data$x))/M
for (j in 1:M){
	theta[j]=min(data$x)+delta.x*(j-0.5)
}

kron = function(a,b) {
	if (a==b) {return(1.)} else {return(0.)}
}

T=array(0.,c(N,M))
a=array(0.,c(10,11))
F=array(0.,c(11,11,N))
breaks=seq(min(data$x)-1e-9,max(data$x)+1e-9,length.out=11)
bin = findInterval(data$x,breaks)
data = cbind(data,bin)

for (k in 1:N) {
	for (j in 1:M) {
		T[k,j]=max(0,1-abs(data$x[k]-theta[j])/delta.x)+(2*kron(k,1)-kron(k,2))*max(0,1-abs(data$x[k]-theta[j]+delta.x)/delta.x)+(2*kron(k,M)-kron(k,M-1))*max(0,1-abs(data$x[k]-theta[j]-delta.x)/delta.x)
	}
}

dum=stats.bin(data$x,data$D,breaks=breaks)
for (j in 1:M) {
	if (j!=9) {
		for (i in -5:5) {
			a[j,i+6]=dum$stats["Q1",j]+i*max(data$sigma[which(data$bin==j)])/10
			a[9,i+6]=0
		}
	}
}

for (l in 1:M) {
	for (i in -5:5) {
		for (k in 1:N) {
			for (j in 1:M) {
				if (j==l) {F[l,i+6,k]=F[l,i+6,k]+T[k,j]*a[j,i+6]} else {F[l,i+6,k]=F[l,i+6,k]+T[k,j]*a[j,0+6]}
			}
		}
	}
}

Chisq=array(0.,c(M,11))
for (j in 1:M) {
	for (i in -5:5) {
		for (k in 1:N) {
			Chisq[j,i+6] = Chisq[j,i+6] + (F[j,i+6,k]-data$D[k])^2/(data$sigma[k])^2
		}
	}
}

af = rep(NA,M)
for (i in (1:M)){af[i] = a[which(Chisq==min(Chisq[,i]),arr.ind=TRUE)]}

#2
T = t(T)%*%T
dirty = function(i,j){
    tmp = 0
    for(k in (1:M)){
        tmp = tmp + (T[k,i] * T[k,j])/(data$sigma[k]**2) 
        }
         return(tmp)
        }
H=array(0.,c(M,M))
for (i in 1:M){
    for (j in 1:M){
        H[i,j] = 2 * dirty(i,j)}}

lambda=eigen(H)$val
e=eigen(H)$vec
l=seq(1,10,1)

x11()
plot(l,lambda, main="Eigenvalues", xlab="l", ylab="lambda")
x11()
plot(theta,e[1,], main = "Eigenvectors")
for (l in 2:M) {
	lines(theta+0.0001*l,e[l,],type="p",pch=l)
}
x11()
plot(theta,e[1,]/lambda[1], main = "Eigenvectors --- weighted")
for (l in 2:M) {
	lines(theta+0.0001*l,e[l,]/lambda[l],type="p",pch=l)
}

T_e = solve(e) %*% T %*%e

#########

Ff = function(T,a){
    tmp = 0
    for (j in (1:length(T))){
        tmp = tmp + T[j]*a[j]
        return(tmp)
    }}
afm = array(NA,c(M,M))
for (Mn in (1:M)){
    print('we are at')
    print(Mn)
    N=length(data$x)
    theta=rep(0,Mn)
    delta.x=(max(data$x)-min(data$x))/Mn
    for (j in 1:Mn){
        theta[j]=min(data$x)+delta.x*(j-0.5)
    }

kron = function(a,b) {
        if (a==b) {return(1.)} else {return(0.)}
    }

T=array(0.,c(N,Mn))
    a=array(0.,c(10,11))
    F=array(0.,c(11,11,N))
    breaks=seq(min(data$x)-1e-9,max(data$x)+1e-9,length.out=11)
    bin = findInterval(data$x,breaks)
    data = cbind(data,bin)

dum=stats.bin(data$x,data$D,breaks=breaks)
    for (j in 1:Mn) {
        if (j!=9) {
            for (i in -5:5) {
                a[j,i+6]=dum$stats["Q1",j]+i*max(data$sigma[which(data$bin==j)])/Mn
                a[9,i+6]=0
            }
        }
    }

for (k in 1:N) {
        for (j in 1:Mn) {
            T[k,j]=max(0,1-abs(data$x[k]-theta[j])/delta.x)+(2*kron(k,1)-kron(k,2))*max(0,1-abs(data$x[k]-theta[j]+delta.x)/delta.x)+(2*kron(k,M)-kron(k,M-1))*max(0,1-abs(data$x[k]-theta[j]-delta.x)/delta.x)
        }
    }
    T = t(T)%*%T
    dirty = function(i,j){
        tmp = 0
        for(k in (1:Mn)){
            tmp = tmp + (T[k,i] * T[k,j])/(data$sigma[k]**2) 
            }
             return(tmp)
            }
    H=array(0.,c(Mn,Mn))
    for (i in 1:Mn){
        for (j in 1:Mn){
            H[i,j] = 2 * dirty(i,j)}}
    lambda=eigen(H)$val
    e=eigen(H)$vec
    l=seq(1,10,1)

Chisq=array(0.,c(Mn,11))
    for (j in 1:Mn) {
        for (i in -5:5) {
            for (k in 1:Mn) {
                Chisq[j,i+6] = Chisq[j,i+6] + (Ff(T[k,],a)-data$D[k])^2/(data$sigma[k])^2
            }

}
    }
    for (i in (1:Mn)){afm[Mn,i] = a[which(Chisq==min(Chisq[,i]),arr.ind=TRUE)[Mn]]}
}
print('The best fit c')
print(afm)

SS2012-StatisticalMethods.git

add sheet 9