講義のーと : データ解析のための統計モデリング. 第２回

Size: px

Start display at page:

Download "講義のーと : データ解析のための統計モデリング. 第２回"

よしじろうさかど
5 years ago
Views:

1 Title 講義のーと : データ解析のための統計モデリング Author(s) 久保, 拓弥 Issue Date 2008 Doc URL Type learningobject Note この講義資料は, 著者のホームページードできます Note(URL) Additional There Information are other files related to this item in HUSCAP File Information kubostat2008b.pdf ( 第 2 回 ) Instructions for use Hokkaido University Collection of Scholarly and Aca

2 ( :11 ) 1 ( ) 5 (+2) 2 ( ) kubo@ees.hokudai.ac.jp R &?! 2 2.? 6 3. : 9 4. λ R? R ( )?

3 ( :11 ) 2 : ( ) / 7 1:?? R & 5?! R 3 R 4 data 5 vector object ( ) 50 R 6 > data [1] [26] data vector object 0 1, 2 ( ) vector [1] [26] 1 26

4 ( :11 ) 3 vector object data R data R 7 8 data 7. R > R 50 data R 9 > length(data) # data? [1] 50 > summary(data) Min. 1st Qu. Median Mean 3rd Qu. Max length() data 50 summary() summary(data) Min. Max. data 1st Qu., Median, 3rd Qu. data 25%, 50%, 75% 8. print(data) 9. R # # 10 Median ( ) 10. (Mean) 3.56 ( ) (variance) > var(data) [1] data ( ) (standard deviation; SD) R > sqrt(var(data)) [1] > sd(data) [1]

5 ( :11 ) ? R 11 > summary(as.factor(data)) > table(data) data (histogram) 11. summary() :-) seq( ) R seq(-0.5, 8.5, 1) > hist(data, breaks = seq(-0.5, 8.5, 1)) Histogram of data Frequency data R ( 50 ) 3.56 (probablistic)

6 ( :11 ) 5? ( ) (Poisson distribution) ? > y <- 0:9 > prob <- dpois(y, lambda = 3.56) > plot(y, prob, type = "b", lty = 2) 14. data ( ) prob y > data.frame(y, prob) # y prob y prob y ( y y {0, 1, 2, 3, }) y ( prob) R y prob dpois(y, lambda =

7 ( :11 ) 6 )? ? ( 3.56 ) (data 50 )? hist() > # data, y, prob > hist(data, breaks = seq(-0.5, 8.5, 1)) lines() 50 > lines(y, prob * 50, type = "b", lty = 2) # 50 Histogram of data R Frequency data? ( 3.56 ) 3.56?? ( )?

8 ( :11 ) 7 2.? ? (probablistic density function) p(y λ) = λy exp( λ) y! p(y λ) 19 ( ) λ y ( : 2 y = 2 ) 19. p(y) y λ exp( y) = e y e (e = ) y! y 4! y = 0, 1, 2, 3,, ( ) p(y λ) = 1 λ (λ 0) y=1 λ : λ = = 21 λ > y <- 0:20 > plot(y, dpois(y, lambda = 3.5), type = "b", lty = 2, pch = 21, + ylab = "prob") > lines(y, dpois(y, lambda = 7.7), type = "b", lty = 2, pch = 23) > lines(y, dpois(y, lambda = 15.1), type = "b", lty = 2, pch = 24) > legend("topright", legend = c(3.5, 7.7, 15.1), pch = c(21, 23, 24), + title = "lambda", cex = 0.7) ( )? 21.! 22. plot() lines() 23. R > +

9 ( :11 ) 8 prob lambda y data (y ) 0, 1, 2, y ( ) 24 > mean(data) # data [1] 3.56 > var(data) # data [1] (event)??

10 ( :11 ) (negative binomial distribution) ? R : data > data [1] [26] ( ) 27 (random number) (sample) sampling 27.

11 ( :11 ) 10 = prob y (, ) : ( )? (estimation) ( 3) prob ? y 3: 1. ( ) X ( X ) 2. X X (parameter) data

12 ( :11 ) 11 R R rpois() data > data <- rpois(50, lambda = 3.5) > data [1] [26] > y <- rpois(50, lambda = 3.5); print(y); print(mean(y)) [1] [26] [1] > y <- rpois(50, lambda = 3.5); print(y); print(mean(y)) [1] [26] [1] 3.14 > y <- rpois(50, lambda = 3.5); print(y); print(mean(y)) [1] [18] [35] [1] λ data

13 ( :11 ) data 2. data 3.56 ( ) λ ?? (maximum likelihood estimation) 31 data i y i {y 1 = 2, y 2 = 2, y 3 = 4,, y 49 = 2, y 50 = 3} λ {y i } = {y 1, y 2, y 3,, y 50 } (likelihood) p(y λ) = λy exp( λ) y! 30. λ 31. ( ) λ L(λ {y i }) = (y 1 2 ) (y 2 2 ) (y 50 3 ) L(λ {y i }) = p(y 1 λ) p(y 2 λ) p(y 3 λ) p(y 50 λ) λ 50 = p(y i λ) i=1 = 50 i=1 λ y i exp( λ) y i! λ λ 33. L(λ {y i }) λ (log likelihood) 34 { } 50 y i log L(λ {y i }) = y i log(λ) λ log(k) i=1 λ log L(λ {y i }) λ k=1 ˆλ 34. log(l) L L log(l)

