git.schokokeks.org
Repositories
Help
Report an Issue
SS2012-StatisticalMethods.git
Code
Commits
Branches
Tags
Suche
Strukturansicht:
0d1a76b
Branches
Tags
master
SS2012-StatisticalMethods.git
9
solution9.R
add sheet 9
Michi
commited
0d1a76b
at 2012-07-01 21:25:27
solution9.R
Blame
History
Raw
#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)