总共有四个文件用于距包估计

  1. source("~/R/optim/test_optim/test_problem_2_quadratic_on_sphere.R") 是二维中最好用估计距包的版本，想从具体问题出发构建距包估计。
  2. /home/ghy/R/optim/matrix_hull.R 是我要解决多维多维问题的版本
  3. /home/ghy/R/matrix_hull.R 是最老的版本
  4. ~/R/optim/estimation_of_conditional_ess_inf.R 刚开始专注本质下却界估计的版本，包含距包和凸包

凸包距包的定义

数学严格定义，集合S的距包为所有象限集集合的交集，即 \[Mhull(S) = \bigcap \{A: S \subseteq A^c,I_A = I(x <c)I(x>d)I(y<s)I(y>t),\; c,d,s,t \in R \cup\{\infty, -\infty\}\}\] 对比凸包的严格定义，我们知道： \[Chull(S) = \bigcap \{A: S \subseteq A^c, I_A = I(ax + by <c)\; or\; I(ax + by >c)\; or\; I(ay + by<c)\;or \; I(ay +by >c),\; a,b,c \in R\}\]

凸包，距包做为一种点集合的包面，它的核心新知在于任意给定一个样本点，它一定在包面的内部，并且过任何一个样本点的一维直线（所谓一维直线是指只有一个坐标轴可以变动的曲线）在包面内部是联通集。

上面定义的方法可以直接推广到多维的情况。

上面数学定义中的距包是直线距包，普通距包是把这直线距包的极点用直线连接起来的封闭曲线。

凸包，距包的图示

图示就是距包，凸包的简单展示。

source("~/R/optim/test_optim/test_problem_2_quadratic_on_sphere.R")
source("/home/ghy/R/optim/matrix_hull.R")
set.seed(1)
dat = matrix(rnorm(n*m, sd = sigma), ncol = n) %>%
    (function(dat){cbind(g_value = apply(dat, 1, g),f_value =  apply(dat, 1, f),  dat)}) %>%
    (function(dat){colnames(dat) = c("g_value", "f_value", paste0("x", 1:n)); dat}) %>% 
    as.data.frame()
d = dat
matrix.hull.plot(dat);vert.matrix.hull.plot(dat);vert.matrix.hull.plot2(dat)

double.2d.quick.matrix.hull.plot(dat)

ind1 = matrix.hull.index(dat)
dat1 = vert.matrix.hull(dat[-ind1,])

vert.matrix.hull.plot(d) + 
  geom_path(aes(x = g_value,y = f_value),data = dat1, color = "blue")

新的算法，专业用处理距包，构建算法之后，我们对比原来算法

多维凸包和距包

如果能搞定三维的推广，那么高维的问题能解决。

我们的方法具有明显的优势，我们的方法并不依赖于数据的维数。

下面的例子告诉我，我的凸包下部函数可以写得更简单，直接从第一个点开始取就行了。

球面上的每个点都是凸包的极点

下面的代码画出了三维的凸包！

## 自带的包只能处理2维
m = 100
dat = matrix(rnorm(m*3), ncol = 3)
geometry::convhulln(dat) %>% head
##      [,1] [,2] [,3]
## [1,]   36   47   30
## [2,]   37   47   30
## [3,]   74   36    2
## [4,]   74   36   47
## [5,]   74   88    2
## [6,]   74   88   47
hpts <- chull(dat)
plot(x = dat[,1], y = dat[,2])
lines(dat[hpts, ])

##  geometry包能够处理多维的凸包
library(geometry)
library(rgl)
## 生成球面上的点 
ps <- matrix(rnorm(300), ncol=3)  # generate points on a sphere
ps <- sqrt(3)*ps/drop(sqrt((ps^2) %*% rep(1, 3)))

ts.surf <- t(convhulln(ps))  # see the qhull documentations for the options
## Not run: 
# rgl.triangles(ps[ts.surf,1],ps[ts.surf,2],ps[ts.surf,3],col="red",alpha=.5)
# for(i in 1:(8*36)) rgl.viewpoint(i/8)
# convhulln(ps, options = "Fv")

