Introdução

Motivação

Um grande número de problemas complexos no mundo real pode ser encarado como um problema de otimização, em que o objetivo é maximizar ou minimizar um determinado objetivo. Ao invés de tomar decisões com base na experiência e intuição humana, há uma tendência cada vez maior de adotar ferramentas computacionais, baseadas em métodos de otimização, na análise de dados do mundo real a fim de tomar decisões com uma quantidade maior de informações.

Otimização é um tópico central da área de Pesquisa Operacional, na qual foram desenvolvidas várias técnicas clássicas de otimização, como a programação linear (proposta na década de 1940) e branch and bound (desenvolvida na década de 1960) (Schrijver 1998). Nas últimas décadas, houve o surgimento de novos algoritmos de otimização, muitas vezes denominados “otimização moderna” (Michalewicz et al. 2006), “heurísticas modernas” (Michalewicz e Fogel 2013) ou “metaheurísticas” (Luke 2012). Como em Cortez (2014), será adotado o primeiro termo, otimização moderna, para descrever esses algoritmos.

Em contraste com os métodos clássicos, os métodos de otimização modernos são, em geral, generalistas, ou seja, aplicáveis a uma ampla gama de problemas, dado ser necessário pouco conhecimento do domínio. Por exemplo, o problema de otimização não precisa ser diferenciável, o que é exigido por métodos clássicos como o do Método do Gradiente. Ao adotar métodos heurísticos modernos, o usuário precisa especificar apenas duas questões principais: a representação da solução, que define o espaço de busca e seu tamanho e a função objetivo, que define o quão boa é uma determinada solução, permitindo comparar diferentes soluções. Em particular, os métodos modernos são úteis para resolver problemas complexos, para os quais não tenho sido desenvolvido nenhum algoritmo de otimização epecializado. Por outro lado, ao contrário dos métodos clássicos, a os métodos modernos de otimização não garantem que a solução ideal será sempre encontrada. Entretanto, muitas vezes os métodos modernos alcançam soluções de alta qualidade com um uso muito mais razoável de recursos computacionais, como por exemplo, memória e esforço de processamento. (Michalewicz e Fogel 2013).

Representação de uma solução:

Cortez (2014) salienta que, quando são usados métodos modernos de otimização. é importante a decisão sobre a representação de possíveis soluçõe. Tal decisão define o espaço de busca e seu tamanho e influencia a forma de busca de novas soluções. Há várias possibilidades para representar uma solução. Números binários, inteiros, caracteres, valores reais e vetores ordenados, matrizes, árvores, além de, virtualmente, qualquer forma de representação baseada em computador (por exemplo, programa de computador). Uma determinada representação pode incluir uma mistura de codificações diferentes (por exemplo, valores binários e reais). Além disso, uma representação pode ser fixa (por exemplo, vetores binários fixos) ou de comprimento variável (por exemplo, árvores). Alguns desses tipos de representação são vinculados historicamente a métodos de otimização específicos. Por exemplo, codificações binárias são a base dos algoritmos genéticos (Holland 1992). A Busca Tabu foi projetada para funcionar em espaços discretos (Glover e Laguna, 1998). As codificações de valor real são adotadas por vários algoritmos evolutivos (por exemplo, estratégia de evolução) (Bäck e Schwefel 1993), evolução diferencial e enxames de partículas. A programação genética é usada para otimizar estruturas de árvores (Banzhaf et al. 1998). Salienta-se que muitas vezes esses métodos de otimização podem ser adaptados a outras representações. Por exemplo, um novo enxame de partículas foi proposto para otimização discreta em Chen et al. (2010).

Neste material, as representações adotadas estão restritas àquelas com implementações disponíveis nos pacotes do R utilizados.

Função objetivo:

