Załadownie potrzebnych bibliotek, własnych skryptów oraz przygotowanie głównych parametrów zadania.

library(tidyverse)
library(ggpubr)
library(kableExtra)

source('Funkcje dodatkowe.R', encoding = 'UTF-8')

set.seed(1111)
n =100
sigma2 = 0.5 
B = c(80,100,0.005)
x = seq(n)*10/n

Punkt 1 - Wygenerowanie obserwacji

Do wygenerowania obserwacji \(Y_1,...,Y_n\) użyłam swojej własnej funkcji g, zdefiniowanej w skrypcie Funkcje dodatkowe.R.

Y = g(x,B) + rnorm(n,0,sigma2)

Punkt 2 - Wzrór na funkcję \(g(x,\beta)\) oraz jej gradient

Funkcja \(g(x,\beta)\) dana jest wzorem:

\[ g(\boldsymbol x_i, \boldsymbol{\beta}) = \beta_1\boldsymbol x_{i,1}+\beta_2\boldsymbol x_{i,2}e^{-\beta_3\boldsymbol x_{i,3}}+\varepsilon_i, \]

gdzie \(\boldsymbol x_{i,1}\equiv\boldsymbol 1\), \(\boldsymbol x_{i,2}\equiv \boldsymbol 1\), a \(\boldsymbol x_{i,3}\) to obserwacje zmiennej \(X\) z zadania.

Gradient funkcji \(g\) dany jest wzorem:

\[ \boldsymbol{G} = \left[\frac{\partial g(\boldsymbol{x_i}, \boldsymbol{\beta})}{\partial\beta_j}\right] = \left[\boldsymbol 1, e^{-\beta_3\boldsymbol x_{i,3}}, -\beta_2\boldsymbol x_{i,3}e^{-\beta_3\boldsymbol x_{i,3}}\right] \]

Punkt 3 - Wyznaczenie estymatorów

model = nlmgn(Y, x, c(79,101,0.004), delta = 0.0001)

Algorytm Gaussa-Newtona zawarłam we własnej funkcji nlmgn. W wyniku działania algorytmu zostały wyznaczone estymatory, których wartości zostały umieszczone na wykresie z punktu 6.

Punkt 4 - Ilustracja kolejnych kroków algorytmu

kable(model$tab, digits = c(0,2,2,3,3,2), format = 'html',
      col.names = c("s", 
                    "$\\widehat\\beta_1^{(s)}$", 
                    "$\\widehat\\beta_2^{(s)}$", 
                    "$\\widehat\\beta_3^{(s)}$", 
                    "$(\\widehat\\sigma^2)^{(s)}$", 
                    "$L\\Big({\\widehat\\beta^{(s)}}, (\\widehat\\sigma^2)^{(s)}\\Big)$" )) %>% kable_styling()
s \(\widehat\beta_1^{(s)}\) \(\widehat\beta_2^{(s)}\) \(\widehat\beta_3^{(s)}\) \((\widehat\sigma^2)^{(s)}\) \(L\Big({\widehat\beta^{(s)}}, (\widehat\sigma^2)^{(s)}\Big)\)
0 79.00 101.00 0.004 0.486 -104.33
1 1751.50 -1571.22 0.072 273920.009 -766.42
2 169.55 10.72 0.072 0.600 -114.89
3 169.55 10.71 0.058 0.301 -80.29
4 169.28 10.99 0.057 0.300 -80.14
5 169.29 10.98 0.057 0.300 -80.14
6 169.28 10.99 0.057 0.300 -80.14
7 169.28 10.99 0.057 0.300 -80.14
8 169.28 10.99 0.057 0.300 -80.14

Punkt 5 - Porównanie \(\hat\beta\) z \(\beta\)

Porównanie estymowanych wartości współczynników regresji można dokonać obserwując wykres z punktu 6, gdzie zostało zamieszczone równanie zarówno dla \(Y_i\) jak i dla \(\hat{Y_i}\). Można tutaj zauważyć, że co prawda, estymowane wartości tych współczynników dość znacznie odbiegają od wartości rzeczywistych, to jednak estymowana krzywa regresji tylko nieznacznie różni się swoim kształtem od rzeczywistej krzywej.