14 ( :11 ) 13 data λ > lambda <- seq(2, 5, 0.1) > likelihood <- function(lambda) sum(dpois(data, lambda, log = TRUE)) > plot( + lambda, + sapply(lambda, likelihood), + type = "l", + xlab = "lambda", + ylab = "log likelihood" + ) log likelihood lambda ( ) λ log L(λ {y i }) λ 50 { yi } = λ 1 = i=1 50 i=1 y i 50 λ log L(λ {y i }) λ ˆλ 50 i=1 y i 50 = 0 ˆλ ˆλ = 50 i=1 y i 50 = y i = 35 {y i } ˆλ 35. (maximum likelihood estimator) {y 1 = 2, y 2 = 2, y 3 = 4,, y 49 = 2, y 50 = 3} y i ˆλ = 3.56 (maximum likelihood estimate) data 3.56

15 ( :11 ) (quasi likelihood) overdispersion 5. (pesudo likelihood) 0, 1, 2, 0 N y y {0, 1, 2,, N} (binomial distribution) ( ) N p(y N, q) = q y (1 q) N y y Bayesian y {0, 1, 2,, N} > y <- 0:10 > plot(y, dbinom(y, 10, prob = 0.1), type = "b", lty = 2, pch = 21, + ylab = "prob") > lines(y, dbinom(y, 10, prob = 0.3), type = "b", lty = 2, pch = 23) > lines(y, dbinom(y, 10, prob = 0.6), type = "b", lty = 2, pch = 24) > legend("topright", legend = c(0.1, 0.3, 0.6), pch = c(21, 23, 24), + title = "q", cex = 0.7) prob y q N y 4 : 11/10 ( ) (GLM) 2

16 ( :11 ) (normal distribution; :, Gaussian distribution) 37 {1.5, 3.2, 7.7, } µ ± σ > y <- seq(-5, 5, 0.1) > plot(y, dnorm(y, mean = 0, sd = 1), type = "l", + ylab = "prob.density") # > lines(y, dnorm(y, mean = 0, sd = 3), lty = 2) # > lines(y, dnorm(y, mean = 3, sd = 1), lty = 3) # 37.? 38. λ = prob.density y µ σ p(y µ, σ) = 1 µ)2 exp( (y ) 2πσ 2 2σ cm 5 cm cm 178 cm 40.

17 ( :11 ) 16 > pnorm(178, mean = 170, sd = 5) - pnorm(177, mean = 170, sd = 5) [1] > dnorm(177.5, mean = 170, sd = 5) * 1.0 # [1] cm > pnorm(178, mean = 170, sd = 5) - pnorm(178, mean = 170, sd = 5) [1] 0 > dnorm(178, mean = 170, sd = 5) [1] {y i } (i = 1, 2,, N) L(y i µ, σ) = = N p(y µ, σ) y i=1 N 1 exp( (y i µ) 2 ) y 2πσ 2 2σ 2 i=1 y ( : 0.1 cm ) 41 p(y µ, σ) ( ) drnorm() log L(y i µ, σ) = N 0.5N log(2πσ 2 i=1 ) i µ) 2 + N y 2σ 2 σ µ 43 log L(y i µ, σ) 43. N y N i=1 (y i µ) 2 N 2σ 2 i=1 (y i µ) 2 (

18 ( :11 ) 17 ) ˆµ σ 7.? 1.? : {, }, : {0.56, 1.33, 12.4, 9.84, } 2.? {0, 1,, N}, {0, 1,, }, [y min, y max ], [, ], 3. ( )?,,, n, : : < : ({0, 1, 2,, N}) : [, + ] : [0, + ] : [0, + ] : [0, 1]

19 ( :11 ) 18 R ( ) rbinom() glm(family = binomial) rbinom() glm(family = binomial) rpois() glm(family = poisson) rnbinom() glm.nb() ( ) rgamma() glm(family = gamma) rnorm() glm(family = gaussian) glm() 3 : (GLM) 1 N p {0, 1, 2,, N} N = 1 {0, 1} :,, / 7 4:

k2 ( :35 ) ( k2) (GLM) web web 1 :

k2 ( :35 ) ( k2) (GLM) web web 1 : 2012 11 01 k2 (2012-10-26 16:35 ) 1 6 2 (2012 11 01 k2) (GLM) kubo@ees.hokudai.ac.jp web http://goo.gl/wijx2 web http://goo.gl/ufq2 1 : 2 2 4 3 7 4 9 5 : 11 5.1................... 13 6 14 6.1......................