Outro passo importante em otimização é a definição da função objetivo, a qual deve representar o objetivo (ou objetivos) que se deseja seja maximizado ou minimizado. Essa função permite comparar diferentes soluções, fornecendo uma classificação (função objetivo ordinal) ou um score da medida de qualidade (função numérica) (Michalewicz et al., 2006). No caso das funções numéricas, a forma pode ser convexa ou não convexa, com vários máximos (ou mínimos) locais (Fig. 1.1). A otimização de funções objetivo convexas é muito mais fácil de resolver, existindo algoritmos especializados (por exemplo, mínimos quadrados, programação linear) que são bastante eficazes para lidar com tais problemas (Boyd e Vandenberghe 2004). No entanto, muitos problemas práticos não são convexos, muitas vezes incluindo ruído e funções complexas, com descontinuidades. Os problemas de otimização podem também ser dinâmicos, mudando com o tempo. Para todos esses problemas complexos, uma alternativa interessante é usar algoritmos de otimização modernos que pesquisam apenas um subconjunto do espaço de pesquisa, mas com tendência a alcançar soluções quase ótimas em um espaço de tempo razoável.

Por padrão, algumas implementações de métodos de otimização executam apenas uma minimização de uma função objetivo numérica. Nesses casos, uma abordagem simples é transformar a função de maximização \(\max\{f(s)\}\) em uma minimização equivalente \(\min\{f^*(s)\}\), em que \(f^*(s) = f(s)\) e \(s\) denota a solução.

Em vários campos de aplicação existe mais de um objetivo para ser otimizado. Freqüentemente, esses objetivos são conflitantes e é necessário definir compensações, uma vez que a otimização de soluções sob um único objetivo pode levar a resultados inaceitáveis em relação aos objetivos restantes. Nesses casos, uma abordagem muito melhor é adotar uma otimização multiobjetivo. Nesse material dedicaremos atenção às otimizações mono-objetivas.

Restrições:

Existem dois tipos principais de restrições (Michalewicz 2008): hard e soft. As restrições hard não podem ser violadas e devem-se a fatores não controlados, tais como leis ou restrições físicas. As restrições soft estão relacionadas a outros objetivos do usuário (geralmente não prioritários), tais como aumentar a eficiência da produção e, ao mesmo tempo, reduzir os custos ambientais. Restrições soft podem ser tratadas através da adoção de uma abordagem multi-objetivo, enquanto as restrições hard podem originar soluções inviáveis que precisam ser tratadas pelo procedimento de otimização. Para lidar com soluções inviáveis. Nesse caso, podem ser adotados vários métodos (Michalewicz et al. 2006): pena de morte, penalização ponderada, reparação e geração de apenas soluções viáveis. A pena de morte é um método simples, que envolve atribuir um valor de penalidade muito grande, de forma que soluções inviáveis sejam rapidamente descartadas pela pesquisa. No entanto, esse método não é muito eficiente e muitas vezes coloca o esforço da busca no descarte de soluções e não em encontrar o valor ideal. Penalização ponderada é uma solução melhor, sendo também fácil de ser implementada. Por exemplo, muitas vezes, a forma de uma função de avaliação pode ser definida na forma \(f(s) = \text{objetivo}(s) - \text{penalidade}(s)\) (Rocha et al. 2001). Por exemplo, para um determinado negócio, uma possível função de avaliação poderia ser \(f = w_1 \times \text{Lucro} - w_2 \times \text{Custo}\). O principal problema com a penalidade ponderada é que muitas vezes é difícil encontrar os pesos ideais, em particular quando há várias restrições envolvidas. A abordagem de reparação transforma uma solução inviável em viável. Freqüentemente, isso é conseguido usando informações dependentes de domínio ou aplicando uma pesquisa local (por exemplo, procurando uma solução viável na vizinhança do espaço de solução). Por fim, as abordagens que apenas geram soluções viáveis são baseadas em decodificadores e operadores especiais. Os decodificadores funcionam apenas em um espaço de busca viável, adotando uma representação indireta, enquanto os operadores especiais usam o conhecimento do domínio para criar novas soluções a partir dos anteriores.

Métodos de otimização:

Há diferentes maneiras para classificar os métodos de otimização. Em Cortez (2014) são adotados três fatores de análise: o tipo de busca guiada; se a busca é determinística ou estocástica e se o método é inspirado por processos físicos ou biológicos. O tipo de busca pode ser cega ou orientada. O primeiro é exaustivo e pressupõe o esgotamento de todas as alternativas para encontrar a solução ótima, enquanto o último usa buscas anteriores para orientar a busca atual. Os métodos modernos usam uma busca guiada, que muitas vezes é subdividida em duas categorias principais: local, que pesquisa na vizinhança de uma solução inicial, e busca global, que usa uma população de soluções. Na maioria dos problemas práticos, com espaços de busca de alta dimensão, a busca exaustiva (cega pura) não é viável, exigindo muito esforço computacional. Por sua vez, a busca local (ou de único estado) apresenta, em geral, uma convergência muito mais rápida, quando comparada com os métodos de busca global. No entanto, se o espaço de busca for muito complexa ou com muito ruído, com vários mínimos locais, em que os métodos locais podem travar com muita facilidade, como por exemplo na figura abaixo:

Função objetivo côncava

Função objetivo côncava

Nesses casos, pode-se executar o algoritmo várias vezes, com diferentes soluções iniciais, para melhorar a convergência. Embora os algoritmos baseados em população tendam a exigir mais esforço computacional que métodos locais, eles realizam buscas simultâneas em regiões distintas do espaço de busca, funcionando melhor como métodos de otimização global.

Os diferentes tipos de busca podem também ser combinados. Por exemplo, pode ser definida uma busca em duas fases, em que emprega-se, por exemplo, um método global em uma primeira etapa, para identificar rapidamente regiões interessantes no espaço de busca e, em seguida, aprimoram-se as melhores soluções, empregando um método de pesquisa local. Há muitas outras maneiras de combinar tipos de busca.

Vários métodos de otimização modernos empregam algum grau de aleatoriedade, pertencendo, portanto, à família dos métodos de otimização estocástica, como o simulated annealing e os algoritmos genéticos (Luke 2012). Além disso, vários desses métodos são inspirados na natureza (por exemplo, algoritmos genéticos, otimização por enxame de partículas, colônia de formigas) (Holland 1975; Eberhart et al. 2001). A abaixo mostra a taxonomia completa dos métodos de otimização (com os respectivos pacotes R) apresentados em Cortez (2014) .

Função objetivo côncava

Função objetivo côncava

Os métodos de otimização modernos compartilham algumas características comuns. O algoritmo abaixo mostra (em pseudocódigo) um método de otimização moderno genérico que é aplicável a todos os métodos discutidos em Cortez (2014). Este algoritmo global recebe duas entradas, a função objetivo (\(f\)) e um conjunto de parâmetros de controle (\(C\)), que inclui não apenas os parâmetros internos do método (por exemplo, temperatura inicial, tamanho da população), mas também está relacionado com a representação da solução (por exemplo, limites inferior e superior, tipo de representação e comprimento). Em todos os métodos de otimização modernos, há uma configuração inicial seguida por um loop principal que termina quando um determinado critério de parada é atendido. Podem ser adotados ou combinados critérios distintos:

PSEUDO CÓDIGO

• Medidas computacionais máximas (Exemplo: número de iterações, cálculos de função de avaliação, tempo decorrido, etc).

• medidas-alvo (Exemplo: critério de parada se o valor ótimo for maior ou igual a um determinado limite).

• medidas de convergência (Exemplo: número de iterações sem qualquer melhoria na solução ou melhoria média atingida nas últimas iterações).

• medidas de distribuição (Exemplo: medir como as últimas soluções testadas estão espalhadas no espaço de busca e parar se a dispersão for menor que um determinado limite).

O que distingue os métodos está relacionado a dois aspectos principais. Primeiro, se em cada iteração há um único estado (com base local) ou uma população de soluções. Em segundo lugar, a forma como novas soluções são criadas (mudança de função) e utilizadas para guiar na busca (seleção de função). No pseudocódigo genérico, o número de iterações (i) também é incluído como uma entrada das funções change, best e select porque se presume que o comportamento dessas funções pode ser dinâmico, mudando conforme a evolução do método de busca.

Alguns problemas:

Uso da função optim(){stats}:

Considere uma função objetivo \(f(\mathbf{x})\), em que \(\mathbf{x}\) é um vetor. Problemas de otimização referem-se à busca de \(\mathbf{x}^*\), tal que \(f(\mathbf{x}^*)\) é um máximo (ou mínimo) local. No caso de maximização, tem-se então que:

\[ \mathbf{x}^* = \arg \max f(\mathbf{x}) \]

