AST 430 all
Monowar Hossain
2024-12-02
Lecture
Lecture 1
meanfun<- function(x){
m<- sum(x)/length(x)
return(m)
}
x<-c(4,5,6,7)
meanfun(x)## [1] 5.5 library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors#install.packages("palmerpenguins")
library(palmerpenguins)
glimpse(penguins)## Rows: 344
## Columns: 8
## $ species <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adel…
## $ island <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgerse…
## $ bill_length_mm <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, …
## $ bill_depth_mm <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, …
## $ flipper_length_mm <int> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186…
## $ body_mass_g <int> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, …
## $ sex <fct> male, female, female, NA, female, male, female, male…
## $ year <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007…library(gapminder)
glimpse(gapminder)## Rows: 1,704
## Columns: 6
## $ country <fct> "Afghanistan", "Afghanistan", "Afghanistan", "Afghanistan", …
## $ continent <fct> Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asia, Asia, …
## $ year <int> 1952, 1957, 1962, 1967, 1972, 1977, 1982, 1987, 1992, 1997, …
## $ lifeExp <dbl> 28.801, 30.332, 31.997, 34.020, 36.088, 38.438, 39.854, 40.8…
## $ pop <int> 8425333, 9240934, 10267083, 11537966, 13079460, 14880372, 12…
## $ gdpPercap <dbl> 779.4453, 820.8530, 853.1007, 836.1971, 739.9811, 786.1134, …#mtcars
mtca<-as_tibble(mtcars)
mtca## # A tibble: 32 × 11
## mpg cyl disp hp drat wt qsec vs am gear carb
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 21 6 160 110 3.9 2.62 16.5 0 1 4 4
## 2 21 6 160 110 3.9 2.88 17.0 0 1 4 4
## 3 22.8 4 108 93 3.85 2.32 18.6 1 1 4 1
## 4 21.4 6 258 110 3.08 3.22 19.4 1 0 3 1
## 5 18.7 8 360 175 3.15 3.44 17.0 0 0 3 2
## 6 18.1 6 225 105 2.76 3.46 20.2 1 0 3 1
## 7 14.3 8 360 245 3.21 3.57 15.8 0 0 3 4
## 8 24.4 4 147. 62 3.69 3.19 20 1 0 4 2
## 9 22.8 4 141. 95 3.92 3.15 22.9 1 0 4 2
## 10 19.2 6 168. 123 3.92 3.44 18.3 1 0 4 4
## # ℹ 22 more rows#filter()
#select()
#mutate()
#arrange()
#summarise()
#group_by()
filter(.data = penguins,species == "Chinstrap")## # A tibble: 68 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Chinstrap Dream 46.5 17.9 192 3500
## 2 Chinstrap Dream 50 19.5 196 3900
## 3 Chinstrap Dream 51.3 19.2 193 3650
## 4 Chinstrap Dream 45.4 18.7 188 3525
## 5 Chinstrap Dream 52.7 19.8 197 3725
## 6 Chinstrap Dream 45.2 17.8 198 3950
## 7 Chinstrap Dream 46.1 18.2 178 3250
## 8 Chinstrap Dream 51.3 18.2 197 3750
## 9 Chinstrap Dream 46 18.9 195 4150
## 10 Chinstrap Dream 51.3 19.9 198 3700
## # ℹ 58 more rows
## # ℹ 2 more variables: sex <fct>, year <int>filter(.data = penguins,species == "Chinstrap" |species == "Gentoo")## # A tibble: 192 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Gentoo Biscoe 46.1 13.2 211 4500
## 2 Gentoo Biscoe 50 16.3 230 5700
## 3 Gentoo Biscoe 48.7 14.1 210 4450
## 4 Gentoo Biscoe 50 15.2 218 5700
## 5 Gentoo Biscoe 47.6 14.5 215 5400
## 6 Gentoo Biscoe 46.5 13.5 210 4550
## 7 Gentoo Biscoe 45.4 14.6 211 4800
## 8 Gentoo Biscoe 46.7 15.3 219 5200
## 9 Gentoo Biscoe 43.3 13.4 209 4400
## 10 Gentoo Biscoe 46.8 15.4 215 5150
## # ℹ 182 more rows
## # ℹ 2 more variables: sex <fct>, year <int>filter(.data = penguins,sex == "female" & island == "Dream")## # A tibble: 61 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Dream 39.5 16.7 178 3250
## 2 Adelie Dream 39.5 17.8 188 3300
## 3 Adelie Dream 36.4 17 195 3325
## 4 Adelie Dream 42.2 18.5 180 3550
## 5 Adelie Dream 37.6 19.3 181 3300
## 6 Adelie Dream 36.5 18 182 3150
## 7 Adelie Dream 36 18.5 186 3100
## 8 Adelie Dream 37 16.9 185 3000
## 9 Adelie Dream 36 17.9 190 3450
## 10 Adelie Dream 37.3 17.8 191 3350
## # ℹ 51 more rows
## # ℹ 2 more variables: sex <fct>, year <int>select(.data = penguins,species,starts_with("bill"))## # A tibble: 344 × 3
## species bill_length_mm bill_depth_mm
## <fct> <dbl> <dbl>
## 1 Adelie 39.1 18.7
## 2 Adelie 39.5 17.4
## 3 Adelie 40.3 18
## 4 Adelie NA NA
## 5 Adelie 36.7 19.3
## 6 Adelie 39.3 20.6
## 7 Adelie 38.9 17.8
## 8 Adelie 39.2 19.6
## 9 Adelie 34.1 18.1
## 10 Adelie 42 20.2
## # ℹ 334 more rowspenguins %>%
filter(sex == "female") %>%
filter(island == "Dream") ## # A tibble: 61 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Dream 39.5 16.7 178 3250
## 2 Adelie Dream 39.5 17.8 188 3300
## 3 Adelie Dream 36.4 17 195 3325
## 4 Adelie Dream 42.2 18.5 180 3550
## 5 Adelie Dream 37.6 19.3 181 3300
## 6 Adelie Dream 36.5 18 182 3150
## 7 Adelie Dream 36 18.5 186 3100
## 8 Adelie Dream 37 16.9 185 3000
## 9 Adelie Dream 36 17.9 190 3450
## 10 Adelie Dream 37.3 17.8 191 3350
## # ℹ 51 more rows
## # ℹ 2 more variables: sex <fct>, year <int>library(tidyverse)
#install.packages("palmerpenguins")
library(palmerpenguins)
#glimpse(penguins)
library(gapminder)
#glimpse(gapminder)Practice
Practice 1
- Create a subset from penguins that only contains Gentoo penguins with a bill depth greater than or equal to 15.5 millimeters.
penguins %>%
filter(species == "Gentoo") %>%
filter(bill_depth_mm >= 15.5)## # A tibble: 40 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Gentoo Biscoe 50 16.3 230 5700
## 2 Gentoo Biscoe 49 16.1 216 5550
## 3 Gentoo Biscoe 49.3 15.7 217 5850
## 4 Gentoo Biscoe 46.3 15.8 215 5050
## 5 Gentoo Biscoe 59.6 17 230 6050
## 6 Gentoo Biscoe 48.4 16.3 220 5400
## 7 Gentoo Biscoe 44.4 17.3 219 5250
## 8 Gentoo Biscoe 48.7 15.7 208 5350
## 9 Gentoo Biscoe 49.6 16 225 5700
## 10 Gentoo Biscoe 50.5 15.9 222 5550
## # ℹ 30 more rows
## # ℹ 2 more variables: sex <fct>, year <int>- Create a subset from penguins that contains observations for male penguins recorded at Dream or Biscoe Islands.
penguins %>%
filter( sex == "male") %>%
filter( island %in% c("Dream","Biscoe"))## # A tibble: 145 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Biscoe 37.7 18.7 180 3600
## 2 Adelie Biscoe 38.2 18.1 185 3950
## 3 Adelie Biscoe 38.8 17.2 180 3800
## 4 Adelie Biscoe 40.6 18.6 183 3550
## 5 Adelie Biscoe 40.5 18.9 180 3950
## 6 Adelie Dream 37.2 18.1 178 3900
## 7 Adelie Dream 40.9 18.9 184 3900
## 8 Adelie Dream 39.2 21.1 196 4150
## 9 Adelie Dream 38.8 20 190 3950
## 10 Adelie Dream 39.8 19.1 184 4650
## # ℹ 135 more rows
## # ℹ 2 more variables: sex <fct>, year <int>- Create a subset from penguins that contains penguins that are either Gentoo or Chinstrap type and have a body mass greater than 4500 g.
penguins %>%
filter( species %in% c("Gentoo","Chinstrap")) %>%
filter( body_mass_g > 4500)## # A tibble: 108 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Gentoo Biscoe 50 16.3 230 5700
## 2 Gentoo Biscoe 50 15.2 218 5700
## 3 Gentoo Biscoe 47.6 14.5 215 5400
## 4 Gentoo Biscoe 46.5 13.5 210 4550
## 5 Gentoo Biscoe 45.4 14.6 211 4800
## 6 Gentoo Biscoe 46.7 15.3 219 5200
## 7 Gentoo Biscoe 46.8 15.4 215 5150
## 8 Gentoo Biscoe 40.9 13.7 214 4650
## 9 Gentoo Biscoe 49 16.1 216 5550
## 10 Gentoo Biscoe 45.5 13.7 214 4650
## # ℹ 98 more rows
## # ℹ 2 more variables: sex <fct>, year <int>Practice 2
1
- Starting with the penguins data, only keep the body_mass_g variable.
penguins %>%
select(body_mass_g)## # A tibble: 344 × 1
## body_mass_g
## <int>
## 1 3750
## 2 3800
## 3 3250
## 4 NA
## 5 3450
## 6 3650
## 7 3625
## 8 4675
## 9 3475
## 10 4250
## # ℹ 334 more rows2
- Starting with the penguins data, keep columns from bill_length_mm to body_mass_g, and year.
penguins %>%
select(bill_length_mm:body_mass_g,year)## # A tibble: 344 × 5
## bill_length_mm bill_depth_mm flipper_length_mm body_mass_g year
## <dbl> <dbl> <int> <int> <int>
## 1 39.1 18.7 181 3750 2007
## 2 39.5 17.4 186 3800 2007
## 3 40.3 18 195 3250 2007
## 4 NA NA NA NA 2007
## 5 36.7 19.3 193 3450 2007
## 6 39.3 20.6 190 3650 2007
## 7 38.9 17.8 181 3625 2007
## 8 39.2 19.6 195 4675 2007
## 9 34.1 18.1 193 3475 2007
## 10 42 20.2 190 4250 2007
## # ℹ 334 more rows3
- Starting with the penguins data, keep all columns except island.
penguins %>%
select(-"island")## # A tibble: 344 × 7
## species bill_length_mm bill_depth_mm flipper_length_mm body_mass_g sex
## <fct> <dbl> <dbl> <int> <int> <fct>
## 1 Adelie 39.1 18.7 181 3750 male
## 2 Adelie 39.5 17.4 186 3800 female
## 3 Adelie 40.3 18 195 3250 female
## 4 Adelie NA NA NA NA <NA>
## 5 Adelie 36.7 19.3 193 3450 female
## 6 Adelie 39.3 20.6 190 3650 male
## 7 Adelie 38.9 17.8 181 3625 female
## 8 Adelie 39.2 19.6 195 4675 male
## 9 Adelie 34.1 18.1 193 3475 <NA>
## 10 Adelie 42 20.2 190 4250 <NA>
## # ℹ 334 more rows
## # ℹ 1 more variable: year <int>4
- From penguins, keep the species column and any columns that end with “mm”.
penguins %>%
select(ends_with("mm"))## # A tibble: 344 × 3
## bill_length_mm bill_depth_mm flipper_length_mm
## <dbl> <dbl> <int>
## 1 39.1 18.7 181
## 2 39.5 17.4 186
## 3 40.3 18 195
## 4 NA NA NA
## 5 36.7 19.3 193
## 6 39.3 20.6 190
## 7 38.9 17.8 181
## 8 39.2 19.6 195
## 9 34.1 18.1 193
## 10 42 20.2 190
## # ℹ 334 more rowsPractice 3
Practice 3A
1a
- In a piped sequence, starting from penguins:
- Only keep observations for female penguins observed on Dream Island, THEN…
penguins %>%
filter(sex == "female") %>%
filter(island == "Dream")## # A tibble: 61 × 8
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Dream 39.5 16.7 178 3250
## 2 Adelie Dream 39.5 17.8 188 3300
## 3 Adelie Dream 36.4 17 195 3325
## 4 Adelie Dream 42.2 18.5 180 3550
## 5 Adelie Dream 37.6 19.3 181 3300
## 6 Adelie Dream 36.5 18 182 3150
## 7 Adelie Dream 36 18.5 186 3100
## 8 Adelie Dream 37 16.9 185 3000
## 9 Adelie Dream 36 17.9 190 3450
## 10 Adelie Dream 37.3 17.8 191 3350
## # ℹ 51 more rows
## # ℹ 2 more variables: sex <fct>, year <int>1b
- Keep variables species, and any variable starting with “bill”
penguins %>%
filter(sex == "female") %>%
filter(island == "Dream") %>%
select(species,starts_with("bill"))## # A tibble: 61 × 3
## species bill_length_mm bill_depth_mm
## <fct> <dbl> <dbl>
## 1 Adelie 39.5 16.7
## 2 Adelie 39.5 17.8
## 3 Adelie 36.4 17
## 4 Adelie 42.2 18.5
## 5 Adelie 37.6 19.3
## 6 Adelie 36.5 18
## 7 Adelie 36 18.5
## 8 Adelie 37 16.9
## 9 Adelie 36 17.9
## 10 Adelie 37.3 17.8
## # ℹ 51 more rows2
- Add a column to penguins that contains a new column flipper_m, which is the flipper_length_mm (flipper length in millimeters) converted to units of meters.
penguins %>%
mutate(flipper_length_m = flipper_length_mm/1000)## # A tibble: 344 × 9
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Torgersen 39.1 18.7 181 3750
## 2 Adelie Torgersen 39.5 17.4 186 3800
## 3 Adelie Torgersen 40.3 18 195 3250
## 4 Adelie Torgersen NA NA NA NA
## 5 Adelie Torgersen 36.7 19.3 193 3450
## 6 Adelie Torgersen 39.3 20.6 190 3650
## 7 Adelie Torgersen 38.9 17.8 181 3625
## 8 Adelie Torgersen 39.2 19.6 195 4675
## 9 Adelie Torgersen 34.1 18.1 193 3475
## 10 Adelie Torgersen 42 20.2 190 4250
## # ℹ 334 more rows
## # ℹ 3 more variables: sex <fct>, year <int>, flipper_length_m <dbl>Practice 3B
1
- The year column in penguins is currently an integer. Add a new column named year_fct that is the year converted to a factor (hint: as.factor())
penguins %>%
mutate(year_fct = as.factor(year))## # A tibble: 344 × 9
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Torgersen 39.1 18.7 181 3750
## 2 Adelie Torgersen 39.5 17.4 186 3800
## 3 Adelie Torgersen 40.3 18 195 3250
## 4 Adelie Torgersen NA NA NA NA
## 5 Adelie Torgersen 36.7 19.3 193 3450
## 6 Adelie Torgersen 39.3 20.6 190 3650
## 7 Adelie Torgersen 38.9 17.8 181 3625
## 8 Adelie Torgersen 39.2 19.6 195 4675
## 9 Adelie Torgersen 34.1 18.1 193 3475
## 10 Adelie Torgersen 42 20.2 190 4250
## # ℹ 334 more rows
## # ℹ 3 more variables: sex <fct>, year <int>, year_fct <fct>2a
- Starting with penguins, do the following within a single mutate() function:
- Convert the species variable to a character
penguins %>%
mutate(species_cha = as.character(species))## # A tibble: 344 × 9
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Torgersen 39.1 18.7 181 3750
## 2 Adelie Torgersen 39.5 17.4 186 3800
## 3 Adelie Torgersen 40.3 18 195 3250
## 4 Adelie Torgersen NA NA NA NA
## 5 Adelie Torgersen 36.7 19.3 193 3450
## 6 Adelie Torgersen 39.3 20.6 190 3650
## 7 Adelie Torgersen 38.9 17.8 181 3625
## 8 Adelie Torgersen 39.2 19.6 195 4675
## 9 Adelie Torgersen 34.1 18.1 193 3475
## 10 Adelie Torgersen 42 20.2 190 4250
## # ℹ 334 more rows
## # ℹ 3 more variables: sex <fct>, year <int>, species_cha <chr>2b
- Add a new column (called flipper_cm with flipper length in centimeters)
penguins %>%
mutate(flipper_length_cm = flipper_length_mm/100)## # A tibble: 344 × 9
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Torgersen 39.1 18.7 181 3750
## 2 Adelie Torgersen 39.5 17.4 186 3800
## 3 Adelie Torgersen 40.3 18 195 3250
## 4 Adelie Torgersen NA NA NA NA
## 5 Adelie Torgersen 36.7 19.3 193 3450
## 6 Adelie Torgersen 39.3 20.6 190 3650
## 7 Adelie Torgersen 38.9 17.8 181 3625
## 8 Adelie Torgersen 39.2 19.6 195 4675
## 9 Adelie Torgersen 34.1 18.1 193 3475
## 10 Adelie Torgersen 42 20.2 190 4250
## # ℹ 334 more rows
## # ℹ 3 more variables: sex <fct>, year <int>, flipper_length_cm <dbl>2c
- Convert the island column to lowercase
penguins %>%
mutate(island_low = tolower(island))## # A tibble: 344 × 9
## species island bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
## <fct> <fct> <dbl> <dbl> <int> <int>
## 1 Adelie Torgersen 39.1 18.7 181 3750
## 2 Adelie Torgersen 39.5 17.4 186 3800
## 3 Adelie Torgersen 40.3 18 195 3250
## 4 Adelie Torgersen NA NA NA NA
## 5 Adelie Torgersen 36.7 19.3 193 3450
## 6 Adelie Torgersen 39.3 20.6 190 3650
## 7 Adelie Torgersen 38.9 17.8 181 3625
## 8 Adelie Torgersen 39.2 19.6 195 4675
## 9 Adelie Torgersen 34.1 18.1 193 3475
## 10 Adelie Torgersen 42 20.2 190 4250
## # ℹ 334 more rows
## # ℹ 3 more variables: sex <fct>, year <int>, island_low <chr>Practice 4
- Starting with penguins, create a summary table containing the maximum and minimum flippers length (call the columns flip_max and flip_min) for Adelie penguins, grouped by island.
penguins %>%
filter(species == "Adelie") %>%
group_by(island) %>%
summarise(
flip_max= max(flipper_length_mm, na.rm = T),
flip_min= min(flipper_length_mm ,na.rm = T)
)## # A tibble: 3 × 3
## island flip_max flip_min
## <fct> <int> <int>
## 1 Biscoe 203 172
## 2 Dream 208 178
## 3 Torgersen 210 176- Starting with penguins, in a piped sequence:
- Add a new column called bill_ratio that is the ratio of bill length to bill depth (hint: mutate())
- Only keep columns species and bill_ratio
- Create a summary table containing the mean of the bill_ratio variable, by species (name the column in the summary table bill_ratio_mean)
penguins %>%
mutate(bill_ratio = bill_length_mm/bill_depth_mm) %>%
select(species,bill_ratio) %>%
group_by(species) %>%
summarise(bill_ratio_mean= mean(bill_ratio, na.rm = T) )## # A tibble: 3 × 2
## species bill_ratio_mean
## <fct> <dbl>
## 1 Adelie 2.12
## 2 Chinstrap 2.65
## 3 Gentoo 3.18Lecture 2
penguins %>%
count(year,species )## # A tibble: 9 × 3
## year species n
## <int> <fct> <int>
## 1 2007 Adelie 50
## 2 2007 Chinstrap 26
## 3 2007 Gentoo 34
## 4 2008 Adelie 50
## 5 2008 Chinstrap 18
## 6 2008 Gentoo 46
## 7 2009 Adelie 52
## 8 2009 Chinstrap 24
## 9 2009 Gentoo 44penguins %>%
count(year,species ) %>%
pivot_wider(names_from = species,values_from = n)## # A tibble: 3 × 4
## year Adelie Chinstrap Gentoo
## <int> <int> <int> <int>
## 1 2007 50 26 34
## 2 2008 50 18 46
## 3 2009 52 24 44penguins %>%
count(year,species ) %>%
pivot_wider(names_from = species,values_from = n) %>%
pivot_longer(cols = Adelie:Gentoo,names_to = "species" ,values_to = "n" )## # A tibble: 9 × 3
## year species n
## <int> <chr> <int>
## 1 2007 Adelie 50
## 2 2007 Chinstrap 26
## 3 2007 Gentoo 34
## 4 2008 Adelie 50
## 5 2008 Chinstrap 18
## 6 2008 Gentoo 46
## 7 2009 Adelie 52
## 8 2009 Chinstrap 24
## 9 2009 Gentoo 44penguins %>%
count(species ,island ,year) %>%
pivot_wider(names_from = c(species,island ),values_from = n)## # A tibble: 3 × 6
## year Adelie_Biscoe Adelie_Dream Adelie_Torgersen Chinstrap_Dream
## <int> <int> <int> <int> <int>
## 1 2007 10 20 20 26
## 2 2008 18 16 16 18
## 3 2009 16 20 16 24
## # ℹ 1 more variable: Gentoo_Biscoe <int>lecture 3
l_f<-function(data,mu,sigma){
n<-length(data)
r <- -(n/2)*log(2*pi* (sigma)^2)-sum((data-mu)^2)/(2*sigma^2)
return(r)
}
data<-c(5,6,7,8,9,15,14)
mu<-10
sigma<-3
l_f(data = data,mu = mu,sigma = sigma)## [1] -19.45619dnorm(x=data,mean=10,sd = 3)## [1] 0.03315905 0.05467002 0.08065691 0.10648267 0.12579441 0.03315905 0.05467002prod(dnorm(x=data,mean=10,sd = 3))## [1] 3.550459e-09#that is like likelihood
log(prod(dnorm(x=data,mean=10,sd = 3)))## [1] -19.45619# this is loglikelihoodAssignment-02: Tidyverse and ggplot Applications.
1.
From the gapminder dataset, filter the data to include only records from the year 2007 and countries in Asia.
library(tidyverse)
library(gapminder)
#glimpse(gapminder)
library(kableExtra)## Warning: package 'kableExtra' was built under R version 4.4.2##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowsq1<-gapminder %>%
filter(continent == "Asia" ,year >= 2007)
kable(head(q1))| country | continent | year | lifeExp | pop | gdpPercap |
|---|---|---|---|---|---|
| Afghanistan | Asia | 2007 | 43.828 | 31889923 | 974.5803 |
| Bahrain | Asia | 2007 | 75.635 | 708573 | 29796.0483 |
| Bangladesh | Asia | 2007 | 64.062 | 150448339 | 1391.2538 |
| Cambodia | Asia | 2007 | 59.723 | 14131858 | 1713.7787 |
| China | Asia | 2007 | 72.961 | 1318683096 | 4959.1149 |
| Hong Kong, China | Asia | 2007 | 82.208 | 6980412 | 39724.9787 |
2.
Using mtcars, create a new dataframe containing only the columns mpg, hp, wt, and gear.
library(kableExtra)
q2<-mtcars %>%
select(mpg, hp, wt, gear)
kable(head(q2))| mpg | hp | wt | gear | |
|---|---|---|---|---|
| Mazda RX4 | 21.0 | 110 | 2.620 | 4 |
| Mazda RX4 Wag | 21.0 | 110 | 2.875 | 4 |
| Datsun 710 | 22.8 | 93 | 2.320 | 4 |
| Hornet 4 Drive | 21.4 | 110 | 3.215 | 3 |
| Hornet Sportabout | 18.7 | 175 | 3.440 | 3 |
| Valiant | 18.1 | 105 | 3.460 | 3 |
3.
Arrange the mtcars data in descending order by mpg and hp.
q3<- mtcars %>%
arrange(desc(mpg),desc(hp))
kable(head(q3))| mpg | cyl | disp | hp | drat | wt | qsec | vs | am | gear | carb | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Toyota Corolla | 33.9 | 4 | 71.1 | 65 | 4.22 | 1.835 | 19.90 | 1 | 1 | 4 | 1 |
| Fiat 128 | 32.4 | 4 | 78.7 | 66 | 4.08 | 2.200 | 19.47 | 1 | 1 | 4 | 1 |
| Lotus Europa | 30.4 | 4 | 95.1 | 113 | 3.77 | 1.513 | 16.90 | 1 | 1 | 5 | 2 |
| Honda Civic | 30.4 | 4 | 75.7 | 52 | 4.93 | 1.615 | 18.52 | 1 | 1 | 4 | 2 |
| Fiat X1-9 | 27.3 | 4 | 79.0 | 66 | 4.08 | 1.935 | 18.90 | 1 | 1 | 4 | 1 |
| Porsche 914-2 | 26.0 | 4 | 120.3 | 91 | 4.43 | 2.140 | 16.70 | 0 | 1 | 5 | 2 |
4.
In gapminder, create a new column gdp_total that multiplies gdpPercap by pop for each observation.
q4<- gapminder %>%
mutate(gdp_total =gdpPercap*pop)
kable(head(q4))| country | continent | year | lifeExp | pop | gdpPercap | gdp_total |
|---|---|---|---|---|---|---|
| Afghanistan | Asia | 1952 | 28.801 | 8425333 | 779.4453 | 6567086330 |
| Afghanistan | Asia | 1957 | 30.332 | 9240934 | 820.8530 | 7585448670 |
| Afghanistan | Asia | 1962 | 31.997 | 10267083 | 853.1007 | 8758855797 |
| Afghanistan | Asia | 1967 | 34.020 | 11537966 | 836.1971 | 9648014150 |
| Afghanistan | Asia | 1972 | 36.088 | 13079460 | 739.9811 | 9678553274 |
| Afghanistan | Asia | 1977 | 38.438 | 14880372 | 786.1134 | 11697659231 |
5.
For the mtcars dataset, calculate the average mpg grouped by the cyl (cylinder) variable.
q5<- mtcars %>%
group_by(cyl) %>%
summarise(ave_mpg= mean(mpg))
kable(head(q5))| cyl | ave_mpg |
|---|---|
| 4 | 26.66364 |
| 6 | 19.74286 |
| 8 | 15.10000 |
6.
Using gapminder, find the mean life expectancy (lifeExp) by continent for the year 2007.
q6<- gapminder %>%
filter(year == 2007) %>%
group_by(continent) %>%
summarise('Mean Life Expectancy'= mean(lifeExp))
kable(head(q6))| continent | Mean Life Expectancy |
|---|---|
| Africa | 54.80604 |
| Americas | 73.60812 |
| Asia | 70.72848 |
| Europe | 77.64860 |
| Oceania | 80.71950 |
7.
Add a new column to mtcars called mpg_level, which is “high” if mpg is above 20 and “low” otherwise.
q7<- mtcars %>%
mutate(mpg_level = case_when(mpg >20 ~ "High",mpg <= 20 ~ "Low"))
kable(head(q7))| mpg | cyl | disp | hp | drat | wt | qsec | vs | am | gear | carb | mpg_level | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Mazda RX4 | 21.0 | 6 | 160 | 110 | 3.90 | 2.620 | 16.46 | 0 | 1 | 4 | 4 | High |
| Mazda RX4 Wag | 21.0 | 6 | 160 | 110 | 3.90 | 2.875 | 17.02 | 0 | 1 | 4 | 4 | High |
| Datsun 710 | 22.8 | 4 | 108 | 93 | 3.85 | 2.320 | 18.61 | 1 | 1 | 4 | 1 | High |
| Hornet 4 Drive | 21.4 | 6 | 258 | 110 | 3.08 | 3.215 | 19.44 | 1 | 0 | 3 | 1 | High |
| Hornet Sportabout | 18.7 | 8 | 360 | 175 | 3.15 | 3.440 | 17.02 | 0 | 0 | 3 | 2 | Low |
| Valiant | 18.1 | 6 | 225 | 105 | 2.76 | 3.460 | 20.22 | 1 | 0 | 3 | 1 | Low |
8.
Count the number of records for each continent in gapminder.
q8<- gapminder %>%
count(continent)
kable(head(q8))| continent | n |
|---|---|
| Africa | 624 |
| Americas | 300 |
| Asia | 396 |
| Europe | 360 |
| Oceania | 24 |
9.
Filter mtcars to only include cars with hp greater than 100, then count how many have gear equal to 4.
q9<- mtcars %>%
filter(hp >100 ) %>%
count(gear)
kable(head(q9))| gear | n |
|---|---|
| 3 | 14 |
| 4 | 5 |
| 5 | 4 |
10.
In gapminder, create a new column lifeExp_category that classifies countries with lifeExp below 50 as “Low”, between 50 and 70 as “Medium,” and above 70 as “High”.
q10<- gapminder %>%
mutate(lifeExp_category= case_when(
lifeExp < 50 ~ "Low",
lifeExp >= 50 & lifeExp <=70 ~ "Medium" ,
lifeExp > 70 ~ "High") )
kable(head(q10))| country | continent | year | lifeExp | pop | gdpPercap | lifeExp_category |
|---|---|---|---|---|---|---|
| Afghanistan | Asia | 1952 | 28.801 | 8425333 | 779.4453 | Low |
| Afghanistan | Asia | 1957 | 30.332 | 9240934 | 820.8530 | Low |
| Afghanistan | Asia | 1962 | 31.997 | 10267083 | 853.1007 | Low |
| Afghanistan | Asia | 1967 | 34.020 | 11537966 | 836.1971 | Low |
| Afghanistan | Asia | 1972 | 36.088 | 13079460 | 739.9811 | Low |
| Afghanistan | Asia | 1977 | 38.438 | 14880372 | 786.1134 | Low |
11.
Create a scatter plot of hp vs wt from mtcars. Label the axes as “Horsepower” and “Weight”.
ggplot(mtcars) +
geom_point(aes(x = hp,y = wt)) +
labs(xlab= "Horsepower", ylab= "Weight",title = 'Scatter plot of WEIGHT & Horsepower of the Cars')
12.
Using gapminder, create a histogram of gdpPercap for the year 2007.
gapminder %>%
filter(year == 2007) %>%
ggplot()+
geom_histogram(aes(x = gdpPercap))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
13.
Create a box plot of lifeExp by continent from gapminder for the year 2007.
gapminder %>%
filter(year == 2007) %>%
ggplot()+
geom_boxplot(mapping = aes(x = continent,y = lifeExp))
14.
Using mtcars, create a bar plot of the count of cars for each gear, filled by cyl.
ggplot(mtcars) +
geom_bar(aes(x = factor(gear),fill = factor(cyl)))
15.
Using gapminder, create a scatter plot of lifeExp vs gdpPercap faceted by continent for the year 2007.
gapminder %>%
filter(year ==2007) %>%
ggplot()+
geom_point(aes(x = lifeExp,y = gdpPercap)) +
facet_wrap(~continent)
16.
Starting with gapminder, filter the data to include only countries with a lifeExp greater than 70 in the year 2000, then select only the columns country, continent, and lifeExp. Use pipes (%>%) to achieve this in a single line.
q16<- gapminder %>%
filter(year > 2000 & lifeExp >70) %>%
select(country, continent,lifeExp)
kable(head(q16))| country | continent | lifeExp |
|---|---|---|
| Albania | Europe | 75.651 |
| Albania | Europe | 76.423 |
| Algeria | Africa | 70.994 |
| Algeria | Africa | 72.301 |
| Argentina | Americas | 74.340 |
| Argentina | Americas | 75.320 |
17.
Using gapminder, create a summary table that shows the minimum, maximum, and average gdpPercap by continent for the year 1997.
q17<- gapminder %>%
filter(year == 1997) %>%
group_by(continent) %>%
summarise(minimum=min(gdpPercap), maximum = max(gdpPercap), average =mean(gdpPercap))
kable(head(q17))| continent | minimum | maximum | average |
|---|---|---|---|
| Africa | 312.1884 | 14722.84 | 2378.760 |
| Americas | 1341.7269 | 35767.43 | 8889.301 |
| Asia | 415.0000 | 40300.62 | 9834.093 |
| Europe | 3193.0546 | 41283.16 | 19076.782 |
| Oceania | 21050.4138 | 26997.94 | 24024.175 |
18.
Using mtcars, filter for cars with cyl equal to 6, then arrange by hp in ascending order and display the first 5 rows.
q18<- mtcars %>%
filter(cyl == 6) %>%
arrange(hp)
kable(head(q18))| mpg | cyl | disp | hp | drat | wt | qsec | vs | am | gear | carb | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Valiant | 18.1 | 6 | 225.0 | 105 | 2.76 | 3.460 | 20.22 | 1 | 0 | 3 | 1 |
| Mazda RX4 | 21.0 | 6 | 160.0 | 110 | 3.90 | 2.620 | 16.46 | 0 | 1 | 4 | 4 |
| Mazda RX4 Wag | 21.0 | 6 | 160.0 | 110 | 3.90 | 2.875 | 17.02 | 0 | 1 | 4 | 4 |
| Hornet 4 Drive | 21.4 | 6 | 258.0 | 110 | 3.08 | 3.215 | 19.44 | 1 | 0 | 3 | 1 |
| Merc 280 | 19.2 | 6 | 167.6 | 123 | 3.92 | 3.440 | 18.30 | 1 | 0 | 4 | 4 |
| Merc 280C | 17.8 | 6 | 167.6 | 123 | 3.92 | 3.440 | 18.90 | 1 | 0 | 4 | 4 |
19.
In gapminder, create a summary that shows the average lifeExp for each continent in 2007, but only include continents where the average lifeExp is above 60.
q19<- gapminder %>%
filter(year ==2007 ) %>%
group_by(continent) %>%
summarise(average_lifeExp=mean(lifeExp)) %>%
filter(average_lifeExp >60 )
kable(head(q19))| continent | average_lifeExp |
|---|---|
| Americas | 73.60812 |
| Asia | 70.72848 |
| Europe | 77.64860 |
| Oceania | 80.71950 |
20.
Add a column to mtcars called weight_class which is “light” if wt is less than 3, “medium” if between 3 and 4, and “heavy” if above 4.
q20<- mtcars %>%
mutate(weight_class= case_when(wt < 3 ~ "light",
wt >=3 & wt <=4 ~ "medium",
wt > 4 ~ "heavy") )
kable(head(q20))| mpg | cyl | disp | hp | drat | wt | qsec | vs | am | gear | carb | weight_class | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Mazda RX4 | 21.0 | 6 | 160 | 110 | 3.90 | 2.620 | 16.46 | 0 | 1 | 4 | 4 | light |
| Mazda RX4 Wag | 21.0 | 6 | 160 | 110 | 3.90 | 2.875 | 17.02 | 0 | 1 | 4 | 4 | light |
| Datsun 710 | 22.8 | 4 | 108 | 93 | 3.85 | 2.320 | 18.61 | 1 | 1 | 4 | 1 | light |
| Hornet 4 Drive | 21.4 | 6 | 258 | 110 | 3.08 | 3.215 | 19.44 | 1 | 0 | 3 | 1 | medium |
| Hornet Sportabout | 18.7 | 8 | 360 | 175 | 3.15 | 3.440 | 17.02 | 0 | 0 | 3 | 2 | medium |
| Valiant | 18.1 | 6 | 225 | 105 | 2.76 | 3.460 | 20.22 | 1 | 0 | 3 | 1 | medium |
21.
In gapminder, group by continent and calculate both the average and the maximum pop (population) for each continent.
q21<- gapminder %>%
group_by(continent) %>%
summarise(Average =mean(pop), Maximum =max(pop))
kable(head(q21))| continent | Average | Maximum |
|---|---|---|
| Africa | 9916003 | 135031164 |
| Americas | 24504795 | 301139947 |
| Asia | 77038722 | 1318683096 |
| Europe | 17169765 | 82400996 |
| Oceania | 8874672 | 20434176 |
22.
Using gapminder, create a column called gdp_level that classifies gdpPercap as “Low” (<1000), “Medium” (1000-10000), and “High” (>10000).
q22<- gapminder %>%
mutate(gdp_level= case_when(
gdpPercap < 1000 ~ "Low",
gdpPercap >= 1000 & gdpPercap <= 10000 ~ "Medium",
gdpPercap > 10000 ~ "High") )
kable(head(q22))| country | continent | year | lifeExp | pop | gdpPercap | gdp_level |
|---|---|---|---|---|---|---|
| Afghanistan | Asia | 1952 | 28.801 | 8425333 | 779.4453 | Low |
| Afghanistan | Asia | 1957 | 30.332 | 9240934 | 820.8530 | Low |
| Afghanistan | Asia | 1962 | 31.997 | 10267083 | 853.1007 | Low |
| Afghanistan | Asia | 1967 | 34.020 | 11537966 | 836.1971 | Low |
| Afghanistan | Asia | 1972 | 36.088 | 13079460 | 739.9811 | Low |
| Afghanistan | Asia | 1977 | 38.438 | 14880372 | 786.1134 | Low |
23.
For gapminder, filter to keep only the data for Japan. Then, use
pivot_wider to create a table where each year is a separate
column and the values represent lifeExp.
q23<- gapminder %>%
filter(country == "Japan") %>%
select(year,lifeExp ) %>%
pivot_wider(names_from = year,values_from = lifeExp)
kable(head(q23))| 1952 | 1957 | 1962 | 1967 | 1972 | 1977 | 1982 | 1987 | 1992 | 1997 | 2002 | 2007 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 63.03 | 65.5 | 68.73 | 71.43 | 73.42 | 75.38 | 77.11 | 78.67 | 79.36 | 80.69 | 82 | 82.603 |
24.
Create a density plot for mpg in mtcars, setting different fill colors for each cyl.
ggplot(mtcars)+
geom_density(aes(x = mpg,fill = factor(cyl)))
25.
Using gapminder, create a line plot showing the lifeExp over the years for the country “India.”
gapminder %>%
filter(country == "India") %>%
ggplot(aes(x = year,y = lifeExp))+
geom_point()+
geom_line()
26.
Create a violin plot for gdpPercap by continent in gapminder, for the year 2007.
gapminder %>%
filter(year == 2007) %>%
ggplot()+
geom_violin(mapping = aes(x = continent ,y = gdpPercap))
27.
Using mtcars, create a scatter plot of hp
vs mpg, where the size of the points represents
wt, and the color represents cyl.
ggplot(mtcars) +
geom_point(aes(x = hp,y = mpg,size = wt, colour = cyl))
28.
Create a faceted plot from gapminder showing gdpPercap vs lifeExp, faceted by continent and year (use only 2002 and 2007).
gapminder %>%
filter(year== 2002 | year == 2007) %>%
ggplot()+
geom_point(aes(x = gdpPercap,y = lifeExp))+
facet_grid(continent ~ year)
29.
Create a scatter plot of lifeExp vs gdpPercap in gapminder for the year 2007. Add text annotations to label the points for countries where lifeExp is above 80.
gapminder %>%
filter(year ==2007) %>%
ggplot(aes(x = gdpPercap,y = lifeExp))+
geom_point()+
geom_text(data = gapminder %>%
filter(year ==2007) %>%
filter(lifeExp > 80),
mapping = aes( label = country))
30.
In gapminder, group the data by both continent and year. For each group, calculate the total population (pop) and average life expectancy (lifeExp). Display only the results for years 1952, 1977, and 2007.
q30<- gapminder %>%
group_by(continent,year) %>%
summarise('Total population'= sum(pop) ,
'Average life expectancy'= mean(lifeExp)) %>%
filter(year %in% c(1952,1977,2007))## `summarise()` has grouped output by 'continent'. You can override using the
## `.groups` argument.kable(q30)| continent | year | Total population | Average life expectancy |
|---|---|---|---|
| Africa | 1952 | 237640501 | 39.13550 |
| Africa | 1977 | 433061021 | 49.58042 |
| Africa | 2007 | 929539692 | 54.80604 |
| Americas | 1952 | 345152446 | 53.27984 |
| Americas | 1977 | 578067699 | 64.39156 |
| Americas | 2007 | 898871184 | 73.60812 |
| Asia | 1952 | 1395357351 | 46.31439 |
| Asia | 1977 | 2384513556 | 59.61056 |
| Asia | 2007 | 3811953827 | 70.72848 |
| Europe | 1952 | 418120846 | 64.40850 |
| Europe | 1977 | 517164531 | 71.93777 |
| Europe | 2007 | 586098529 | 77.64860 |
| Oceania | 1952 | 10686006 | 69.25500 |
| Oceania | 1977 | 17239000 | 72.85500 |
| Oceania | 2007 | 24549947 | 80.71950 |
Lecture 04
Profile Likelihood
Suppose that the following sample is drawn from normal distribution with mean, μ and variance σ2 where both the parameter are unknown.
23.64 31.75 24.93 32.06 25.40 22.06 23.19 21.30 18.07 22.89 24.46 31.05 29.64 24.81 28.48 17.77 27.84 18.33 17.75 19.85
# lik_fun <- function(sigma, data){
# -sum(log(dnorm(data, mu, sigma)))}
plik <- function(mu, data){
sig <- optim(par = 1,
fn = function(sigma, data){
-sum(log(dnorm(data, mu, sigma)))},
data = data)$par
-sum(dnorm(data, mu, sig, log=TRUE))
}
set.seed(25)
x <- round(rnorm(n = 50, mean = 35, sd = 4),2)
optim(15, plik, data = x)## Warning in optim(15, plik, data = x): one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly## $par
## [1] 34.14404
##
## $value
## [1] 138.9369
##
## $counts
## function gradient
## 32 NA
##
## $convergence
## [1] 0
##
## $message
## NULLm <- seq(15, 50, length.out = 250)
mplik <- NULL
for(i in 1:250){
mplik[i] = -plik(m[i], x)
}## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(sigma, data) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directlymu_lik = as.data.frame(cbind(m, mplik)) %>%
filter(mplik == max(mplik)) %>%
select(m) %>%
as.numeric()
mu_lik## [1] 34.11647ggplot()+
geom_line(aes(x = m, y = mplik))+
theme_bw()+
xlab(expression(mu)) + ylab(expression(l(mu))) +
geom_vline(aes(m, mplik), xintercept = mu_lik, color = "red")## Warning: `geom_vline()`: Ignoring `mapping` because `xintercept` was provided.
### Home Work
y<-c(39.91,43.97,23.19,38.87, 39.81, 22.70, 27.72, 44.67, 28.64,38.58,37.38, 49.90,41.48, 41.73, 51.45, 35.72, 33.41, 76.60, 32.02,30.35,38.29, 38.71, 31.39, 39.64,26.49)
plik_gamma<-function(alpha,data){
beta<-optim(par=1,fn = function(beta_1){
-sum(dgamma(data,shape = alpha,scale =beta_1,log = T))}
)$par
-sum(dgamma(data,alpha,beta,log = T))
}
optim(mean(y),plik_gamma,data=y)## Warning in optim(mean(y), plik_gamma, data = y): one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in optim(par = 1, fn = function(beta_1) {: one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly## $par
## [1] 37.88467
##
## $value
## [1] 100.5386
##
## $counts
## function gradient
## 50 NA
##
## $convergence
## [1] 0
##
## $message
## NULLNewtons Method
Example 1
Suppose that a random sample of size 15 has been drawn from N(θ, 1) distribution which yields $ x_i^2 =4434.48$ x2 and $ x_i =257. 59$. Find an approximation of the MLE of θ, \(\hat\theta\). Assume the initial estimate \(\theta_0\) = 10.
theta0 <- 10
eps <- 10
result <- c()
while( eps > 0.001) {
l1 <- 257.59 - 15* theta0
l2 <- -15
theta1 <- theta0 - l1/l2
eps <- abs(theta1 - theta0)
result <- rbind(result, c(theta0, theta1, l1, eps ))
theta0 <- theta1
}
result## [,1] [,2] [,3] [,4]
## [1,] 10.00000 17.17267 107.59 7.172667
## [2,] 17.17267 17.17267 0.00 0.000000Example 2
Suppose that a random sample of size 15 has been drawn from N(0, σ2) distribution which yields \(\sum x_i^2\) = 124.88. Find an approximation of the MLE of \(σ_0^2\) with \(σ_0^2\) = 2.5.
theta0<- 2.5
eps<- 2.5
n<-15
result<- c()
while( eps > .001) {
l1<- -n/(2*theta0)+ 124.88 /(2* theta0^2)
l2<- n/(2*theta0^2)- 124.88 / theta0^3
theta1<- theta0- l1/l2
eps<- abs(theta1- theta0)
result<-rbind(result,c(theta0,theta1,l1,eps))
theta0<-theta1
}
result## [,1] [,2] [,3] [,4]
## [1,] 2.500000 3.529162 6.990400e+00 1.029162e+00
## [2,] 3.529162 4.819141 2.888103e+00 1.289979e+00
## [3,] 4.819141 6.247261 1.132290e+00 1.428120e+00
## [4,] 6.247261 7.495147 3.993398e-01 1.247886e+00
## [5,] 7.495147 8.174777 1.108349e-01 6.796305e-01
## [6,] 8.174777 8.319985 1.689693e-02 1.452075e-01
## [7,] 8.319985 8.325326 5.795058e-04 5.341750e-03
## [8,] 8.325326 8.325333 7.431753e-07 6.868025e-06Example 3
Suppose that a random sample of size 15 has been drawn from N(μ, \(σ^2\)) distribution, where both μ and \(σ^2\) are unknown. If \(\sum x_i\) = 871.67 and $ x_i^2$= 30736.31. Estimate (μ, \(σ^2\)) using Newtons Method. Start with θ0 = (\(μ_0\) = 30, \(σ_0^2\) = 2.5).
dat3 <- read.table("D:/4th year/AST 430 Statistical Computing IX/AST-430 Batch 27/NewtonMethod3.txt")
n <- nrow(dat3)
result <- c()
mu0 <- 30
sigma2_0 <- 2.5
eps <- 5
while (eps > 0.001){
#g matrix
g1 <- sum(dat3 - mu0) / sigma2_0
g2 <- -n/(2*sigma2_0) + sum((dat3 - mu0)^2) / (2 * sigma2_0 ^2)
g <- matrix(c(g1,g2), nrow = 2)
# H matrix
h11 <- -n/ sigma2_0
h12 = h21 <- -sum(dat3-mu0)/(sigma2_0^2)
h22 <- n/(2 * sigma2_0^2) - sum((dat3 - mu0)^2) / (sigma2_0^3)
H <-matrix(c(h11, h12, h21, h22), nrow = 2, byrow = T)
#Theta
theta0 <- matrix(c(mu0, sigma2_0), nrow = 2)
theta1 <- theta0 - solve(H) %*% g
eps <- abs(max(theta1 - theta0))
result <- rbind(result,c(theta1[1,1], theta1[2,1], eps))
mu0 <- theta1[1,1]
sigma2_0 <- theta1[2,1]
}
result ## [,1] [,2] [,3]
## [1,] 37.28061 1.259971 7.280609e+00
## [2,] 35.48041 1.580313 3.203427e-01
## [3,] 35.14595 2.299423 7.191101e-01
## [4,] 34.99273 3.337266 1.037843e+00
## [5,] 34.92093 4.773570 1.436304e+00
## [6,] 34.88811 6.659163 1.885594e+00
## [7,] 34.87397 8.926468 2.267305e+00
## [8,] 34.86858 11.248606 2.322137e+00
## [9,] 34.86697 12.988199 1.739593e+00
## [10,] 34.86669 13.683765 6.955667e-01
## [11,] 34.86668 13.766197 8.243199e-02
## [12,] 34.86668 13.767202 1.005118e-03
## [13,] 34.86668 13.767203 1.467956e-07#Printing the results
result %>%
kable(col.names = c("$\\hat{\\mu}_0$","$\\hat{\\sigma}_0ˆ2$",
"$\\varepsilon$"),
escape = F, booktabs=T,
digits=3, linesep = "")| \(\hat{\mu}_0\) | \(\hat{\sigma}_0ˆ2\) | \(\varepsilon\) |
|---|---|---|
| 37.281 | 1.260 | 7.281 |
| 35.480 | 1.580 | 0.320 |
| 35.146 | 2.299 | 0.719 |
| 34.993 | 3.337 | 1.038 |
| 34.921 | 4.774 | 1.436 |
| 34.888 | 6.659 | 1.886 |
| 34.874 | 8.926 | 2.267 |
| 34.869 | 11.249 | 2.322 |
| 34.867 | 12.988 | 1.740 |
| 34.867 | 13.684 | 0.696 |
| 34.867 | 13.766 | 0.082 |
| 34.867 | 13.767 | 0.001 |
| 34.867 | 13.767 | 0.000 |
Residual maximum likelihood (REML)
library(mvtnorm)## Warning: package 'mvtnorm' was built under R version 4.4.2 set.seed(165)
y <- matrix(
data = rnorm(n = 25, mean = 20, sd = 5),
ncol = 1)
x <- matrix(
data = rep(x = 1, times = 25),
ncol = 1)
A <- diag(x = length(y)) - x %*% solve(t(x) %*% x) %*% t(x)
#residual likelihood
reml_lik <- function(par, y, A){
w <- t(A) %*% y
mu <- rep(x = 0, times = length(w))
cmat <- t(A) %*% (par * diag(length(y))) %*% A
sum(dmvnorm(x = as.numeric(w),
mean = mu, sigma = cmat,log = TRUE))
}
reml_lik(par = 25, y = y, A = A[,-length(y)])## [1] -72.02488optimise(f = reml_lik,interval = c(0, 1e+100),
y = y, A = A[, -length(y)], maximum = TRUE)## $maximum
## [1] 26.98599
##
## $objective
## [1] -71.98891Principal Component Analysis (PCA)
crmat <- matrix(c(
1, 0.505, 0.569, 0.602, 0.621, 0.603,
0.505, 1, 0.422, 0.467, 0.482, 0.450,
0.569, 0.422, 1, 0.926, 0.877, 0.878,
0.602, 0.467, 0.926, 1, 0.874, 0.894,
0.621, 0.482, 0.877, 0.874, 1, 0.937,
0.603, 0.450, 0.878, 0.894, 0.937, 1
), nrow = 6, ncol = 6, byrow = TRUE)
crmat## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 1.000 0.505 0.569 0.602 0.621 0.603
## [2,] 0.505 1.000 0.422 0.467 0.482 0.450
## [3,] 0.569 0.422 1.000 0.926 0.877 0.878
## [4,] 0.602 0.467 0.926 1.000 0.874 0.894
## [5,] 0.621 0.482 0.877 0.874 1.000 0.937
## [6,] 0.603 0.450 0.878 0.894 0.937 1.000n <- 276
m <- 3
# Eigen Values
p <- 6
eg <- eigen(crmat)
eg## eigen() decomposition
## $values
## [1] 4.45644850 0.78240991 0.45842506 0.16883257 0.07908774 0.05479622
##
## $vectors
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] -0.3507933 0.3956532 0.84666989 0.05182494 -0.01451715 -0.02561535
## [2,] -0.2862040 0.8146406 -0.50202982 0.02010032 -0.01468155 -0.04236155
## [3,] -0.4399784 -0.2632446 -0.11051604 0.50466755 -0.59955835 -0.33278818
## [4,] -0.4468851 -0.1972029 -0.09896293 0.47037458 0.59797472 0.41567414
## [5,] -0.4488806 -0.1614667 -0.06586761 -0.54823826 -0.36866442 0.57586240
## [6,] -0.4474933 -0.2134503 -0.06906640 -0.46947138 0.38290503 -0.61838409lmda <- eg$values
lmda## [1] 4.45644850 0.78240991 0.45842506 0.16883257 0.07908774 0.05479622eL <- eg$vectors[, 1:m]
eL## [,1] [,2] [,3]
## [1,] -0.3507933 0.3956532 0.84666989
## [2,] -0.2862040 0.8146406 -0.50202982
## [3,] -0.4399784 -0.2632446 -0.11051604
## [4,] -0.4468851 -0.1972029 -0.09896293
## [5,] -0.4488806 -0.1614667 -0.06586761
## [6,] -0.4474933 -0.2134503 -0.06906640L <- matrix(NA, nrow = p, ncol = m)
for(i in 1:m) {
L[, i] <- sqrt(lmda[i]) * eL[, i]
}
L## [,1] [,2] [,3]
## [1,] -0.7405352 0.3499708 0.57325558
## [2,] -0.6041852 0.7205817 -0.33990980
## [3,] -0.9288076 -0.2328502 -0.07482720
## [4,] -0.9433880 -0.1744338 -0.06700492
## [5,] -0.9476006 -0.1428236 -0.04459705
## [6,] -0.9446719 -0.1888052 -0.04676285# Communalities
cm <- apply(X = L^2,MARGIN = 1, FUN = sum)
cm## [1] 0.9994939 0.9998164 0.9225019 0.9248977 0.9203343 0.9302392# Specific Variance
sv <- 1 - cm
sv## [1] 0.0005060766 0.0001835913 0.0774981111 0.0751022503 0.0796656638
## [6] 0.0697608285# Proportion explained
prop <- diag(t(L) %*% L) / p
prop## [1] 0.74274142 0.13040165 0.07640418# Cumulative proportion
cumsum(prop)## [1] 0.7427414 0.8731431 0.9495472varimax(L)## $loadings
##
## Loadings:
## [,1] [,2] [,3]
## [1,] -0.354 0.244 0.903
## [2,] -0.234 0.949 0.211
## [3,] -0.921 0.166 0.218
## [4,] -0.903 0.214 0.252
## [5,] -0.887 0.229 0.284
## [6,] -0.907 0.192 0.264
##
## [,1] [,2] [,3]
## SS loadings 3.454 1.123 1.120
## Proportion Var 0.576 0.187 0.187
## Cumulative Var 0.576 0.763 0.950
##
## $rotmat
## [,1] [,2] [,3]
## [1,] 0.8546542 -0.3385584 -0.3936299
## [2,] 0.4824845 0.7979214 0.3612896
## [3,] 0.1917680 -0.4986980 0.8452960# factanal performs MLE factor analysis on Covariance matrix
fit <- factanal(covmat = crmat, n.obs = n, factors = m, rotation = "none")
fit##
## Call:
## factanal(factors = m, covmat = crmat, n.obs = n, rotation = "none")
##
## Uniquenesses:
## [1] 0.512 0.317 0.117 0.005 0.019 0.088
##
## Loadings:
## Factor1 Factor2 Factor3
## [1,] 0.622 0.284 0.143
## [2,] 0.484 0.660 0.121
## [3,] 0.937
## [4,] 0.992 -0.105
## [5,] 0.920 0.367
## [6,] 0.926 0.232
##
## Factor1 Factor2 Factor3
## SS loadings 4.187 0.520 0.236
## Proportion Var 0.698 0.087 0.039
## Cumulative Var 0.698 0.784 0.824
##
## The degrees of freedom for the model is 0 and the fit was 8e-04