2017-10-01 6 views
0

Je suis complètement déconcerté par le code R simple suivant. Dans la première partie sera égal à v (c'est ce que je veux).Même logique mais résultats différents d'une simple optimisation dans R

Mais étrangement dans la deuxième partie je change les valeurs d'entrée, mais suivre exactement la même logique que dans la première partie mais cette fois x et v ne correspondent plus! Je me demande profondément où est le problème?

Première partie:

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" 

Deuxième partie:

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 

Répondre

1

La différence est que pour votre première partie, la fonction a a deux racines, et la fonction optimize trouve un d'entre eux (m3=10.31207). A cette valeur de m3, le fait que a(m3)==0 implique que la somme normalisée des carrés (SS) de m1, m2 et m3 est égal à 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 
> 

Par la définition de la variance de l'échantillon, la v variables est égal à un tiers de la SS, de sorte que vous obtenez v==x. Par contre, dans la deuxième partie, votre fonction a n'a pas de racine. Il atteint un minimum à m3=14.375, mais à cette valeur de m3, la valeur de a(m3)==3.87366 n'est pas zéro, donc la somme normalisée des carrés n'est pas égale à 3*x, et donc il n'y a aucune raison de s'attendre à ce que v (un tiers de la SS) devrait être égal à 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 
>