Punkt 6 - Stworzenie wykresów

hatB = model$B
hatsigma = model$sigma2

lab1 = bquote(y[i]~ "=" ~ .(B[1]) ~ "+" ~ .(B[2]) ~ e^{.(B[3])~x[i]}~ ", "~sigma^2~" = "~.(sigma2))
lab2 = bquote(hat(y[i])~"=" ~ .(round(hatB[1],1)) ~ "+" ~ .(round(hatB[2],1)) ~ e^{.(round(hatB[3],3))~x[i]}~ ", "~hat(sigma)^2~" = "~.(round(hatsigma, 3)))

wykres = ggplot(tibble(x=x, Y=Y), aes(x, Y)) + 
  geom_point(pch=21, fill = "blue", alpha=0.5) +
  stat_function(fun=g, args = list(B=B), color = "red") + 
  stat_function(fun=g, args = list(B=hatB), color = "blue") +
  annotate("text", x=6, y=180, label = lab1, color = "red", size = 5) + 
  annotate("text", x=3, y=176.2, label = lab2, color = "blue", size = 5) 

wykres = 
  ggpar(wykres,
        title = "Wykres ilustrujący oryginalny kształt krzywej oraz jej kształt estymowany",
        ylab = "y",
        caption = "A.Fiołka")
  
wykres

Dodatkowo został tu umieszczony jeszcze jeden wykres przedstawiający pierwsze pięć kroków działania algorytmu.

kolory = c("red", "blue", "green", "orange", "yellow")
wykres = ggplot(tibble(x=x, Y=Y), aes(x, Y)) + 
  geom_point(pch=21, fill = "blue", alpha=0.5) +
  stat_function(fun=g, args = list(B=B), color = "black", lty = 2) +
  stat_function(aes(color="s0"), fun=g, args=list(B=model$Bk[[1]])) +
  stat_function(aes(color="s1"), fun=g, args=list(B=model$Bk[[2]])) +
  stat_function(aes(color="s2"), fun=g, args=list(B=model$Bk[[3]])) +
  stat_function(aes(color="s3"), fun=g, args=list(B=model$Bk[[4]])) +
  stat_function(aes(color="s4"), fun=g, args=list(B=model$Bk[[5]])) +
  annotate("text", x=5, y=180, label = lab1, color = "black", size = 5) + 
  scale_colour_manual("s", values = kolory) + 
  coord_cartesian(ylim = c(min(Y), max(Y)))
wykres = 
  ggpar(wykres,
        title = "Wykres przedstawiający pierwsze pięć kroków algorytmu Gausa-Newtona",
        caption = "A.Fiołka")
wykres

Punkt 7 - Oszacowanie macierzy kowariancji dla \(\hat\beta\)

kable(model$CovB, digits = 3, format='html') %>% kable_styling()
20.904 -20.326 0.142
-20.326 19.780 -0.138
0.142 -0.138 0.001

Punkt 8 - Ponowne uruchomienie algorytmu

model2 = nlmgn(Y, x, c(82,99,0.006), delta = 0.0001)
model3 = nlmgn(Y, x, c(10,10,1), delta = 0.0001)

hatB2 = model2$B
hat2sigma = model2$sigma2
hatB3 = model3$B
hat3sigma = model3$sigma2


lab1 = bquote(y[i]~ "=" ~ .(B[1]) ~ "+" ~ .(B[2]) ~ e^{.(B[3])~x[i]}~ ", "~sigma^2~" = "~.(sigma2))
lab2 = bquote(hat(y[i])~"=" ~ .(round(hatB2[1],1)) ~ "+" ~ .(round(hatB2[2],1)) ~ e^{.(round(hatB2[3],3))~x[i]}~ ", "~hat(sigma)^2~" = "~.(round(hat2sigma, 3)))
lab3 = bquote(hat(y[i])~"=" ~ .(round(hatB3[1],1)) ~ "+" ~ .(round(hatB3[2],1)) ~ e^{.(round(hatB3[3],3))~x[i]}~ ", "~hat(sigma)^2~" = "~.(round(hat3sigma, 3)))

