Browse code

add sheet 9

Michi authored on01/07/2012 21:25:27
Showing1 changed files
1 1
new file mode 100644
... ...
@@ -0,0 +1,177 @@
1
+#Sheet 9 Michael Klauser
2
+library(gplots)
3
+library(fields)
4
+Filepath="ProjectedDensity.dat"
5
+data=scan(Filepath, what=list(x=0., D=0., sigma=0.))
6
+data=data.frame(x=data$x, D=data$D, sigma=data$sigma)
7
+
8
+plotCI(data$x, data$D, uiw=data$sigma, xlab="Angular Position", ylab="Projected Density")
9
+
10
+M=10
11
+N=length(data$x)
12
+theta=rep(0,M)
13
+delta.x=(max(data$x)-min(data$x))/M
14
+for (j in 1:M){
15
+	theta[j]=min(data$x)+delta.x*(j-0.5)
16
+}
17
+
18
+kron = function(a,b) {
19
+	if (a==b) {return(1.)} else {return(0.)}
20
+}
21
+
22
+T=array(0.,c(N,M))
23
+a=array(0.,c(10,11))
24
+F=array(0.,c(11,11,N))
25
+breaks=seq(min(data$x)-1e-9,max(data$x)+1e-9,length.out=11)
26
+bin = findInterval(data$x,breaks)
27
+data = cbind(data,bin)
28
+
29
+for (k in 1:N) {
30
+	for (j in 1:M) {
31
+		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)
32
+	}
33
+}
34
+
35
+dum=stats.bin(data$x,data$D,breaks=breaks)
36
+for (j in 1:M) {
37
+	if (j!=9) {
38
+		for (i in -5:5) {
39
+			a[j,i+6]=dum$stats["Q1",j]+i*max(data$sigma[which(data$bin==j)])/10
40
+			a[9,i+6]=0
41
+		}
42
+	}
43
+}
44
+
45
+for (l in 1:M) {
46
+	for (i in -5:5) {
47
+		for (k in 1:N) {
48
+			for (j in 1:M) {
49
+				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]}
50
+			}
51
+		}
52
+	}
53
+}
54
+
55
+Chisq=array(0.,c(M,11))
56
+for (j in 1:M) {
57
+	for (i in -5:5) {
58
+		for (k in 1:N) {
59
+			Chisq[j,i+6] = Chisq[j,i+6] + (F[j,i+6,k]-data$D[k])^2/(data$sigma[k])^2
60
+		}
61
+	}
62
+}
63
+
64
+af = rep(NA,M)
65
+for (i in (1:M)){af[i] = a[which(Chisq==min(Chisq[,i]),arr.ind=TRUE)]}
66
+
67
+#2
68
+T = t(T)%*%T
69
+dirty = function(i,j){
70
+    tmp = 0
71
+    for(k in (1:M)){
72
+        tmp = tmp + (T[k,i] * T[k,j])/(data$sigma[k]**2) 
73
+        }
74
+         return(tmp)
75
+        }
76
+H=array(0.,c(M,M))
77
+for (i in 1:M){
78
+    for (j in 1:M){
79
+        H[i,j] = 2 * dirty(i,j)}}
80
+
81
+
82
+lambda=eigen(H)$val
83
+e=eigen(H)$vec
84
+l=seq(1,10,1)
85
+
86
+x11()
87
+plot(l,lambda, main="Eigenvalues", xlab="l", ylab="lambda")
88
+x11()
89
+plot(theta,e[1,], main = "Eigenvectors")
90
+for (l in 2:M) {
91
+	lines(theta+0.0001*l,e[l,],type="p",pch=l)
92
+}
93
+x11()
94
+plot(theta,e[1,]/lambda[1], main = "Eigenvectors --- weighted")
95
+for (l in 2:M) {
96
+	lines(theta+0.0001*l,e[l,]/lambda[l],type="p",pch=l)
97
+}
98
+
99
+T_e = solve(e) %*% T %*%e
100
+
101
+
102
+#########
103
+
104
+Ff = function(T,a){
105
+    tmp = 0
106
+    for (j in (1:length(T))){
107
+        tmp = tmp + T[j]*a[j]
108
+        return(tmp)
109
+    }}
110
+afm = array(NA,c(M,M))
111
+for (Mn in (1:M)){
112
+    print('we are at')
113
+    print(Mn)
114
+    N=length(data$x)
115
+    theta=rep(0,Mn)
116
+    delta.x=(max(data$x)-min(data$x))/Mn
117
+    for (j in 1:Mn){
118
+        theta[j]=min(data$x)+delta.x*(j-0.5)
119
+    }
120
+
121
+    kron = function(a,b) {
122
+        if (a==b) {return(1.)} else {return(0.)}
123
+    }
124
+
125
+    T=array(0.,c(N,Mn))
126
+    a=array(0.,c(10,11))
127
+    F=array(0.,c(11,11,N))
128
+    breaks=seq(min(data$x)-1e-9,max(data$x)+1e-9,length.out=11)
129
+    bin = findInterval(data$x,breaks)
130
+    data = cbind(data,bin)
131
+
132
+
133
+
134
+    dum=stats.bin(data$x,data$D,breaks=breaks)
135
+    for (j in 1:Mn) {
136
+        if (j!=9) {
137
+            for (i in -5:5) {
138
+                a[j,i+6]=dum$stats["Q1",j]+i*max(data$sigma[which(data$bin==j)])/Mn
139
+                a[9,i+6]=0
140
+            }
141
+        }
142
+    }
143
+
144
+    for (k in 1:N) {
145
+        for (j in 1:Mn) {
146
+            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)
147
+        }
148
+    }
149
+    T = t(T)%*%T
150
+    dirty = function(i,j){
151
+        tmp = 0
152
+        for(k in (1:Mn)){
153
+            tmp = tmp + (T[k,i] * T[k,j])/(data$sigma[k]**2) 
154
+            }
155
+             return(tmp)
156
+            }
157
+    H=array(0.,c(Mn,Mn))
158
+    for (i in 1:Mn){
159
+        for (j in 1:Mn){
160
+            H[i,j] = 2 * dirty(i,j)}}
161
+    lambda=eigen(H)$val
162
+    e=eigen(H)$vec
163
+    l=seq(1,10,1)
164
+
165
+    Chisq=array(0.,c(Mn,11))
166
+    for (j in 1:Mn) {
167
+        for (i in -5:5) {
168
+            for (k in 1:Mn) {
169
+                Chisq[j,i+6] = Chisq[j,i+6] + (Ff(T[k,],a)-data$D[k])^2/(data$sigma[k])^2
170
+            }
171
+
172
+        }
173
+    }
174
+    for (i in (1:Mn)){afm[Mn,i] = a[which(Chisq==min(Chisq[,i]),arr.ind=TRUE)[Mn]]}
175
+}
176
+print('The best fit c')
177
+print(afm)