Grande parte dos problemas estatísticos de estimação são problemas de maximização. Por exemplo, se \(f\) é a função de verossimilhança e \(\mathbf{x}\) é o vetor cos valores do parâmtero, então \(\mathbf{x}^*\) é o estimador de máxima verossimilhança (EMV), que tem importantes propriedades assintóticas.

Nosso foco nesta seção será o uso da função optim() em problemas de minimização. Se você quiser maximizar, basta multiplicar a função objetivo por \(-1\). O método default para a função optim() é a rotina de otimização denominada algoritmo simplex Nelder-Mead (não usa derivadas). A sintaxe básica é optim(init, f), em que init é o vetor com os valores inciais que você deve especificar e f é a função objetivo. O argumento method apresenta as seguintes alternativas:

A função optim() pode também usar a derivada, armazenanda na função df. A sintaxe então é optim(init, f, df, method = "CG").

Salienta-se que os valores iniciais podem ter um grande impacto no ponto de convergência de algumas funções, principalmente aquelas com muitos mínimos,

Exemplos unidimensionais:

Exemplo 1

Suponha \(f(x) = \operatorname{e}^{-(x-2)^2}\), cuja derivada é \(f'(x) = -2(x-2)f(x)\). Como $f(x) é estritamente positivo, \(f'(x) = 0\) para \(x=2\). Assim, \(x^* = 2\) é um ponto de máximo. Como desejamos encontrar o máximo por meio da função optim(), precisamos usar \(-f(x)\) como função objetivo, assim:

# como o objetivo é maximizar, usa-se -f(x).
f <- function(x) -exp(-( (x-2)^2 ))

cujo gráfico pode ser observado abaixo:

Usando a função optim(), com valor incial e método de otimização default temos:

  # sem a derivada, usando valor inicial 1
optim(1, f)
## Warning in optim(1, f): one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## $par
## [1] 1.9
## 
## $value
## [1] -0.9900498
## 
## $counts
## function gradient 
##       10       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

O ponto \(x^*\) encontrado pelo algoritmo é 1.9 que está disponível em optim(1, f)$par. Por outro lado, -optim(1, f)$value fornece o valor do máximo (\(f(1.9)\)). Sabe-se entretanto que o ponto de máximo exato é \(x^* = 2\). Há também uma mensagem de alerta, indicando que é necessário o uso de método unidimensional ("Brent"). Esse método exige a delimitação do espaço de busca, que, inspecionando o gráfico, escolhemos como sendo o intervalo \((0, 4)\):

 # método "Brent", com busca ente 0 e 4
optim(1, f, method = "Brent", lower = 0, upper = 4)$par
## [1] 2

atingindo agora o valor exato de \(x^*=2\), com o máximo da função sendo \(f(2) = 1\). Podemos usar a derivada da função objetivo (df)

  # usando a derivada da função objetivo
df <- function(x) -2*(x-2)*f(x)
  # método"CG"
optim(1, f, df, method = "CG")$par
## [1] 2

Note que o método default ("Nelder-Mead") não utiliza a função da derivada e aparenta ser mais sensível ao valor inicial.

Exemplo 2:

Suponha a função \(f(x) = sen(cos(x))\).

f <- function(x) sin(x*cos(x))

Essa função tem muitos ótimos locais (mínimos e máximos). Observe o gráfico abaixo:

Vamos verificar como a função optim() encontra os mínimos dessa função que aparentam estar nas vizinhanças de 2.1, 4.1, 5.8, 7.5 e muitos outros. Observe abaixo os resultados obtidos para quatro buscas com valores iniciais distintos:

 # valor incial 2
a1 <- optim(2, f)$par
a1
## [1] 2.316016
  # valor inicial 4
a2 <- optim(4, f)$par
a2
## [1] 4.342236
  #valor incial 6
a3 <- optim(6, f)$par
a3
## [1] 5.688647
  # valor incial 8
a4 <- optim(8, f)$pa
a4
## [1] 7.132227

Os valores da função objetivo nesses pontos de mínimos locais são:

f(c(a1, a2, a3, a4))
## [1] -1 -1 -1 -1

Há vários pontos com a função objetivo atingindo o mesmo valor mínimo e aparenta que o algoritmo encontra com bastante precisão o mínimo mais próximo ao valor incial utilizado.

Exemplos bidimensionais:

Exemplo 3:

Seja a função \(f(x,y) = (1-x)^2 + 100(y - x^2)^2\), que é denominada função de Rosenbrock. Observe abaixo o gráfico dessa função:

  # função objetivo
f <- function(x1, y1) (1 - x1)^2 + 100*(y1 - x1^2)^2

x <- seq(-2, 2, by = 0.15)
y <- seq(-1, 3, by = 0.15)
z <- outer(x, y, f)

  # gráfico da superfície
persp(x, y, z, phi = 45, theta = -45, col = "yellow", shade = 0.00000001, ticktype = "detailed")

Essa função é estritamente positiva e é zero quando \(y = x^2\) e \(x = 1\), portanto \((1, 1)\) é um mínimo. Vamos verificar o uso de optim() para obter esse ponto de mínimo. É importante salientar qie seu uso em otimização multidimensional exige que a entrada da função objetivo seja um único vetor, ou seja:

  # função objetivo
f <- function(x) (1 - x[1])^2 + 100*(x[2] - x[1]^2)^2

implicando assim o uso de vetor como valor inicial, adotaremos \((0, 0)\) no presente caso:

  # valor inicial deve ser um vetor
optim( c(0,0), f )$par
## [1] 0.9999564 0.9999085

Obtém-se coordenadas bastante próximas daquelas do valor exato do mínimo \((1, 1)\).

Exemplo 4:

Seja a função \(f(x, y) = (x^2 + y - 11)^2 + (x + y^2 - 7)^2\), denominada função de Himmelblau. Observe abaixo o gráfico dessa função para \(x \in (-4, 4)\):

 # função objetivo
f <- function(x1, y1) (x1^2 + y1 - 11)^2 + (x1 + y1^2 - 7)^2

  # coordenadas da superfície
x <- seq(-4.5, 4.5, by = 0.2)
y <- seq(-4.5, 4.5, by = 0.2)
z <- outer(x, y, f)

  # gráfico da superfície
persp(x, y, z, phi = -45, theta = 45, col = "yellow", shade = 0.65, 
    ticktype = "detailed", cex.axis = 0.75)

Inspecionando o gráfico, percebem-se quatro “buracos” com aparência de serem valores mínimos, nas vizinhanças dos pontos \((-4,-4)\), \((2,-2)\), \((2, 2)\) e \((-4, 4)\).

Adotaremos como valores iniciais da função optim() os pares ordenados daqueles pontos observados na inspeção do gráfico da superfície, ou sejam, \((-4,-4)\), \((2,-2)\), \((2, 2)\) e \((-4, 4)\):

  # função objetivo (para uso no optim())
