-
Notifications
You must be signed in to change notification settings - Fork 2
/
ModalRegression.R
48 lines (42 loc) · 1.03 KB
/
ModalRegression.R
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
## 1d modal regression
RegMS1d = function(X, Y, G.x=X, G.y=Y, h.x, h.y, iter=100, tolerance=1e-8){
x0 = G.x
y0 = G.y
#initialize points
m=length(y0)
#size of grid
i0=0
tol0=1
while(i0<iter && tol0> tolerance){
y0.tmp = y0
y0 = sapply(1:m, function(x){
y0.tmp2 = sum(Y*dnorm(y0[x]-Y, sd=h.y)*dnorm(x0[x]-X,sd=h.x))/sum(dnorm(y0[x]-Y, sd=h.y)*dnorm(x0[x]-X,sd=h.x))
return(y0.tmp2)
}
)
i0=i0+1
tol0 = max(abs(y0.tmp-y0))
} # Partial mean shift
return(y0)
}
### 2d modal regression
RegMS2d = function(X,Y, G.x=X, G.y=Y, h.x, h.y, iter=100,tolerance=1e-8){
x0 = G.x
y0 = G.y
#initialize points
m=length(y0)
#size of grid
i0=0
tol0=1
while(i0<iter && tol0> tolerance){
y0.tmp = y0
y0 = sapply(1:m, function(x){
y0.tmp2 = sum(Y*dnorm(y0[x]-Y, sd=h.y)*dnorm(x0[x,1]-X[,1],sd=h.x)*dnorm(x0[x,2]-X[,2],sd=h.x))/sum(dnorm(y0[x]-Y, sd=h.y)*dnorm(x0[x,1]-X[,1],sd=h.x)*dnorm(x0[x,2]-X[,2],sd=h.x))
return(y0.tmp2)
}
)
i0=i0+1
tol0 = max(abs(y0.tmp-y0))
} # Partial mean shift
return(y0)
}