Computação Estatística II
Atividade de avaliação 2: Métodos de Otimização em Estatística
03/05/2021
Questão 1
As técnicas de otimização são utilizadas para encontrar determinados valores que maximizam ou minimizam uma função de interesse. Tais métodos tornaram-se muito relevantes na estimação estatística, ajustes de modelo, etc. Um problema básico que pode ser considerado é encontrar estimativas dos dois parâmetros num modelo linear simples relacionado a variável dependente, \(y\), e uma variável independente, \(x\). O modelo pode ser formulado como: \[y_i= \alpha + \beta x_{i} + \epsilon_{i}\] onde \(x_{i}\), \(y_{i}\) com \(i= 1,...,n\), são os valores da variável independente e dependente consideradas no modelo, e são “erros” ou termos residuais com zero valores esperados, contando por quanto uma observação,\(y_{i}\) defere do seu valor predito \(\alpha + \beta x_{i}\) O problema de encontrar estimativas para \(\alpha\) e \(\beta\) pode ser abordado de diversas maneiras. Talvez o mais comum seja encontrar algum critério de adequação que meça o quão próximo o modelo se adequa com os dados observados e, em seguida, escolher valores para os dois parâmetros que minimizam a medida de ajuste escolhida. Um critério óbvio de adequação para o modelo de regressão linear simples é a soma dos quadrados dos termos de erro em (1.1), ou seja, \[S= \sum_{i=1}^n \epsilon_{i}^2\] \(S\) mede o quão bem os valores observados da variável dependente se ajustam aos previstos pelo modelo, com menores valores de \(S\) indicando um melhor ajuste. Consequentemente, escolher como estimativas para \(\alpha\) e \(\beta\) valores que minimizam \(S\) é um procedimento intuitivamente razoável e, é claro, nada menos do que a bem conhecida técnica de estimativa de mínimos quadrados. Outro problema de estimativa comum em estatística surge quando desejamos estimar o parâmetro (ou parâmetros) de uma função de densidade de probabilidade dada uma amostra aleatória retirada da função de densidade. Por exemplo, podemos ter uma amostra de \(n\) valores, \(x_{1},...,x_{n}\), de uma função de densidade exponencial da forma \[f(x)= \lambda e^{-\lambda x}\] e desejamos estimar. Um método útil nesta situação é formar a função de densidade de probabilidade conjunta das observações, ou seja, \[\ell(x_{1},...,x_{1}; \lambda) = \prod_{i=1}^n \lambda e^{-\lambda x_{i}}\] e escolher como estimativa de \(\lambda\) o valor que maximiza \(\ell\), ou seja, uma estimativa de máxima verossimilhança.
Os dois problemas de estimativa descritos acima podem ser formulados em termos de otimização de alguma função numérica relacionados a uma série de parâmetros, e existem ainda muitos outros problemas estatísticos podem ser formulados de maneira semelhante. O presente trabalho visa elaborar um breve estudo a partir da resolução de alguns problemas de otimização em estatística a partir do software R. Além dos métodos descritos acima, é possível ainda modelarmos em R outros problemas como, por exemplo, minimizar a soma residual dos quadrados no âmbito da obtenção das estimativas dos mínimos quadrados no modelo de regressão linear pelo método dos mínimos quadrados (Shalabh, 2016). Além disso, a otimização é fundamental para estimar modelos de regressão não linear (MRNL) através da otimização direta que consiste em derivar a soma dos quadrados dos erros em relação a cada coeficiente ou parâmetro desconhecido, igualamos a zero a equação resultante e resolvemos simultaneamente as equações normais resultantes (Gujarati & José, 2004).
plot(PIB,IPB,pch=16,col=4)
abline(a=0.25352, b=-926.09039)
Regra <- lm(IPB~PIB)
Regra
Call:
lm(formula = IPB~PIB)
Coefficients:
(Intercept) PIB
-926.0904 0.2535
sum(residuals(Regra)^2)/(Regra$df.residual)
4412.537
SQE <- sum(residuals(Regra)^2)
SQT <. sum(IPB-mean(IPB)^2)
1- (SQE/SQT)
0.9647772Questão 2
a) No Método da Secante o cálculo de \(x_{n+1}\) é feito à custa de \(x_{n}\) e \(x_{n-1}\). Necessita-se de dois valores iniciais para desencadear o processo iterativo da secante.
Temos: \[x_{n+1} = x_{n}- \frac{f(x_{n})}{\frac{f(x_{n})-f(x_{n-1})}{x_{n}-x{n-1}}}\]
Em seguida temos o programa para o método em questão.
R: Método da secante
sec = function(f, a, b, tot = 10^-7, maxiter = 100){
iter = 0 #x0=a e x1=b
while(abs(b-a)>tot){
iter = iter+1
if(iter > maxiter) return("Atingido o limite de iterações.")
D = (f(b)-f(a))/(b-a)
a = b
b = b-f(b)/D
}
return(list(x = b, iter = iter))
}Por forma a obter a raiz da função \(cosx-x\) usa-se sec.
ftn <- function(x){cos(x)-x}
sec(ftn, 1, 2, tot = 10^-7, maxiter = 100)## $x
## [1] 0.7390851
##
## $iter
## [1] 5b) Trata-se seguidamente de comparar o desempenho de três métodos distintos na determinação da raíz da função \(log(x)-e^{-x}\) com \(x_{0}=1\) e \(x_{1}=2\).
Inicia-se com contextualização do método de Newton-raphson e a condição suficiente de convergência.
Sendo \(f(x)= 0\) uma equação com uma raiz real separada no intervalo \(X=[a, b]\), onde \(f(a)\) e \(f(b)\) têm sinais contrário,
- \(f\), \(f'\), \(f''\) funções contínuas em \(X\),
- \(f'(x)\) não se anula em \(X\),
- \(f''(x)\) não se anula em \(X\),
então a fómula iteradora \(x_{n+1}= x_{n}- \frac{f(x_{n})}{f'(x_{n})}\), gera uma sucessão convergente para a raiz. Salienta-se que a condição enunciada é apenas suficiente, sendo que a condição pode não se cumprir e o método mesmo assim ser convergente. 2
R: Método de Newton-Raphson
metodo.nr <- function(f, a, b, tol = 1e-5, n = 1000) {
require(numDeriv) # pacote para derivada de f'(x)
x0 <- a # Valor inicial
k <- n # Contagem de iteracoes
# Verificar os limtes superior e inferior se = 0
fa <- f(a)
if (fa == 0.0) {
return(a)
}
fb <- f(b)
if (fb == 0.0) {
return(b)
}
for (i in 1:n) {
dx <- genD(func = f, x = x0)$D[1] #f'(x0)
x1 <- x0 - (f(x0) / dx) # Calculo do proximo valor x1
k[i] <- x1 # Armazena x1 no vetor k
# Uma vez que a diferenca x0-x1 satisfatoria, print output
if (abs(x1 - x0) < tol) {
root.approx <- tail(k, n=1)
res <- list('Numero de Iter' = i, 'Iter' = k, 'Raiz' = root.approx)
return(res)
}
# Se a Newton-Raphson nao atingir convergencia, fazer x1 = x0 e continuar
x0 <- x1
}
print('Iteracoes demasiadas no metodo!') ### em caso de repeticoes 'sem fim'!
}Entenda-se que o método do ponto fixo em termos de qualidades intrínsecas do algoritmo são fracas, onde a convergência não está garantida.Porém, é dos métodos mais usados no ramo da matemátoca pela sua riqueza no conceito base.
Se a função \(g:D \subseteq \mathbb{R} \rightarrow \mathbb{R}\) for contínua no intervalo \(J = [a,b] \subseteq D\), com \(g(J) \subset J\), então \(g\) tem um ponto fixo em \(J\)
R: Método do ponto fixo
pontofixo <- function(ftn, x0, tol = 1e-9, max.iter = 100) {
xold <- x0
xnew <- ftn(xold)
iter <- 1
cat("Na iteracao 1 o valor de x é:", xnew, "\n")
# continuar iteracoes até o valor de x atingir tol
while ((abs(xnew-xold) > tol) && (iter < max.iter)) {
xold <- xnew;
xnew <- ftn(xold);
iter <- iter + 1
cat("Na iteraÇão", iter, "o valor de x é:", xnew, "\n")
}
# output depende do sucesso do algoritmo
if (abs(xnew-xold) > tol) {
cat("O algoritimo falhou em convergir\n")
return(NULL)
} else {
cat("Algoritimo convergiu!\n")
return(xnew)
}
}R: Método da secante O programa para este método encontra-se delinhado na alínea a.
De seguida estrutura-se a função em R:
func.log <- function(x) {
log(x)-exp(-x)
}
x0 <-1
x1 <-2Estamos em condições de efetuar as comparações de desempenho entre os diferentes métodos.
sec(func.log, x0, x1)## $x
## [1] 1.3098
##
## $iter
## [1] 6metodo.nr(func.log, x0, x1)## $`Numero de Iter`
## [1] 4
##
## $Iter
## [1] 1.268941 1.309108 1.309799 1.309800
##
## $Raiz
## [1] 1.3098funcaog <- function(x){exp(exp(-x))}
pontofixo(funcaog,x0)## Na iteracao 1 o valor de x é: 1.444668
## Na iteraÇão 2 o valor de x é: 1.265952
## Na iteraÇão 3 o valor de x é: 1.32574
## Na iteraÇão 4 o valor de x é: 1.304222
## Na iteraÇão 5 o valor de x é: 1.311778
## Na iteraÇão 6 o valor de x é: 1.309101
## Na iteraÇão 7 o valor de x é: 1.310047
## Na iteraÇão 8 o valor de x é: 1.309712
## Na iteraÇão 9 o valor de x é: 1.30983
## Na iteraÇão 10 o valor de x é: 1.309789
## Na iteraÇão 11 o valor de x é: 1.309803
## Na iteraÇão 12 o valor de x é: 1.309798
## Na iteraÇão 13 o valor de x é: 1.3098
## Na iteraÇão 14 o valor de x é: 1.309799
## Na iteraÇão 15 o valor de x é: 1.3098
## Na iteraÇão 16 o valor de x é: 1.3098
## Na iteraÇão 17 o valor de x é: 1.3098
## Na iteraÇão 18 o valor de x é: 1.3098
## Na iteraÇão 19 o valor de x é: 1.3098
## Na iteraÇão 20 o valor de x é: 1.3098
## Na iteraÇão 21 o valor de x é: 1.3098
## Algoritimo convergiu!## [1] 1.3098Os métodos possuem um número baixo de iterações para obtenção da raiz. Porém, o método Newton-Raphson foi mais eficaz, realizando apenas 4 iterações perante 6 iterações do método da secante. Já o metodo do ponto fixo teve duas desvantagens: o calculo de g(x) e necessitou de pelo menos 9 iterações para encontrar o valor com 4 casas decimais.
Questão 3
Na resolução desta questão faz-se o uso da biblioteca ROI. Trata-se de um framework sofisticado para resolução de problemas de otimização em R.
library(ROI)
library(MASS)
######DEFININDO FUNCOES
#FUNCAO OBJETIVO
f <- function(x){x[1]^2+x[2]^2+x[3]^2}
#FUNCOES DE RESTRICAO
g <- function(x){x[1]^2-x[2]^2-x[3]^2}
h <- function(x){2*x[1]+x[3]-2}
### Ajustando parametros com comando OP
nlp <- OP(F_objective(f, n=3), #funcao objetivo + numero de variaveis
F_constraint(F = list(g,h), dir = c("==","=="),
J = NULL, rhs = c(0,0)), maximum = FALSE)
##funcoes restricao para minimo
### Verificar "Solver" (resolvedor) adequado
ROI_applicable_solvers(nlp)## [1] "nlminb"# Ponto inicial para tentativas numéricas
start=c(0,0,0)
# Solucao para (x,y,z)
sol = ROI_solve(nlp, solver = "nlminb",start=start)
s <- solution(sol)
solucao <- data.frame(
var = c("x", "y", "z"),
valor = round(s,3),
fracao = c("2/3","0","2/3"))
solucao## var valor fracao
## 1 x 0.667 2/3
## 2 y 0.000 0
## 3 z 0.667 2/3# Distancia minima calculada com valores substituidos em f(x,y,z)
fractions(f(s))## [1] 8/9Questão 4
###funcao exponencial sem constante
fx1 <- function(x) {
fx <- ifelse(x<0,0,ifelse(x>1,0,exp(-x^2/2)))
return(fx)
}
##encontrando o valor de K -> integral da FDP precisa ser = 1
integral <- (integrate(fx1, lower = 0, upper = +Inf))
k <- 1/integral[["value"]][1]
##adicionando k na expressao
fx2 <- function(x) {
fx <- ifelse(x<0,0,ifelse(x>1,0,k*exp(-x^2/2)))
return(fx)
}
## Verificando valor de k -> integral da FDP com k precisa ser 1
integrate(fx2, lower = 0, upper = +Inf)## 1 with absolute error < 3e-11## Se desejar, pode-se plotar o valor teorico da FDP
plot(fx2, from = -0.1, to = 1.1)
###Aplicando o metodo da reijeicao - funcao
rejeicaoK <- function(fx, a, b, K) {
while (TRUE) {
x <- runif(1, a, b) #valores de x
y <- runif(1, 0, K) #funcao envelope
if (y < fx(x)) return(x)
}
}
set.seed(22)
nreps <- 200000 #repeticoes para metodo da rejeicao
Observacoes <- rep(0, nreps)
for(i in 1:nreps) {
Observacoes[i] <- rejeicaoK(fx2, 0, 1, 1.2) #definindo teto para envelope = 1.2
}
#Se necessario, podemos usar o grafico da funcao de densidade dos valores simulados
hist(Observacoes, breaks = seq(-0.1, 1.1, by=0.01), freq = FALSE,ylim = c(0, 1.2), main = "")
### Utilizando Esperanca E(X)=Int(x*f(x))
fx3 <- function(x) {
fx <- ifelse(x<0,0,ifelse(x>1,0,x*(k*exp(-x^2/2))))
return(fx)
}
#RESPOSTAS
### Utilizando metodo da rejeicao
mean(Observacoes)## [1] 0.459472### Utilizando Esperanca E(X)=Int(x*f(x))
integrate(fx3 , lower = -Inf , upper = +Inf)## 0.4598622 with absolute error < 4.1e-09Questão 5
Existem vários métodos numéricos por forma a obter a função distribuição normal.Opta-se pelo método do trapézio e o método de simpson.
R: Método do trapézio
### Criar funcao da questao 5 com variaveis
funcao5 <- function(x,med=0,desvpad=1){
fx=(1/(sqrt(2*pi*desvpad^2)))*exp(-((x-med)^2/2*desvpad^2))
return(fx)
}
### Funcao Algoritmo trapezoidal para integração numérica
trapez <- function(func,med=0,desvpad=1,inf=-5,sup=5,size=100000) {
vetor=c(seq(inf,sup,size))
result=0
for (i in 2:length(vetor)){
result=
result+
abs(((func(vetor[i],med,desvpad)+
func(vetor[i-1],med,desvpad))/2)*(vetor[i]-vetor[i-1]))
}
return(result)
}
### Exibindo valores - confirmando resultados da integracao
trapez(funcao5,0,1,-100,0,size=1e-4) ##Sobre a média## [1] 0.5trapez(funcao5,0,1,-1,1,size=1e-4) ## intervalo de -1 a 1 desvio padrao## [1] 0.6826895trapez(funcao5,0,1,-3,3,size=1e-4) ## intervalo de -3 a 3 desvios padrao## [1] 0.9973002Sabe-se que por norma as aproximações pela regra do trapézio são, em geral, bastante grosseiras com ordens de interpolação polinomial elavados, muito devido ao fenómeno de Runge.
R: Método de simpson
simpson_n <- function(ftn, a, b, n=100){
n <-max(c(2*(n %% 2), 4))
h <-(b-a)/n
x.vec1 <- seq(a+h, b-h, by=2*h)
x.vec2 <- seq(a+h*2, b-2*h, by=2*h)
f.vec1 <- sapply(x.vec1, ftn)
f.vec2 <- sapply(x.vec2, ftn)
S <- h/3*(ftn(a) + ftn(b) + 4*sum(f.vec1) + 2*sum(f.vec2))
return(S)
}
phi <- function(x) return(exp(-x^2/2)/sqrt(2*pi))
Phi <- function(z){
if(z<0){
return(0.5 - simpson_n(phi, z, 0))
} else{
return(0.5 + simpson_n(phi, 0, z))
}
}
z <- seq(-5, 5, by =0.1)
phi.z <-sapply(z, phi)
Phi.z <-sapply(z, Phi)
plot(z, Phi.z, type = "l", col="blue")
lines(z, phi.z, col="red")
legend(-4, 0.95, legend=c("FDP", "FDA"),
col=c("red", "blue"),lty=1:1,cex=0.8)
Phi(0)## [1] 0.5Phi(1)-Phi(-1)## [1] 0.682711Phi(3)-Phi(-3)## [1] 0.9969588Questão 6
### Definicao da funcao
normalex5 <- function(x,mu,sigma){
fx=(1/(sqrt(2*pi*sigma^2)))*exp(-((x-mu)^2)/(2*sigma^2))
return(fx)
}
### Uso da regra 1/3 para integracao numerica de Simpson
Int.Simpson <-function(zp){
f_a=normalex5(-10,0,1)
f_b=normalex5(zp,0,1)
fab2=normalex5((-10+zp)/(2),0,1)
integral=(f_a+4*fab2+f_b)*((zp-(-10))/6)
return(integral)
}
### Metodo secante com adicao de raiz
metodo.secante <- function(func, xa, xb ,raiz, tol = 1e-5) {
i=0
x1n=0
fx=func(xa)-raiz
fx1=func(xb)-raiz
while (abs(fx) > tol){
x1n=xa-(fx * (xa-xb)/(fx-fx1))
xb=xa
xa=x1n
fx=func(xa)-raiz
fx1=func(xb)-raiz
i=i+1
}
resultado <- c(i,xa)
return(resultado)
}
### Visualizacao da curva da integral no intervalo estabelecido para cada p
d <- seq(-10,10,0.001)
test <- sapply(d, Int.Simpson)
par(mfrow = c(2, 2), mar = c(2,2,2,2))
plot(d,test-0.9, type = "l", ylab = "", main = "p=0.9")
abline(0,0, col = "darkgreen", lty = 3)
abline(v=seq(5,7),lty=2,col="blue")
plot(d,test-0.95, type = "l", ylab = "", main = "p=0.95")
abline(0,0, col = "darkgreen", lty = 3)
abline(v=seq(5,7),lty=2,col="blue")
plot(d,test-0.975, type = "l", ylab = "", main = "p=0.975")
abline(0,0, col = "darkgreen", lty = 3)
abline(v=seq(5,7),lty=2,col="blue")
plot(d,test-0.99, type = "l", ylab = "", main = "p=0.99")
abline(0,0, col = "darkgreen", lty = 3)
abline(v=seq(5,7),lty=2,col="blue")
#Conforme graficos, valor da raiz ficara entre 6 e 7 para todos p's, logo:
r09 <- metodo.secante(Int.Simpson,6,7,0.9)
r095 <- metodo.secante(Int.Simpson,6,7,0.95)
r0975 <- metodo.secante(Int.Simpson,6,7,0.975)
r099 <- metodo.secante(Int.Simpson,6,7,0.99)
dataf <- data.frame(
p = c("0.9","0.95","0.975","0.99"),
Iteracoes = c(r09[1],r095[1],r0975[1],r099[1]),
Valor = c(r09[2],r095[2],r0975[2],r099[2]))
dataf## p Iteracoes Valor
## 1 0.9 4 6.443708
## 2 0.95 4 6.501067
## 3 0.975 4 6.528943
## 4 0.99 4 6.545428Conforme esperado, o valor aumenta com o aumento de \(p\).
Questão 7
Usa-se o método de newton-raphson para a função bivariada, com solução inicial \(\theta = (0.9, 0.9)\)
R: Método Newton-Raphson
# Instalacao do biblioteca "optimx", que possui ferramentas para otimizacao em R.
#install.packages("optimx")
require(optimx)
#Funcao de Rosenbrock (Questao 7)
fq7 <- function(x){
return(100*(x[2] - x[1]*x[1])^2 + (1-x[1])^2)
}
#1a Derivada -> gradiente
fq7d <- function(x){
return(c(-400*x[1]*(x[2] - x[1]*x[1]) - 2*(1 - x[1]), 200*(x[2] - x[1]*x[1])))
}
#Solucao inicial apresentada no enunciado
x0 <- c(0.9, 0.9)
# Faz-se uso da funcao optim do pacote optimx para otimizacao em R.
# O metodo escolhido para resolucao o metodo Quase-Newton BFGS (que vem de seus
# autores Broyden, Fletcher, Goldfarb e Shanno). O metodo utiliza valores da
# funcao e gradientes para contruir uma 'foto' da superficie a ser otimizada.
optim(x0,fq7,fq7d,method="BFGS")## $par
## [1] 1 1
##
## $value
## [1] 7.103903e-20
##
## $counts
## function gradient
## 35 15
##
## $convergence
## [1] 0
##
## $message
## NULLA solução da função optim, ‘par’ apresenta os pontos mais próximos encontrados pelo método utilizado. Ainda, ‘counts’ apresenta o número de chamadas para fq7 e seu gradiente, respectivamente. O valor ‘convergence’ com o resultado 0 informa que o método convergiu, ou seja, foi efetivo em chegar a um valor próximo ao real. 3 4
Bibliografia
Everitt, B.S. (1987). Introduction to Optimization Methods and their Application in Statistics, Dordrecht Springer, Netherlands.
Gujarati, D. N., & José,M.(2004). Econometria Básica. Elsevier,Campus.
Serranho, P., (2017) Matemática Aplicada e Análise Numérica, Secção de Matemática, Departamento de Ciências e Tecnologia, Universidade Aberta.
Shalabh, H. C. (2016). Introduction to statistics and data analysis-with exercises, solutions an. Springer International Publish.
Ver exemplo completo em: Gujarati & José, 2004, pp. 172↩
Ver Serranho,P., Matemática Aplicada e Análise Numérica, Secção de Matemática, Departamento de Ciências e Tecnologia, Universidade Aberta, 2017↩
Código fonte optim encontra-se em: https://rdrr.io/rforge/optreplace/src/R/optim.R↩
Ver documentação disponível em: https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/optim↩