R 中简单优化的逻辑相同但结果不同
Same logic but different results from a simple optimization in R
我完全被以下简单的 R 代码搞糊涂了。在 第一部分 x
将等于 v
(这就是我想要的)。
但奇怪的是,在 第二部分 我更改了输入值,但遵循与第一部分完全相同的逻辑 但是这次 x
和 v
不再匹配!我很想知道问题出在哪里?
第一部分:
m1 = 5
m2 = 1.3*m1
A = m1 + m2
x = 5
a <- function(m3){
abs((m1 - (A + m3)/3)^2 + (1.3*m1 - (A + m3)/3)^2 + (m3 - (A + m3)/3)^2 - 3*x) }
m3 = optimize(a, interval = c(0, 100), tol = 1e-20)[[1]]
v = var(c(m1, m2, m3))*(2/3) # gives "5" same as "x"
第二部分:
eta.sq = .25
beta = qnorm(c(1e-12, .999999999999))
q = c(0, 25)
mu.sig = solve(cbind(1L, beta), q)
m1 = mu.sig[[1]]
H = (mu.sig[[2]])^2
m2 = 1.3 * m1
A = m1 + m2
x = (H * eta.sq) / (1 - eta.sq) # "x" is: 1.052529
a = function(m3){
abs((m1 - (A + m3)/3)^2 + (1.3*m1 - (A + m3)/3)^2 + (m3 - (A + m3)/3)^2 - 3*x) }
m3 = optimize(a, interval = c(0, 100), tol = 1e-20)[[1]]
v = var(c(m1, m2, m3))*(2/3) # "v" is: 2.343749
不同之处在于,对于您的第一部分,函数 a
有两个根,优化函数找到其中一个 (m3=10.31207
)。在 m3
的这个值上,a(m3)==0
意味着 m1
、m2
和 m3
的归一化平方和 (SS) 等于3*x
:
> a(m3)
[1] 3.348097e-07
> ss <- function(x) { sum((x-mean(x))^2) }
> ss(c(m1, m2, m3))
[1] 15
> 3*x
[1] 15
>
根据样本方差的定义,变量v
等于SS的三分之一,所以得到v==x
.
相比之下,在第二部分中,您的函数 a
没有根。它在 m3=14.375
处达到最小值,但在 m3
处,a(m3)==3.87366
的值不为零,因此归一化平方和不等于 3*x
,因此没有理由期望 v
(SS 的三分之一)应该等于 x
.
> a(m3)
[1] 3.87366
> ss(c(m1, m2, m3))
[1] 7.031247 -- actual SS value...
> 3*x
[1] 3.157587 -- ...couldn't be optimized to equal 3*x
>
我完全被以下简单的 R 代码搞糊涂了。在 第一部分 x
将等于 v
(这就是我想要的)。
但奇怪的是,在 第二部分 我更改了输入值,但遵循与第一部分完全相同的逻辑 但是这次 x
和 v
不再匹配!我很想知道问题出在哪里?
第一部分:
m1 = 5
m2 = 1.3*m1
A = m1 + m2
x = 5
a <- function(m3){
abs((m1 - (A + m3)/3)^2 + (1.3*m1 - (A + m3)/3)^2 + (m3 - (A + m3)/3)^2 - 3*x) }
m3 = optimize(a, interval = c(0, 100), tol = 1e-20)[[1]]
v = var(c(m1, m2, m3))*(2/3) # gives "5" same as "x"
第二部分:
eta.sq = .25
beta = qnorm(c(1e-12, .999999999999))
q = c(0, 25)
mu.sig = solve(cbind(1L, beta), q)
m1 = mu.sig[[1]]
H = (mu.sig[[2]])^2
m2 = 1.3 * m1
A = m1 + m2
x = (H * eta.sq) / (1 - eta.sq) # "x" is: 1.052529
a = function(m3){
abs((m1 - (A + m3)/3)^2 + (1.3*m1 - (A + m3)/3)^2 + (m3 - (A + m3)/3)^2 - 3*x) }
m3 = optimize(a, interval = c(0, 100), tol = 1e-20)[[1]]
v = var(c(m1, m2, m3))*(2/3) # "v" is: 2.343749
不同之处在于,对于您的第一部分,函数 a
有两个根,优化函数找到其中一个 (m3=10.31207
)。在 m3
的这个值上,a(m3)==0
意味着 m1
、m2
和 m3
的归一化平方和 (SS) 等于3*x
:
> a(m3)
[1] 3.348097e-07
> ss <- function(x) { sum((x-mean(x))^2) }
> ss(c(m1, m2, m3))
[1] 15
> 3*x
[1] 15
>
根据样本方差的定义,变量v
等于SS的三分之一,所以得到v==x
.
相比之下,在第二部分中,您的函数 a
没有根。它在 m3=14.375
处达到最小值,但在 m3
处,a(m3)==3.87366
的值不为零,因此归一化平方和不等于 3*x
,因此没有理由期望 v
(SS 的三分之一)应该等于 x
.
> a(m3)
[1] 3.87366
> ss(c(m1, m2, m3))
[1] 7.031247 -- actual SS value...
> 3*x
[1] 3.157587 -- ...couldn't be optimized to equal 3*x
>