f <- function(x) (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2

  # Mínimo nas proximidades de (-4, -4)
ot1 <- optim(c(-4, -4), f)$par
ot1
## [1] -3.779347 -3.283172
  # Mínimo nas proximidades de (2, -2)
ot2 <- optim(c(2, -2), f)$par
ot2
## [1]  3.584370 -1.848105
  # Mínimo nas proximidades de (2, 2)
ot3 <- optim(c(2, 2), f)$par
ot3
## [1] 3.000014 2.000032
  # Mínimo nas proximidades de (-4, 4)
ot4 <- optim(c(-4, 4), f)$par
ot4
## [1] -2.805129  3.131435

Percebe-se que a função objetivo nesses pontos é aproximadamente zero, ou seja, eles são na realidade os pontos de mínimo. :

  # valores da função objetivo bidimensional nos pontos de mínimo
f(ot1); f(ot2); f(ot3); f(ot4)
## [1] 1.00504e-07
## [1] 1.787419e-07
## [1] 3.37533e-08
## [1] 6.008897e-07

A função objetiva é estritamente positiva e atinge seus pontos de mínimo quando \(x^2 + y - 11 = 0\) e \(x^2 + y^2 - 7 = 0\). Podendo-se verificar que nesses pontos as derivadas parciais são aproximadamente zero:

  # derivada da função objetivo
dfx <- function(x) x[1]^2 + x[2] - 11
dfy <- function(x) x[1] + x[2]^2 - 7

  # valores das derivadas parciais nos pontos encontrados
oti <- rbind(ot1, ot2, ot3, ot4)
for(i in 1:4) print(c(dfx(oti[i, ]), dfy(oti[i, ])))
## [1]  0.0002890876 -0.0001301243
## [1] -0.0003989734 -0.0001398646
## [1] 0.0001158080 0.0001426247
## [1] 0.0001837859 0.0007530687

Exemplos de estimação:

Exemplo 5 - Minimização da soma dos quadrados dos resíduos

Seja uma regressão linear por meio da minimização da soma dos quadrados dos resíduos. Seja, por exemplo. o conjunto de dados cars, em que consideraremos a variável dist (distância de frenagem) como resposta e a variável speed (velocidade) como explicativa. O gráfico dos pontos é apresentado abaixo:

Busca-se a reta \(y = \beta_0 + \beta_1 x\) que passe em meio a essa nuvem de pontos \((x_i, y_i)\) e que atenda a um determinado critério. Deseja-se determinar os valores dos parâmetros \(\beta_0\) e \(\beta_1\) que minimizem uma determinada função objetivo.

Critério: Mínima soma de resíduos ao quadrado

Seja \(\hat{y}_i\) o valor na reta ajustada para \(x_i\), \(i = 1, 2, \dots, n\), ou seja, \(\hat{y}_i = \beta_0 + \beta_1 x_i\). Define-se o resíduo da i-ésima observação como \(r_i = y_i - \hat{y}_i\), em que \(y_i\) é i-ésimo valor observado da resposta. Em nosso exemplo, desenham-se todos os resíduos (\(r_i = y_i - \hat{y}_i\)) para a particular reta \(y_i = -15 + 4 x_i\):

Deseja-se assim encontrar a reta que ofereça a mínima soma de quadrados de resíduos de todos os \(n\) pontos, ou seja, \((\hat{\beta}_0, \hat{\beta}_1) = \arg \min \sum_{i=1}^n r_i^2\). A função objetivo a ser minimizada será então:

\[ f(\beta_o, \beta_1) = \sum_{i=1}^n (y_i - (\beta_0 + \beta_1 x_i))^2 \]

em que, as maiores distâncias serão mais penalizadas (os resíduos são elevados ao quadrado).

 # Função objetivo: soma dos quadrados do resíduos
SQR <- function(dados, par){
    with(dados, sum((par[1] + par[2]*dados[, 1] - dados[, 2])^2))
}

  # valores dos parâmetros que minimizam SQR
par.SQR <- optim(par = c(0, 1), fn = SQR, dados = cars)
par.SQR
## $par
## [1] -17.578151   3.932216
## 
## $value
## [1] 11353.52
## 
## $counts
## function gradient 
##       91       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Esse é o conhecido procedimento de regressão linear por mínimos quadrados ordinários e é bem estabelecido na literatura. Ele pode ser obtido obtido também por meio da função lm(), ou seja:

 # parametros de regressão linear por EMQO - função lm()
lm(dist ~ speed, data = cars)
## 
## Call:
## lm(formula = dist ~ speed, data = cars)
## 
## Coefficients:
## (Intercept)        speed  
##     -17.579        3.932

Percebe-se que são alcançados valores muito próximos daqueles que obtivemos com o procedimento de otimização. O gráfico da reta de regressão de mínimos quadrados, \(\hat{y}= -17.578 + 3.932x\), encontra-se desenhada a seguir:

Critério: Mínima soma de resíduos absolutos

Deseja-se nesse caso a reta que ofereça a mínima soma de resíduos absolutos de todos os \(n\) pontos. A função objetivo a ser minimizada será então:

\[ f(\beta_o, \beta_1) = \sum_{i=1}^n |y_i - (\beta_0 + \beta_1 x_i)| \]

em que, todos os desvios terão a mesma penalização (consideram-se os valores absolutos dos resíduos). Em nosso exemplo com o conjunto de dados cars temos:

  # Função objetivo: soma dos valores absolutos dos resíduos
SRA <- function(dados, par){
    with(dados, sum(abs(par[1] + par[2]*dados[, 1] - dados[, 2])))
}
  # valores dos parâmetros que minimizam SRA
par.SRA <- optim(par = c(0, 1), fn = SRA, dados = cars)
par.SRA
## $par
## [1] -11.600071   3.400006
## 
## $value
## [1] 563.8001
## 
## $counts
## function gradient 
##      145       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

Percebe-se diferença entre os parâmetros, já que os critérios não são equivalentes (Por que?)

##                   SQR        SRA
## Intercepto -17.578151 -11.600071
## Inclinação   3.932216   3.400006

Cada ajuste atende a um dos critério, não sendo o resultado ótimo com respeito ao outro critério. Perceba abaixo os valores de cada função objetivo para os dois ajustes realizados:

##           SQR      SRA
## Aj.1 11353.52 578.9935
## Aj.2 11988.27 563.8001

Perceba a influência dos critérios em cada uma das retas ajustadas:

Critério: Mínima soma de distâncias à reta (regressão linear ortogonal)
Exemplo 6 - Máxima verossimilhança

Deseja-se ajustar uma distribuição de Poisson ao seguinte conjunto de dados de contagem:

 # conjunto de dados de contagem
contagem <- c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 17, 42, 43)
freq <- c(1392, 1711, 914, 468, 306, 192, 96, 56, 35, 17, 15, 6, 2, 2, 1, 1)

  # vetor de dados
