CODEWITHSUNDEEP
rfunctionmathoptimizationalgebra
rnorouzian
rnorouzian

Reputation: 7517

Same logic but different results from a simple optimization in R

I'm completely baffled by the following simple R code. In the first part x will equal v (that's what I want).

But then strangely in the second part I change the input values but follow the exact same logic as in the first part HOWEVER this time x and v no longer match! I'm deeply wondering where is the problem?

First Part:

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"

Second Part:

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

Upvotes: 0

Views: 73

Answers (1)

K. A. Buhr
K. A. Buhr

Reputation: 50819

The difference is that for your first part, the function a has two roots, and the optimize function finds one of them (m3=10.31207). At this value of m3, the fact that a(m3)==0 implies that the normalized sum of squares (SS) of m1, m2, and m3 is equal to 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
>

By the definition of the sample variance, the variable v is equal to one-third the SS, so you get v==x.

In contrast, in the second part, your function a has no roots. It attains a minimum at m3=14.375, but at this value of m3, the value of a(m3)==3.87366 is not zero, so the normalized sum of squares is not equal to 3*x, and so there's no reason to expect that v (one-third the SS) should equal 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
>

Upvotes: 1

Related Questions