## 得到极点点集合
## 从中可以看出随着维度的增加，我们的方法具有明显的优势，我们的方法并不依赖于数据的维数。
m = 100
dat = matrix(rnorm(m*6), ncol = 20)
convhulln(dat) %>% as.vector() %>% unique %>% sort
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
## [24] 24 25 26 27 28 29 30
small.m.matrix.hull2(dat) %>% convhulln() %>% as.vector() %>% unique %>% sort
##  [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
## [24] 24 25 26 27 28 29 30

多维距包的算法

二维的距包算法速度很快了，三维或者多维距包我们才用了另外一个算法，速度有一个明显的降低。而先取子集抽样的方法能够大大的降低计算时间。

1. 找一个点作为坐标原点
2. 判断每个点属于哪个象限
3. 选定初始集合，每次增加一个点， 根据原来的距包判断是否讲它包含进来
4. 穷尽所有的点

多维距包的基础函数为 points.inside.matrix.hull.index, 它的最大算法复杂度是 \(n^2\) 比凸包 \(nlog(n) + n^{[\frac{d}{2}]}\)要好。通过一定的处理能把算法复杂度处理成\(n log(h)\), where h is the num of 距包极点个数。

每次随机抽取少量的点，求得多个距包，然后用距包支撑点在做距包！

source("~/R/optim/matrix_hull.R")
set.seed(10)
dat = matrix(rnorm(50*2), ncol = 2)
hull = small.m.matrix.hull(dat) %>% unique
qplot(x = dat[,2], y = dat[,1], geom = "point") + 
  geom_point(mapping = aes(x = x1, y = y), data = hull, size = 8, alpha = 0.2, color = "red")

## 更快的算法  
set.seed(10)
dat = matrix(rnorm(200*2), ncol = 2)
hull = dat[!apply(dat, 1, function(x) if.point.inside.hull(hull = dat, x)),]
qplot(dat[,1], dat[,2], geom = "point") + geom_point(mapping = aes(hull[,1], hull[,2]), color = "red", size = 5, alpha = 0.6)

source("~/R/optim/test_optim/test_problem_2_quadratic_on_sphere.R")
source("~/R/optim/matrix_hull.R")
set.seed(1)
dat = matrix(rnorm(30*2), ncol = 2)
hull = matrix.hull(dat)
# matrix.hull.index(dat)
p = (dat[1:5,] + 0.5) %>% as.data.frame() %>% 
  (function(dat) {names(dat) = c("x", "y"); dat}) %>%
  mutate(z = apply(data.frame(x,y), 1, function(x) if.point.inside.hull(hull, x)))

vert.matrix.hull.plot(dat) + geom_point(aes(x , y, color = z),data = p, size = 3, alpha = 0.8)

# apply(p, 1, function(p) point.inside.hull(hull, p))
# apply(p, 1, function(p) if.point.inside.hull(hull, p))

用凸包作为初始值，然后不停的加点.

If \(p\) is inside the hull, 存在一个同象限的点\(P_1\)包含它。意味着存在一个点 \(P_1\), satisfies that \(P_1 - P >= 0\), 即 \[|sign(P) - sign(P_1)|\leq 1\; \& \; sign(P_1-P) \geq \bf{0}\] 按照某个坐标，将点排序，然后按照顺序添加点。直接判断p是否属于每个点的象限矩形————这样做速度慢，需要调整，下面是一种调整

**一个十分重要的信息，距包是可以用最快的算法的，算法复杂度只有n, 每个点可以直接判断是否属于距包内部。

source("~/R/optim/matrix_hull.R")
set.seed(10)
m = 200; n = 2
dat = matrix(rnorm(m*n), ncol = n)
hull = convex.to.find.matrix.approx.hull(dat)
# small.m.matrix.hull(dat)
# small.m.matrix.hull2(dat)
convex.to.find.matrix.approx.hull.plot(dat)

# p = dat[10,];if.point.inside.hull(hull, p)
d = dat[apply(dat, 1, function(x) if.point.inside.hull(hull = dat, x)),] 
qplot(x = dat[,1], y = dat[,2], geom = "point") + 
  geom_point(aes(d[,1], d[,2]), size = 6, color = "red", alpha = 0.5) +
  geom_point(aes(hull[,1], hull[,2]), size = 3, color = "green", alpha = 0.5)