x <- rep(contagem, freq)

cuja distribuição empírica dos dados é:

  # distribuição empírica dos dados
tabela <- prop.table(table(x))
tabela
## x
##            1            2            3            4            5            6 
## 0.2669735328 0.3281549674 0.1752972766 0.0897583429 0.0586881473 0.0368239356 
##            7            8            9           10           11           12 
## 0.0184119678 0.0107403145 0.0067126966 0.0032604526 0.0028768700 0.0011507480 
##           13           17           42           43 
## 0.0003835827 0.0003835827 0.0001917913 0.0001917913

Perceba que há alguns valores raros (duas contagens maiores que 42 em 5214 observações. Esse fato pode ser melhor observado por meio do gráfico da distribuição empírica abaixo:

Supondo que os dados seguem uma distribuição de Poisson, a função de verossimilhança, dados os valores observados \(\mathbf{x} = (x_1, x_2, \dots, x_n)\) é:

\[\begin{equation} \operatorname{L}(\lambda; \mathbf{x}) = \prod_{i=1}^n f(x_i) = \prod_{i=1}^n \frac{\operatorname{e}^{-\lambda}\lambda^{x_i}}{x_i!} \end{equation}\]

Dessa maneira, a função de log-verossimilhança para o ajuste de Poisson é:

\[\begin{align} l(\lambda; \mathbf{x}) & = \sum_{i=1}^n \left(\log \operatorname{e}^{-\lambda} + \log \lambda^{x_i} - \log (x_i!)\right)\\ & = \sum_{i=1}^n \left( -\lambda +x_i \log \lambda - \log (x_i!) \right) \end{align}\]

O objetivo do procedimento é encontrar o valor de \(\lambda\) que maximiza a soma das log-verossimilhanças das n observações, a qual é a função objetivo desse procedimento de otimização. Assim sendo, para usar a função optim() necessitamos utilizar \(-l(\lambda, \mathbf{x}\) no argumento fn dessa função de otimização, pois, por deafault, ela busca o mínimo da função (usaremos assim -sum(...):

  # função de log-verossimilhança da Poisson (negativo)
l.pois <- function(x, lambda){ 
              - sum(x * log(lambda) - log(factorial(x)) - lambda) 
          }

Como o espaço paramétrico de busca é unidimensional, usa-se o método "Brent":

  # busca de valor de parâmetro unidimensional - método "Brent"
lambda.hat <- optim(par = 2, fn = l.pois, x = x, method = "Brent", lower = 2, upper = 3)$par
lambda.hat
## [1] 2.703682

O valor de \(\lambda\) que maximiza a verossimilhança é 2.7036824. Podemos comparar com o resultado da função fitdistr(), do pacote MASS que utiliza tamabém a máxima verossimilhança para ajustar o modelo:

   # ajuste do modelo de Poisson aos dados - Máxima verossimilhança
MASS::fitdistr(x, "Poisson")
##      lambda  
##   2.70368239 
##  (0.02277154)

obtendo o mesmo valor. Da mesma maneira, sabemos que o estimador de máxima verossimilhança do modelo de Poisson é a media dos dados, ou seja:

mean(x)
## [1] 2.703682

Referências