wykres = ggplot(tibble(x=x, Y=Y), aes(x, Y)) + 
  geom_point(pch=21, fill = "blue", alpha=0.5) +
  stat_function(fun=g, args = list(B=B), color = "red") + 
  stat_function(fun=g, args = list(B=hatB2), color = "blue") +
  stat_function(fun=g, args = list(B=hatB3), color = "purple", lty = 2) +
  annotate("text", x=6, y=180, label = lab1, color = "red", size = 5) + 
  annotate("text", x=3, y=176.2, label = lab2, color = "blue", size = 5) +
  annotate("text", x=3, y=175.5, label = lab3, color = "purple", size = 5) 

wykres = 
  ggpar(wykres,
        title = paste("Wykres ilustrujący oryginalny kształt krzywej oraz kształt\n",
                      "estymowanych krzywych dla różnych parametrów startowych"),
        ylab = "y",
        caption = "A.Fiołka")
  
wykres

Analizując powyższy wykres oraz zawarte na nim równania estymowanych funkcji, które algorytm wyznaczył dla różnych parametrów startowych, należy zauważyć, że wartości estymowanych współczynników regresji, jak również wartość estymowanej wariancji są dokładnie takie same. Jeżeli tylko algorytm pozostaje zbieżny, zawsze zbiega od do takiego samego wyniku, a otrzymane wartości zależą wyłącznie od wartości poszczególnych obserwacji.

Funkcje dodatkowe

Poniżej został umieszczony kod funkcji dodatkowych.

#Pomocnicze gunkcje g oraz jej pochodne cząstkowe 
g = function(x,B) B[1]+B[2]*exp(-B[3]*x)
g = Vectorize(g,"x")
dgB1 = function(x,B) 1
dgB1 = Vectorize(dgB1,"x")
dgB2 = function(x,B) exp(-B[3]*x)
dgB2 = Vectorize(dgB2,"x")
dgB3 = function(x,B) -B[2]*x*exp(-B[3]*x)
dgB3 = Vectorize(dgB3,"x")
hatsigma2 = function(e,B) sum(e^2)/(length(e)-length(B))
Lfun = function(e,sigma2){
  -(length(e)/2)*log(2*pi) - (length(e)/2)*log(sigma2) - sum(e^2/(2*sigma2))
}

#Funkcja stosujaca algorytm Gausa_Newtona do estymacji parametrów regresji nieliniowej
nlmgn = function(y, x, B, delta = 0.0001, M = 1000){
  m = 0
  Bk = list()
  ek = list()
  sigma2k = list()
  Lk = list()
  tab = tibble()
  
  
  Bn = -B
  repeat{
    
    e = y - g(x,B)
    sigma2 = hatsigma2(e, B)
    L = Lfun(e, sigma2)
    
    Bk = append(Bk, list(B))
    ek = append(ek, list(e))
    sigma2k = append(sigma2k, sigma2)
    Lk = append(Lk, list(L))
    tab = tab %>% bind_rows(tibble(
      krok = m, B1 = B[1], B2 = B[2], B3 = B[3], sigma2 = sigma2, L = L))
    
    if (sqrt(sum((B-Bn)^2)) <= delta) break

    G = cbind(dgB1(x,B), dgB2(x,B), dgB3(x,B))
    detG <<- det(t(G) %*% G)
    if (abs(detG) <= 2*.Machine$double.eps){
      message(paste("Gradient matrix jest macierzą osobliwą!\n",
                    "Wyznacznik det(t(G) * G) wynosi ", detG))
      break
    }
    if (is.nan(detG)){
      message("Gradient matrix zawiera wartości Inf!")
      break
    }
    if (m == (M-1)){
      message(paste("Algorytm po wykonaniu", m, "kroków, nie zbiega się!"))
      break
    }
    
    Bn = B
    B = Bn + solve(t(G) %*% G) %*% t(G) %*% e
    CovB = sigma2*solve(t(G) %*% G)
    m = m + 1
  }
  
  list(k = m, B = B, e = e, sigma2 = sigma2, L = L,
    Bk = Bk, ek = ek, sigma2k = sigma2k, Lk = Lk, 
    tab = tab, CovB = CovB
  )
}