Playground
No Noise
n <- 100
x <- seq(0, 1, length.out=n)
b <- 0.2
c <- 0.8
A <- rep(0,n)
B <- ifelse(x<b, 0, ifelse(x > c, 0, 1))
l1 <- c(0,vapply(x[2:n], function(i) {
(mean(B[x < i]) - mean(A[x < i]))^2 * sum(x < i)
}, FUN.VALUE=numeric(1)))
r1 <- vapply(x, function(i) {
(mean(B[x >= i]) - mean(A[x >= i]))^2 * sum(x >= i)
}, FUN.VALUE=numeric(1))
l2 <- c(0,vapply(x[2:n], function(i) {
(mean(B[x < i]) - mean(A[x < i]))
}, FUN.VALUE=numeric(1)))
r2 <- vapply(x, function(i) {
(mean(B[x >= i]) - mean(A[x >= i]))
}, FUN.VALUE=numeric(1))par(mfrow=c(2,2))
plot(x, l1+r1, type="l")
plot(x, abs(l1-r1), type="l")
plot(x, l2+r2, type="l")
plot(x, abs(l2-r2), type="l")
Random Noise
n <- 100
x <- seq(0, 1, length.out=n)
b <- 0.2
c <- 0.8
A <- rep(0,n) + rnorm(n)
B <- ifelse(x<b, 0, ifelse(x > c, 0, 1)) + rnorm(n)
l1 <- c(0,vapply(x[2:n], function(i) {
(mean(B[x < i]) - mean(A[x < i]))^2 * sum(x < i)
}, FUN.VALUE=numeric(1)))
r1 <- vapply(x, function(i) {
(mean(B[x >= i]) - mean(A[x >= i]))^2 * sum(x >= i)
}, FUN.VALUE=numeric(1))
l2 <- c(0,vapply(x[2:n], function(i) {
(mean(B[x < i]) - mean(A[x < i]))
}, FUN.VALUE=numeric(1)))
r2 <- vapply(x, function(i) {
(mean(B[x >= i]) - mean(A[x >= i]))
}, FUN.VALUE=numeric(1))par(mfrow=c(2,2))
plot(x, l1+r1, type="l")
plot(x, abs(l1-r1), type="l")
plot(x, l2+r2, type="l")
plot(x, abs(l2-r2), type="l")
Iterations
Nsim <- 10^3
n <- 101
b <- 0.4
c <- 0.6
x <- seq(0, 1, length.out=n)
A0 <- rep(0.5, n)
B0 <- ifelse(x < b, 0, ifelse(x > c, 2, 1))
calc <- function(f, A, B, L, R) {
s <- vector(,R-L+1)
s[1] <- -Inf # sums are like < i + >= i, so first is all
for (i in (L+1):R) {
s[i-L+1] <- f(L:(i-1), i:R, A, B)
}
s
}
f1 <- function(L.ind, R.ind, A, B) {
f11(L.ind, A, B) + f11(R.ind, A, B)
}
f11 <- function(ind, A, B) {
(mean(B[ind]) - mean(A[ind]))^2 * length(ind)
}
f2 <- function(L.ind, R.ind, A, B) {
abs(f22(L.ind, A, B) - f22(R.ind, A, B))
}
f22 <- function(ind, A, B) {
(mean(B[ind]) - mean(A[ind]))^2 * length(ind)
}
f3 <- function(L.ind, R.ind, A, B) {
abs(f33(L.ind, A, B) - f33(R.ind, A, B))
}
f33 <- function(ind, A, B) {
(mean(B[ind]) - mean(A[ind]))
}
m <- matrix(, nrow=Nsim, ncol=6)
mm <- matrix(, nrow=Nsim, ncol=6)
for (i in 1:Nsim) {
set.seed(i)
A <- A0 + rnorm(n)
B <- B0 + rnorm(n)
r1 <- calc(f1, A, B, 1, length(x))
r2 <- calc(f2, A, B, 1, length(x))
r3 <- calc(f3, A, B, 1, length(x))
a1 <- which.max(r1)
a2 <- which.max(r2)
a3 <- which.max(r3)
b1 <- ifelse(a1 > 50, which.max(calc(f1, A, B, 1, a1-1)), a1 + which.max(calc(f1, A, B, a1+1, length(x))))
b2 <- ifelse(a2 > 50, which.max(calc(f2, A, B, 1, a2-1)), a2 + which.max(calc(f2, A, B, a2+1, length(x))))
b3 <- ifelse(a3 > 50, which.max(calc(f3, A, B, 1, a3-1)), a3 + which.max(calc(f3, A, B, a3+1, length(x))))
m[i,] <- c(sort(c(a1, b1)), sort(c(a2, b2)), sort(c(a3, b3)))
# mm[i,] <- c(a1, a2, b1, b2)
# print(m[i,])
# par(mfrow=c(1,2))
# print(plot(r1, type="l"))
# print(plot(r2, type="l"))
}par(mfrow=c(3,2))
hist(m[,1])
hist(m[,2])
hist(m[,3])
hist(m[,4])
hist(m[,5])
hist(m[,6])
sum((m[,1] - 41)^2) + sum((m[,2] - 62)^2)[1] 358063sum((m[,3] - 41)^2) + sum((m[,4] - 62)^2)[1] 453308sum((m[,5] - 41)^2) + sum((m[,6] - 62)^2)[1] 2370624# sum((mm[,1] > 50)*(mm[,1]-82)^2+(mm[,1] <= 50)*(mm[,1]-21)^2)
# sum((mm[,2] > 50)*(mm[,2]-82)^2+(mm[,2] <= 50)*(mm[,2]-21)^2)