library(tidyverse)
library(here)
library(sloop)
set.seed(36528)Some code from last week:
# some fake data
true_pi <- 0.8
true_mu <- c(0, 2)
true_sigma <- c(1, 0.5)
n <- 250
pop <- rbinom(n, size = 1, prob = 1 - true_pi) + 1
x <- rnorm(n, mean = true_mu[pop], sd = true_sigma[pop])
dmix <- function(x, theta){
pi <- theta[1]
mu_1 <- theta[2]
mu_2 <- theta[3]
sigma_1 <- theta[4]
sigma_2 <- theta[5]
pi * dnorm(x, mean = mu_1, sd = sigma_1) +
(1 - pi) * dnorm(x, mean = mu_2, sd = sigma_2)
}
nllhood_mix <- function(theta, x){
d <- dmix(x, theta)
-1 * sum(log(d))
}Recall:
mle1 <- optim(par = c(1, 0, 1.8, 1.2, 0.7), fn = nllhood_mix, x = x, method = "BFGS")
mle1## $par
## [1] 0.9362468 0.3054438 2.2488241 1.2412696 0.1166007
##
## $value
## [1] 411.6919
##
## $counts
## function gradient
## 125 41
##
## $convergence
## [1] 0
##
## $message
## NULLmle2 <- optim(par = c(0.5, 0, 1, 2, 1), fn = nllhood_mix, x = x, method = "BFGS")
mle2## $par
## [1] -3410.0610 -6066.8409 -230.4462 550.6597 588.7451
##
## $value
## [1] -190.2385
##
## $counts
## function gradient
## 185 100
##
## $convergence
## [1] 1
##
## $message
## NULLYour Turn What do you do with this output? How might you improve what the output looks like to help you as a user?
Our goal: To get a more user friendly output from optim()
A programming paradigm where objects are central. Two important terms:
Class defines what an object is
Methods describe what an object can do
Gives us a way to have the same function operate differently on different classes of object, e.g. summary()
z <- factor(sample(letters[1:5], size = 10, replace = TRUE))
y <- rnorm(n = 10)
summary(z)## a b c d e
## 1 2 3 3 1summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.6029 -0.3881 0.3527 0.2827 1.0383 2.3370There is more than one OOP system in R, we’ll focus on S3.
The most common OO system R. Also the simplest.
S3 allows your functions to return rich results with user-friendly display and programmer-friendly internals. S3 is used throughout base R, so it’s important to master if you want to extend base R functions to work with new types of input.
Primarily to make your functions easier to use (for you or others), by creating a new class and defining some methods:
plot(), summary()More advanced make it easy for others to extend your work for other data types, by supplying generic functions.
As our functions get more complicated they need to return complicated things.
E.g. Linear models
mod <- lm(mpg ~ wt, data = mtcars)
str(mod)## List of 12
## $ coefficients : Named num [1:2] 37.29 -5.34
## ..- attr(*, "names")= chr [1:2] "(Intercept)" "wt"
## $ residuals : Named num [1:32] -2.28 -0.92 -2.09 1.3 -0.2 ...
## ..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
## $ effects : Named num [1:32] -113.65 -29.116 -1.661 1.631 0.111 ...
## ..- attr(*, "names")= chr [1:32] "(Intercept)" "wt" "" "" ...
## $ rank : int 2
## $ fitted.values: Named num [1:32] 23.3 21.9 24.9 20.1 18.9 ...
## ..- attr(*, "names")= chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
## $ assign : int [1:2] 0 1
## $ qr :List of 5
## ..$ qr : num [1:32, 1:2] -5.657 0.177 0.177 0.177 0.177 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:32] "Mazda RX4" "Mazda RX4 Wag" "Datsun 710" "Hornet 4 Drive" ...
## .. .. ..$ : chr [1:2] "(Intercept)" "wt"
## .. ..- attr(*, "assign")= int [1:2] 0 1
## ..$ qraux: num [1:2] 1.18 1.05
## ..$ pivot: int [1:2] 1 2
## ..$ tol : num 1e-07
## ..$ rank : int 2
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 30
## $ xlevels : Named list()
## $ call : language lm(formula = mpg ~ wt, data = mtcars)
## $ terms :Classes 'terms', 'formula' language mpg ~ wt
## .. ..- attr(*, "variables")= language list(mpg, wt)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "mpg" "wt"
## .. .. .. ..$ : chr "wt"
## .. ..- attr(*, "term.labels")= chr "wt"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(mpg, wt)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
## .. .. ..- attr(*, "names")= chr [1:2] "mpg" "wt"
## $ model :'data.frame': 32 obs. of 2 variables:
## ..$ mpg: num [1:32] 21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
## ..$ wt : num [1:32] 2.62 2.88 2.32 3.21 3.44 ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language mpg ~ wt
## .. .. ..- attr(*, "variables")= language list(mpg, wt)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "mpg" "wt"
## .. .. .. .. ..$ : chr "wt"
## .. .. ..- attr(*, "term.labels")= chr "wt"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(mpg, wt)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
## .. .. .. ..- attr(*, "names")= chr [1:2] "mpg" "wt"
## - attr(*, "class")= chr "lm"Complicated objects are built on top of simple objects (vectors, matrices, lists etc.) with the additional of attributes.
E.g. lm() returns a list along with the a class attribute.
The class attribute is special, it’s how R knows what kind of object something is.
Use class() to see the class,
class(mod) # the class of this object## [1] "lm"typeof() to see the base type the object is built on,
typeof(mod) # the base type this object is built on## [1] "list"and attributes() to see all attributes:
attributes(mod)## $names
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"
##
## $class
## [1] "lm"Your turn What class is z? What base object is it built on? What other attributes does this object have?
z <- factor(sample(letters[1:5], size = 10, replace = TRUE))
class(z)## [1] "factor"typeof(z)## [1] "integer"attributes(z)## $levels
## [1] "a" "b" "c" "d" "e"
##
## $class
## [1] "factor"Does everything have a class? class() will return something, even if the object doesn’t have a class attribute:
class(y)## [1] "numeric"attributes(y)## NULLIn S3 methods are associated with generic functions (generics for short).
This is different to most other OO languages out there, which use encapsulated OO - the methods are associated with the object.
When an S3 generic function is called, it looks at the class of the object passed to it, and searches for the appropriate method, of the form generic.class, e.g.
summary(mod) # calls summary.lm##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.5432 -2.3647 -0.1252 1.4096 6.8727
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 37.2851 1.8776 19.858 < 2e-16 ***
## wt -5.3445 0.5591 -9.559 1.29e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared: 0.7528, Adjusted R-squared: 0.7446
## F-statistic: 91.38 on 1 and 30 DF, p-value: 1.294e-10class(mod)## [1] "lm"looks for finds and calls summary.lm()
This process is known as method dispatch. You can see the process with s3_dispatch() in the sloop package:
s3_dispatch(summary(mod))## => summary.lm
## * summary.default=> the method that is called, * a method that is defined but not called.
If you want to know if a function is a generic you can use ftype() in sloop
ftype(summary)## [1] "S3" "generic"ftype(summary.lm)## [1] "S3" "method"Your turn You can use . in a name in R without it being a method. Can you figure out if t.test() is the t method for objects of class test? What about t.data.frame()?
# install.packages("sloop")
library(sloop)
ftype(t.test) # S3 generic, not a method for `t()` for class "test"## [1] "S3" "generic"ftype(t.data.frame) # S3 method for `t()` for class "data.frame"## [1] "S3" "method"fake <- data.frame(x = 1:10)
class(fake)## [1] "data.frame"s3_dispatch(t(fake)) # let's you see how R found the method## => t.data.frame
## -> t.defaultUnanswered question: What if we define a frame class and a t.data() generic, then and a frame method for t.data(), do things break?
print() is a genericftype(print)## [1] "S3" "generic"It gets called when you ask for an object
mod # is the same as##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Coefficients:
## (Intercept) wt
## 37.285 -5.344print(mod)##
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
##
## Coefficients:
## (Intercept) wt
## 37.285 -5.344Behind the scenes it finds the print.lm() function
s3_dispatch(print(mod))## => print.lm
## * print.defaultclass doesn’t get any special treatmentYou can set the class of an object just like any usual assignment:
class(mod) <- "table"
mod## coefficients
## 37.285126, -5.344472
## residuals
## -2.2826106, -0.9197704, -2.0859521, 1.2973499, -0.2001440, -0.6932545, -3.9053627, 4.1637381, 2.3499593, 0.2998560, -1.1001440, 0.8668731, -0.0502472, -1.8830236, 1.1733496, 2.1032876, 5.9810744, 6.8727113, 1.7461954, 6.4219792, -2.6110037, -2.9725862, -3.7268663, -3.4623553, 2.4643670, 0.3564263, 0.1520430, 1.2010593, -4.5431513, -2.7809399, -3.2053627, -1.0274952
## effects
## -113.6497374, -29.1157217, -1.6613339, 1.6313943, 0.1111305, -0.3840041, -3.6072442, 4.5003125, 2.6905817, 0.6111305, -0.7888695, 1.1143917, 0.2316793, -1.6061571, 1.3014525, 2.2137818, 6.0995633, 7.3094734, 2.2421594, 6.8956792, -2.2010595, -2.6694078, -3.4150859, -3.1915608, 2.7346556, 0.8200064, 0.5948771, 1.7073457, -4.2045529, -2.4018616, -2.9072442, -0.6494289
## rank
## 2
## fitted.values
## 23.282611, 21.919770, 24.885952, 20.102650, 18.900144, 18.793255, 18.205363, 20.236262, 20.450041, 18.900144, 18.900144, 15.533127, 17.350247, 17.083024, 9.226650, 8.296712, 8.718926, 25.527289, 28.653805, 27.478021, 24.111004, 18.472586, 18.926866, 16.762355, 16.735633, 26.943574, 25.847957, 29.198941, 20.343151, 22.480940, 18.205363, 22.427495
## assign
## 0, 1
## qr
## -5.656854249, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, 0.176776695, -18.199514334, 5.447820482, 0.148230003, -0.016055881, -0.057356801, -0.061027994, -0.081219555, -0.011466889, -0.004124504, -0.057356801, -0.057356801, -0.172999378, -0.110589098, -0.119767081, -0.389599760, -0.421539139, -0.407037927, 0.170257160, 0.277639553, 0.237256431, 0.121613854, -0.072041573, -0.056439003, -0.130780659, -0.131698458, 0.218900467, 0.181270739, 0.296362637, -0.007795696, 0.065628162, -0.081219555, 0.063792566, 1.176776695, 1.046354399, 1.000000000, 2.000000000, 0.000000100, 2.000000000
## df.residual
## 30
## xlevels
## NULL
## call
## lm(formula = mpg ~ wt, data = mtcars)
## terms
## mpg ~ wt
## model
## 21.000, 21.000, 22.800, 21.400, 18.700, 18.100, 14.300, 24.400, 22.800, 19.200, 17.800, 16.400, 17.300, 15.200, 10.400, 10.400, 14.700, 32.400, 30.400, 33.900, 21.500, 15.500, 15.200, 13.300, 19.200, 27.300, 26.000, 30.400, 15.800, 19.700, 15.000, 21.400, 2.620, 2.875, 2.320, 3.215, 3.440, 3.460, 3.570, 3.190, 3.150, 3.440, 3.440, 4.070, 3.730, 3.780, 5.250, 5.424, 5.345, 2.200, 1.615, 1.835, 2.465, 3.520, 3.435, 3.840, 3.845, 1.935, 2.140, 1.513, 3.170, 2.770, 3.570, 2.780But you probably wouldn’t want to…
If we remove the class of mod
class(mod) <- NULL
mod ## $coefficients
## (Intercept) wt
## 37.285126 -5.344472
##
## $residuals
## Mazda RX4 Mazda RX4 Wag Datsun 710
## -2.2826106 -0.9197704 -2.0859521
## Hornet 4 Drive Hornet Sportabout Valiant
## 1.2973499 -0.2001440 -0.6932545
## Duster 360 Merc 240D Merc 230
## -3.9053627 4.1637381 2.3499593
## Merc 280 Merc 280C Merc 450SE
## 0.2998560 -1.1001440 0.8668731
## Merc 450SL Merc 450SLC Cadillac Fleetwood
## -0.0502472 -1.8830236 1.1733496
## Lincoln Continental Chrysler Imperial Fiat 128
## 2.1032876 5.9810744 6.8727113
## Honda Civic Toyota Corolla Toyota Corona
## 1.7461954 6.4219792 -2.6110037
## Dodge Challenger AMC Javelin Camaro Z28
## -2.9725862 -3.7268663 -3.4623553
## Pontiac Firebird Fiat X1-9 Porsche 914-2
## 2.4643670 0.3564263 0.1520430
## Lotus Europa Ford Pantera L Ferrari Dino
## 1.2010593 -4.5431513 -2.7809399
## Maserati Bora Volvo 142E
## -3.2053627 -1.0274952
##
## $effects
## (Intercept) wt
## -113.6497374 -29.1157217 -1.6613339 1.6313943 0.1111305
##
## -0.3840041 -3.6072442 4.5003125 2.6905817 0.6111305
##
## -0.7888695 1.1143917 0.2316793 -1.6061571 1.3014525
##
## 2.2137818 6.0995633 7.3094734 2.2421594 6.8956792
##
## -2.2010595 -2.6694078 -3.4150859 -3.1915608 2.7346556
##
## 0.8200064 0.5948771 1.7073457 -4.2045529 -2.4018616
##
## -2.9072442 -0.6494289
##
## $rank
## [1] 2
##
## $fitted.values
## Mazda RX4 Mazda RX4 Wag Datsun 710
## 23.282611 21.919770 24.885952
## Hornet 4 Drive Hornet Sportabout Valiant
## 20.102650 18.900144 18.793255
## Duster 360 Merc 240D Merc 230
## 18.205363 20.236262 20.450041
## Merc 280 Merc 280C Merc 450SE
## 18.900144 18.900144 15.533127
## Merc 450SL Merc 450SLC Cadillac Fleetwood
## 17.350247 17.083024 9.226650
## Lincoln Continental Chrysler Imperial Fiat 128
## 8.296712 8.718926 25.527289
## Honda Civic Toyota Corolla Toyota Corona
## 28.653805 27.478021 24.111004
## Dodge Challenger AMC Javelin Camaro Z28
## 18.472586 18.926866 16.762355
## Pontiac Firebird Fiat X1-9 Porsche 914-2
## 16.735633 26.943574 25.847957
## Lotus Europa Ford Pantera L Ferrari Dino
## 29.198941 20.343151 22.480940
## Maserati Bora Volvo 142E
## 18.205363 22.427495
##
## $assign
## [1] 0 1
##
## $qr
## $qr
## (Intercept) wt
## Mazda RX4 -5.6568542 -18.199514334
## Mazda RX4 Wag 0.1767767 5.447820482
## Datsun 710 0.1767767 0.148230003
## Hornet 4 Drive 0.1767767 -0.016055881
## Hornet Sportabout 0.1767767 -0.057356801
## Valiant 0.1767767 -0.061027994
## Duster 360 0.1767767 -0.081219555
## Merc 240D 0.1767767 -0.011466889
## Merc 230 0.1767767 -0.004124504
## Merc 280 0.1767767 -0.057356801
## Merc 280C 0.1767767 -0.057356801
## Merc 450SE 0.1767767 -0.172999378
## Merc 450SL 0.1767767 -0.110589098
## Merc 450SLC 0.1767767 -0.119767081
## Cadillac Fleetwood 0.1767767 -0.389599760
## Lincoln Continental 0.1767767 -0.421539139
## Chrysler Imperial 0.1767767 -0.407037927
## Fiat 128 0.1767767 0.170257160
## Honda Civic 0.1767767 0.277639553
## Toyota Corolla 0.1767767 0.237256431
## Toyota Corona 0.1767767 0.121613854
## Dodge Challenger 0.1767767 -0.072041573
## AMC Javelin 0.1767767 -0.056439003
## Camaro Z28 0.1767767 -0.130780659
## Pontiac Firebird 0.1767767 -0.131698458
## Fiat X1-9 0.1767767 0.218900467
## Porsche 914-2 0.1767767 0.181270739
## Lotus Europa 0.1767767 0.296362637
## Ford Pantera L 0.1767767 -0.007795696
## Ferrari Dino 0.1767767 0.065628162
## Maserati Bora 0.1767767 -0.081219555
## Volvo 142E 0.1767767 0.063792566
## attr(,"assign")
## [1] 0 1
##
## $qraux
## [1] 1.176777 1.046354
##
## $pivot
## [1] 1 2
##
## $tol
## [1] 1e-07
##
## $rank
## [1] 2
##
## attr(,"class")
## [1] "qr"
##
## $df.residual
## [1] 30
##
## $xlevels
## named list()
##
## $call
## lm(formula = mpg ~ wt, data = mtcars)
##
## $terms
## mpg ~ wt
## attr(,"variables")
## list(mpg, wt)
## attr(,"factors")
## wt
## mpg 0
## wt 1
## attr(,"term.labels")
## [1] "wt"
## attr(,"order")
## [1] 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(mpg, wt)
## attr(,"dataClasses")
## mpg wt
## "numeric" "numeric"
##
## $model
## mpg wt
## Mazda RX4 21.0 2.620
## Mazda RX4 Wag 21.0 2.875
## Datsun 710 22.8 2.320
## Hornet 4 Drive 21.4 3.215
## Hornet Sportabout 18.7 3.440
## Valiant 18.1 3.460
## Duster 360 14.3 3.570
## Merc 240D 24.4 3.190
## Merc 230 22.8 3.150
## Merc 280 19.2 3.440
## Merc 280C 17.8 3.440
## Merc 450SE 16.4 4.070
## Merc 450SL 17.3 3.730
## Merc 450SLC 15.2 3.780
## Cadillac Fleetwood 10.4 5.250
## Lincoln Continental 10.4 5.424
## Chrysler Imperial 14.7 5.345
## Fiat 128 32.4 2.200
## Honda Civic 30.4 1.615
## Toyota Corolla 33.9 1.835
## Toyota Corona 21.5 2.465
## Dodge Challenger 15.5 3.520
## AMC Javelin 15.2 3.435
## Camaro Z28 13.3 3.840
## Pontiac Firebird 19.2 3.845
## Fiat X1-9 27.3 1.935
## Porsche 914-2 26.0 2.140
## Lotus Europa 30.4 1.513
## Ford Pantera L 15.8 3.170
## Ferrari Dino 19.7 2.770
## Maserati Bora 15.0 3.570
## Volvo 142E 21.4 2.780s3_dispatch(print(mod))## print.list
## => print.defaultIt just prints like any list. This isn’t very useful output…why?
You can either make the object from it’s base type and add a class (or any other attributes):
y <- 1:6
class(y) <- "my_special_kind_of_vector"
attr(y, "max") <- 6
y## [1] 1 2 3 4 5 6
## attr(,"class")
## [1] "my_special_kind_of_vector"
## attr(,"max")
## [1] 6Or use structure() to do it in one go
y <- structure(1:6, max = 6, class = "my_special_kind_of_vector")
y## [1] 1 2 3 4 5 6
## attr(,"max")
## [1] 6
## attr(,"class")
## [1] "my_special_kind_of_vector"Let’s build a new class to hold optim() output, then write a print() method that prints output in more useful fashion.
General procedure:
Your turn What shall we call the class to hold an optim() result?
my_new_class, optim, my_optim, cl_optim, optim_class, optim_out, optimum
Takes the same structure as existing optim() output - list with certain elements. We aren’t going to check for these elements in our constructor.
https://adv-r.hadley.nz/s3.html#s3-constrcutor
First we’ll write a function that creates an object of class optimum.
By convention, call it new_classname():
new_optimum <- function(l){
# This is where you need to enforce any desired structure
stopifnot(is.list(l))
structure(
l,
class = "optimum"
)
}
new_optimum(list(par = c(1, 2, 3), y = "weird"))## $par
## [1] 1 2 3
##
## $y
## [1] "weird"
##
## attr(,"class")
## [1] "optimum"I’m being a little lazy checking the input, but I’ll assume the incoming list has all the components of an optim() output.
new_optimum(mle1)## $par
## [1] 0.9362468 0.3054438 2.2488241 1.2412696 0.1166007
##
## $value
## [1] 411.6919
##
## $counts
## function gradient
## 125 41
##
## $convergence
## [1] 0
##
## $message
## NULL
##
## attr(,"class")
## [1] "optimum"s3_dispatch(print(new_optimum(mle1)))## print.optimum
## => print.defaulthttps://adv-r.hadley.nz/s3.html#helpers
It’s useful to have a function class_name() that conveniently creates object of this class, which we’ll assume just passes it’s arguments onto optim():
optimum <- function(..., quiet = TRUE){
if(quiet){
suppressWarnings({
opt <- optim(...)
})
} else {
opt <- optim(...)
}
new_optimum(opt)
}optimum(par = c(1, 0, 1.8, 1.2, 0.7),
fn = nllhood_mix, x = x, method = "BFGS") ## $par
## [1] 0.9362468 0.3054438 2.2488241 1.2412696 0.1166007
##
## $value
## [1] 411.6919
##
## $counts
## function gradient
## 125 41
##
## $convergence
## [1] 0
##
## $message
## NULL
##
## attr(,"class")
## [1] "optimum"print.optimum()To define a new method, you just write a regular function with a special name generic.class():
print.optimum <- function(x, ...){
cat("Final parameter values: \n")
cat(x$par)
invisible(x)
}
optimum(par = c(1, 0, 1.8, 1.2, 0.7),
fn = nllhood_mix, x = x, method = "BFGS") ## Final parameter values:
## 0.9362468 0.3054438 2.248824 1.24127 0.1166007Your turn Edit the method to just print the final parameter values.
We should check that value of convergence and print a message if optim() didn’t converge:
print.optimum <- function(x, ...){
if(x$convergence != 0){
cat("Failure to converge \n")
} else {
cat("Coverged! \n")
}
cat("Final parameter values \n")
print(x$par)
invisible(x)
}
optimum(par = c(1, 0, 1.8, 1.2, 0.7),
fn = nllhood_mix, x = x, method = "BFGS") ## Coverged!
## Final parameter values
## [1] 0.9362468 0.3054438 2.2488241 1.2412696 0.1166007optimum(par = c(0.5, 0, 1, 2, 1), fn = nllhood_mix, x = x, method = "BFGS")## Failure to converge
## Final parameter values
## [1] -3410.0610 -6066.8409 -230.4462 550.6597 588.7451Even better we might check some more convergence values and add some color! (might need install.packages("crayon")):
print.optimum <- function(x, ...){
if(x$convergence != 0){
cat(crayon::red("Failure to converge \n"))
if (x$convergence == 1){
cat("Iteration limit `maxit` reached. \n")
} else {
cat(x$message)
}
} else {
cat(crayon::green("Coverged! \n"))
}
cat("Final parameter values \n")
print(x$par)
invisible(x)
}
optimum(par = c(1, 0, 1.8, 1.2, 0.7),
fn = nllhood_mix, x = x, method = "BFGS") ## Coverged!
## Final parameter values
## [1] 0.9362468 0.3054438 2.2488241 1.2412696 0.1166007optimum(par = c(0.5, 0, 1, 2, 1), fn = nllhood_mix, x = x, method = "BFGS")## Failure to converge
## Iteration limit `maxit` reached.
## Final parameter values
## [1] -3410.0610 -6066.8409 -230.4462 550.6597 588.7451Your turn Is coef() a generic? Write a coef() method for optimum that returns the parameter estimates.
ftype(coef) # Yes, it's a generic## [1] "S3" "generic"coef # Our method needs to accept the same arguments as the generic## function (object, ...)
## UseMethod("coef")
## <bytecode: 0x7fd1a7440f08>
## <environment: namespace:stats>coef.optimum <- function(object, ...) {
object$par
}
mle_1a <- optimum(par = c(1, 0, 1.8, 1.2, 0.7),
fn = nllhood_mix, x = x, method = "BFGS")
coef(mle_1a)## [1] 0.9362468 0.3054438 2.2488241 1.2412696 0.1166007Does a similar thing to what we have implemented but it uses the S4 OO system.