uda2008/main.tex 2008/05/

Transcription

1 uda2008/main.tex 2008/05/ EM

2

3 [ ] ˆ R ˆ 1. R 2. ˆ R C Java 2008 R ˆ Mac R Windows PC Mac OK Linux ( PC vmware ˆ 3

4 ˆ ˆ ˆ 2007 [ ] ˆ 2007 ˆ Google ˆ ˆ ˆ ˆ 4

5 [R ] ˆ ˆ ˆ R ˆ R ˆ R S C UNIX AT&T( Lucent Technologies) ( C UNIX ) ˆ CRAN ˆ SAS, SPSS, Mathematica, (matlab) [R ] 5

6 ˆ The R Project for Statistical Computing ˆ Wiki RjpWiki ˆ R PDF R Introduction to R ver Appendix A [R (Windows)] 6

7 ˆ URL /04/03 R R win32.exe ˆ ( R R win32.exe [R (Windows)] ˆ R 7

8 ˆ R R Console GUI Rgui Font MS Gothic Apply Save R 8

9 [R ] (Windows) R (Unix ) OS % R [return] (Windows).RData (Unix ) R > q() [return] Save workspace image? [y/n/c]: y.rdata R R > > a <- 1:10 (1,2,...,10) a > a^2 a [1] > plot(a,a^2) a a^2 9

10 > foo <- function(x) sum(x^2) foo foo(a) [1] 385 for(i in 1:10) {...} i 1,...,10 > x <- rep(0,10); for(i in 1:10) x[i] <- i^2 > x [1]

11 foo2 <- function(x) if(x>0) x else 0 > foo2(3) [1] 3 > foo2(-3) [1] 0 > help(for) for > help("for") > help(":") : > library() > library(mass) MASS > demo() > demo(graphics) > demo(image) [return] emacs ESS emacs (M-x R R ) [ ] 11

12 ˆ Trevor Hastie, Robert Tibshirani, Jerome H. Friedman, The Elements of Statistical Learning: Data Mining, Inference, and Prediction, ˆ R S) W. N. Venables, Brian D. Ripley Modern Applied Statistics with S Springer-Verlag ˆ R R ˆ R R/Rstat.pdf ˆ gakubu html [ ] gakureki-shushou.txt "Gakureki" "Shushou" "Hokkaido" "Aomori" "Iwate" "Miyagi"

13 13

14 14

15 15

16 > dat <- read.table("gakureki-shushou.txt") # > dim(dat) # ( matrix data.frame ) [1] 47 2 > dat[1:5,] # 5 Gakureki Shushou Hokkaido

17 Aomori Iwate Miyagi Akita > plot(dat) # > f <- lm(shushou ~ Gakureki, dat) # > summary(f) # Call: lm(formula = Shushou ~ Gakureki, data = dat) Residuals: Min 1Q Median 3Q Max Coefficients: Estimate Std. Error t value Pr(> t ) (Intercept) < 2e-16 *** Gakureki e-09 *** --- Signif. codes: 0 '***' '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: on 45 degrees of freedom Multiple R-Squared: , Adjusted R-squared: F-statistic: on 1 and 45 DF, p-value: 5.943e-09 > abline(f,col="red") # > plot(dat,type="n"); text(dat,rownames(dat)); abline(f,col="red",lty=2) # ˆ = β 0 + β 1 + ɛ 17

18 ˆ ˆβ 0 = 1.74, ˆβ 1 = ˆ ˆσ 0 = 0.04, ˆσ 1 = ˆ p 0 = , p 1 = ˆ ˆ [ (Windows)] ˆ R Console ˆ ˆ.R ˆ R GUI 18

19 Shushou Shushou Okinawa Saga Fukushima Shimane Yamagata Miyazaki Tottori Fukui Kagoshima Nagano Iwate Nagasaki Kumamoto Kagawa Shiga NiigataGumma Ooita Yamanashi Okayama Aomori Yamaguchi Tochigi Mie Shizuoka Gifu Ibaraki AkitaKochi Wakayama Tokushima Toyama Ehime Ishikawa Aichi Hiroshima Miyagi Hyogo Fukuoka Osaka SaitamaChiba Nara Kyoto Hokkaido Kanagawa Tokyo Gakureki Gakureki 1 19

20 - R Editor ˆ.R 20

21 ˆ ˆ R Console 21

22 [ ] 22

23 ˆ ( ) ˆ ˆ ( ) 23

24 [ (I)] 2004 ˆ (y =1 =0) log (x 369 ) ˆ P (Y = 1 x) P (Y = 0 x) = β 0 + β 1 x β k x k 24

25 [ ] gakureki-rikon-12.txt: Gakureki E09504 ( ) Shushou A05203 Zouka A05201 ( ) Ninzu A06102 ( :person) Kaku A06202 ( ) Tomo F01503 ( ) Tandoku A06205 ( ) 65Sai A ( ) Kfufu A06302 ( ) Ktan A06304 ( ) Konin A06601 Rikon A06602 X2000data.txt 12 > dat <- read.table("gakureki-rikon-12.txt") # > dim(dat) # [1] > pairs(dat,pch=".") # > f <- princomp(dat,cor=t) # > biplot(f) # 25

26 Gakurek Shushou Zouka Ninzu Kaku Tomo Tandoku X65Sai Kfufu Ktan Konin Comp Kagoshima Kochi Kfufu Yamaguchi Ktan Ehime Miyazaki Ooita Wakayama Hokkaido Nagasaki Tandoku Shimane Tokushima Kumamoto Tokyo KagawaHiroshima Fukuoka Rikon Kyoto X65Sai Akita Okayama Kaku Aomori Hyogo Osaka Shushou Tottori Iwate Nagano Yamanashi Mie Saga Gumma Nara Gakureki Tomo Niigata Fukushima Ishikawa Miyagi Kanagawa Toyama Chiba Konin Yamagata Fukui GifuShizuoka Tochigi Ibaraki Aichi Okinawa Zouka Ninzu Saitama Shiga Rikon Comp

27 ˆ vs,,,,, 65,,, ˆ?,,,, ˆ?, ˆ 1 > ## x y > myregplot <- function(x,y,dat) { + e <- formula(paste(y,"~",x)) # y~x + plot(e,dat,type="n") # + text(dat[,x],dat[,y],rownames(dat)) # + f <- lm(e,dat) # + abline(f,col="red",lty=2) # + title(sub=paste(names(f$coef),round(f$coef,4),sep="=",collapse=", ")) + summary(f)$coef # 27

28 + } > myregplot("tomo","shushou",dat) Estimate Std. Error t value Pr(> t ) (Intercept) e-15 Tomo e-05 > myregplot("zouka","shushou",dat[-match(c("okinawa","tokyo"),rownames(dat)),]) Estimate Std. Error t value Pr(> t ) (Intercept) e-51 Zouka e-03 [ ]

29 Shushou Okinawa Saga Fukushima Shimane Miyazaki Tottori Yamaga Fukui Kagoshima Nagano Nagasaki Kumamoto Iwate Kagawa Shiga OoitaOkayama Yamanashi Gumma Niigata Yamaguchi Aomori Ibaraki Tochigi Mie Shizuoka Gifu Wakayama EhimeKochi Aichi Tokushima Ishikawa Akita Toyama Hiroshima Hyogo Miyagi Fukuoka Osaka NaraChiba Saitama KanagawaKyoto Hokkaido Shushou Saga Shimane Fukushima Yamagata Tottori Miyazaki Nagano Fukui Kagoshima Nagasaki IwateKumamoto Kagawa Shiga Ooita Niigata Okayama Yamanashi Gumma Tochigi Yamaguchi Mie Aomori Ibaraki Gifu Shizuoka Kochi Akita Tokushima Wakayama EhimeToyama Ishikawa Aichi Hiroshima Miyagi Hyogo Fukuoka Osaka Nara Chiba Saitama Kyoto Kanagawa Tokyo Hokkaido Tomo (Intercept)=1.0309, Tomo= Zouka (Intercept)=1.4947, Zouka=

30 asahi.com [ ] http: // 47 ( ) R X2000data.txt > X2000.data <- read.table("x2000data.txt") # X2000 > dim(x2000.data) # [1] > dat <- X2000.data[,c("E09504","A05203")] # > names(dat) <- c("gakureki","shushou") # (names colnames > write.table(dat,"test.txt") # X2000item.txt, X2000code.txt, X2000name.txt [ ] txt Windows (shift-jis, CR/LF ) X2000data.txt X2000item.txt Linux Mac nkf 30

31 R > ## cp932 JIS enconding="shift-jis" > X2000.item <- read.table(file("x2000item.txt",encoding="cp932")) > dim(x2000.item) # [1] > X2000.item[1,] # 1 Imi Tani Zenkoku Bunrui A01101 ( ) 100 A 1) > X2000.item[c("E09504","A05203"),] # Gakureki Shushou Imi Tani Zenkoku Bunrui E09504 ( ) E 7) A A 5) [ ] read.table dat data.frame matrix as.matrix > dat[1:3,] # Gakureki Shushou Hokkaido Aomori Iwate > datmat <- as.matrix(dat) # > datmat[1:3,] # Gakureki Shushou Hokkaido

32 Aomori Iwate > c(is.matrix(dat),is.data.frame(dat)) # [1] FALSE TRUE > c(is.matrix(datmat),is.data.frame(datmat)) # [1] TRUE FALSE > t(datmat) %*% datmat # t() %*% Gakureki Shushou Gakureki Shushou > colnames(datmat) <- c("aaa","bbb") # names > datmat[1:3,] # AAA BBB Hokkaido Aomori Iwate [ (II)] 2003 DNA ˆ (DNA ) 32

33 ˆ ˆ [ ] ˆ 12 ˆ dist() 33

34 ˆ hclust() ˆ plot() > dat <- read.table("gakureki-rikon-12.txt") # > x <- scale(dat) # > h1 <- hclust(dist(x)) # > plot(h1) # > h2 <- hclust(dist(t(x))) # > plot(h2) # [ (III)] 2004 Assessing the uncertainty in hierarchical cluster analysis via multiscale bootstrap resampling ( ) ˆ (p=73 n=916 ) 34

35 Cluster Dendrogram Cluster Dendrogram Height Nagasaki Miyazaki Wakayama Ooita Yamaguchi Ehime Kochi Kagoshima Shimane Akita Iwate Niigata Yamagata Toyama Fukui Tottori Fukushima Saga Gumma Yamanashi Mie Ishikawa Nagano Aomori Okayama Kagawa Tokushima Kumamoto Nara Miyagi Shiga Gifu Ibaraki Tochigi Shizuoka Hokkaido Kyoto Fukuoka Hyogo Hiroshima Aichi Saitama Chiba Kanagawa Osaka Okinawa Tokyo Height Ninzu Shushou Tomo X65Sai Kaku Rikon Gakureki Zouka Konin Tandoku Kfufu Ktan dist(x) hclust (*, "complete") dist(t(x)) hclust (*, "complete") 4 35

36 ˆ (pvclust R ) 36

37 1 1.1 X x f(x) g(x) [ 1.1] Y = g(x) E (g(x)) = g(x)f(x) dx (1.1) [ 1.1] a, b R, g 1 (x), g 2 (x) E(ag 1 (X) + bg 2 (X)) = ae(g 1 (X)) + be(g 2 (X)) 37

38 [ 1.2] [ 1.3] [ 1.4] X E(X) (1.1) g(x) X V (X) (1.1) g(x) (a, b), a < b X P (a < X < b) (1.1) g(x) A I(A) = 1 I(A) = 0 (indicator function) [ 1.5] F (x) = x f(s) ds (1.1) x g(s) [ 1.1] (0, 1) X U(0, 1) f(x) = 1, 0 < x < 1, f(x) = 0 X E(X) = xf(x) dx = 1 0 x dx = 1 2 V (X) = (x E(X)) 2 f(x) dx = 1 0 (x 1 2 )2 dx = 1 12 F (x) = x 0 ds = x, 0 < x < 1; F (x) = 0, x 0; F (x) = 1, x 1 38

39 [ 1.2] m X x f(x) g(x) E(g(X)) = g(x)f(x) dx R m X E(X) V (X) E(X 1 ) C(X 1, X 1 ) C(X 1, X m ) X =, E(X) =, V (X) =.. X 1. X m. E(X m ) C(X m, X 1 ) C(X m, X m ) C(X i, X j ) = E [(X i E(X i ))(X j E(X j ))] X i X j ρ(x i, X j ) = C(X i, X j )/ C(X i, X i )C(X j, X j ) X i V (X i ) [ 1.6] X m w R m Y = w X A A E(Y ) = w E(X), V (Y ) = w V (X)w [ 1.2] X m E(X i ) = µ V (X i ) = σ 2 39

40 0 Y = m i=1 X i/m w = ( 1 m,..., 1 m ) E(Y ) = µ, V (Y ) = σ2 m > sx <- 0.2 # SD(X) > m <- 1:10 # m=1,2,...,10 > sy <- sx/sqrt(m) # SD(Y) > plot(m,sy) [ 1.7] X m w R m X i i Y = w X w E(Y ) µ 1 E(Y ) = µ, m w i = 1 i=1 V (Y ) w = Σ 1 A ( A Σ 1 A ) 1 [ ] µ 1 40

41 Σ = V (X) 1 m = (1,..., 1) 1 m A = (E(X), 1 m ) m 2 [ 1.3] m E(X i ) V (X i ) X i X j ρ > ex <- c(0.2,0.15,0.1,0.05) # E(X) > sx <- c(0.2,0.2,0.1,0.05) # SD(X) > rho <- 0.3 Σ Σ 1 A C = Σ 1 A ( A Σ 1 A ) 1 > m <- length(sx) # = 4 > V <- (sx %o% sx) * (diag(m)*(1-rho) + matrix(rho,m,m)) # Sigma > B <- solve(v) # Sigma^(-1) > A <- cbind(ex,1) # A > C <- B %*% A %*% solve(t(a) %*% B %*% A) # C E(Y ) = 0.15 V (Y ) V (Y ) > w <- C %*% c(0.15,1) # > t(w) # [,1] [,2] [,3] [,4] [1,] > t(w) %*% ex # 0.15 [,1] [1,] 0.15 > sqrt(t(w) %*% V %*% w) # SD(Y) [,1] [1,]

42 E(Y ) V (Y ) - > # E(Y) SD(Y) > mysy <- function(ey) { + w <- C %*% c(ey,1) + sqrt(t(w) %*% V %*% w) + } > ey <- seq(-0.1,0.4,length=100) # E(Y) > sy <- sapply(ey,mysy) # SD(Y) > plot(sy,ey,type="l") > points(sx,ex,col="red") [ 1.8] x h(x) 0 E(h(X)) a > 0 E(h(X)) ap (h(x) a) ɛ > 0 E(X) = µ V (X) = σ 2 P ( X µ ɛ) σ2 ɛ 2 42

43 sy ey m sy 5-43

44 [ ] 1.1 f(x) X s 0, s 1,... p(x) E (g(x)) = g(s i )p(s i ) (1.2) i=0 δ(x) f(x) = p(s i )δ(s i ) i=0 (1.1) f A R f(a) = P (X A) E(g(X)) = g(x)f(dx) (1.1) R 44

45 1.2 Y 1, Y 2,... Y n, n = 1, 2,... [ 1.3] Y n Y (convergence in probability) ɛ > 0 lim P ( Y n Y > ɛ) = 0 (1.3) n Y n p Y Y Zn = Y n Y ɛ > 0; lim n P ( Z n > ɛ) = 0 [ ] ˆ, almost surely convergent ( ) P lim Y n(ω) = Y (ω) = 1 n ˆ (convergence in distribution) y; lim n P (Y n y) = P (Y y) 45

46 ˆ ˆ Levy (1937) X k, k = 1, 2,... Y n = n k=1 X k, n = 1, 2,... [ 1.1] X n E(X 1 ) = µ n X n = 1 n X n n p X n µ (1.4) n i=1 X i (sample size) [ ] σ 2 = V (X 1 ) V ( X n ) = σ 2 /n P ( X n µ ɛ ) σ2 nɛ 2 46

47 n 0 [ 1.9] X n E(g(X 1 )) Ȳ n = 1 n n g(x i ) (1.5) i=1 E(Ȳn) = E(g(X 1 )), Ȳ n p E(g(X1 )) n (1.1) (1.5) V (g(x 1 )) V (Ȳn) = V (g(x 1 ))/n [ 1.4] X n (0, 1) X n U(0, 1) X 1,..., X n ( ) x 1,..., x n > ## > x <- runif(10000) # U(0,1) > x[1:5] # [1] > ## > mean(x[1:10]) # 10 [1] > mean(x[1:100]) # 100 [1] > mean(x[1:1000]) # 1000 [1]

48 > mean(x) # (0.5) [1] > y1 <- cumsum(x)/(1:10000) # n > plot(y1,log="x"); abline(h=0.5) > ## > mean((x - mean(x))^2) # (1/12 = 0.833) [1] > ## x 0.1 > mean(x < 0.1) # P(X<0.1)=0.1 [1] > z <- x < 0.1 # 0 1 > as.numeric(z[1:100]) # 100 [1] [38] [75] > cumsum(z[1:100]) # [1] [38] [75] > y2 <- cumsum(z)/(1:10000) # n > plot(y2,log="x"); abline(h=0.1) 48

49 y y Index 図 6 左 Y n = 100 Index 1 n Pn 右 Y n = i=1 Xi 49 1 n Pn i=1 I(Xi < 0.1)

50 [ 1.4] x 1,..., x n ˆF n (x) = 1 n n i=1 I(x i x) = #(x i x) n (1.6) ˆf n (x) = 1 n n δ(x x i ) (1.7) i=1 x i 1/n [ 1.10] x ˆF n (x) X 1,..., X n ˆF n (x) p F (x) n ˆF n (x) ˆF n (x) [ ] x ˆF n (x) p F (x) Glivenko- Cantelli P { } lim sup ˆF n (x) F (x) = 0 n x R = 1 g(x) 50

51 [ 1.5] p = F (q) q = F 1 (p) = inf F (q) p q (inf F ) q F 100p F 1 (1 p) 100p 1 ˆF n (p) 100p p = 0.5 (median) help(quantile) R [ 1.5] 1.4 x 1,..., x n > hist(x,prob=true) # > y <- sort(x) # > y[5000] # (p=0.5) [1] > (y[5000]+y[5001])/2 # (p=0.5) [1] > quantile(x,p=0.5) # median(x) 50% > p <- (1:10000)/10000 # , ,..., , 1 > plot(y,p,pch=".") # [ 1.6] 0 1 X N(0, 1) f(x) = 1 2π exp( x2 2 ) rnorm X 1,..., X N(0, 1) 5% 2.5% ˆF n F 51

52 Histogram of x Density p x y 7 52

53 > x <- rnorm(10000) # N(0,1) > mean(x) # [1] > mean((x-mean(x))^2) # [1] > quantile(x,p=c(0.5,0.95,0.975)) # 5% 2.5% 50% 95% 97.5% > qnorm(c(0.5,0.95,0.975)) # [1] > ### > library(mass) # MASS > truehist(x) # > x0 <- seq(min(x),max(x),length=10000) # x 1000 > lines(x0,dnorm(x0),col=2) # > p <- (1:10000)/10000 # , ,..., , 1 > plot(sort(x),p,pch=".") # > lines(x0,pnorm(x0),col=2) # [ 1.7] library(mass) truehist nbis Scott (1979) > truehist(x) # nbins="scott" > nclass.scott # nbins function (x) { h <- 3.5 * sqrt(stats::var(x)) * length(x)^(-1/3) ceiling(diff(range(x))/h) } 53

54 p x sort(x) 8 54

55 <environment: namespace:grdevices> > truehist(x,nbins=nclass.scott(x)*2) # 2 > truehist(x,nbins=nclass.scott(x)/2) # [ 1.11] Y = (a, b], a < b X 1,..., X n C X F (x) C n(b a) µ σ2 [ 1.12] f(x) 2 h (i) x (a, b) f(x) Y b a h ( ) h E(Y ) f(x) = f (a) + (x a) + O(h 2 ), V (Y ) = f(a) 2 nh + O(n 1 ) x 2 MSE(x) = V (Y ) + (E(Y ) f(x)) 2 = f (a) 2 ( x a + b 2 (ii) (a, b) 1 b a b a MSE(x) dx f (a) 2 h f(a) nh ) 2 + f(a) nh + O ( h 3 + n 1) (iii) 2 MSE(x) dx [f (x) 2 h f(x) ] dx = h2 f (x) 2 dx + 1 nh 12 nh 55

56 x x

57 (iv) h h = [ n 6 ] 1/3 f (x) 2 dx (v) f(x) σ 2 f (x) 2 dx = 1 4 πσ 3 h = (24 π) 1/3 σn 1/3 3.49σn 1/3 Scott (1979) 1.3 (0, 1) U 1, U 2,... U(0, 1) F (x) f(x) X 1, X 2,... (1.1) (1.5) [ 1.13] X = F 1 (U) X F (x) [ 1.8] f(x) = e x, F (x) = 1 e x F 1 (p) = log(1 p) 57

58 > ## X ~ exp(-x) > ## rexp > u <- runif(10000) # U(0,1) > x <- -log(u) # X = F^{-1}(U) > truehist(x) # > x0 <- seq(min(x),max(x),length=1000) > lines(x0,exp(-x0),col=2) # [ 1.14] f(x), g(x) U U(0, 1) V g(v) X f(x) 1. x f(x) cg(x) c 2. U U(0, 1) V g(v) 3. U f(v ) cg(v ) X = V 2 (rejection method) [ 1.9] X N(0, 1) f(x) = 1 2π e x2 2, g(x) = 1 2 e x, c = f(x) cg(x) f(x)/(cg(x)) = e ( x 1)2 2 X 58 2e π

59 > ## X ~ N(0,1) > ## rnorm > u0 <- 2*runif(10000)-1 # U(-1,1) > u <- abs(u0) # U(0,1) > s <- sign(u0) # +1,-1 > v <- -log(runif(10000)) # v ~ exp(-v) > a <- u <= exp(-(v-1)^2/2) # a > x <- (s*v)[a] # > length(x) # [1] 7575 > truehist(x) # > x0 <- seq(min(x),max(x),length=10000) # x 1000 > lines(x0,dnorm(x0),col=2) # [ 1.15] m X N m (µ, Σ) ( f(x µ, Σ) = (2π) m/2 Σ 1/2 exp 1 ) 2 (x µ) Σ 1 (x µ) Σ m m A Σ = AA A A A A m Z N m (0, I m ) X = AZ + µ 59

60 x x 10 60

61 [ 1.10] X N 2 (µ, Σ) [ ] [ ] µ1 1 ρ µ =, Σ = ρ 1 µ 1 = 5, µ 2 = 0, ρ = 0.5 > S <- matrix(c(1,0.5,0.5,1),2,2) # Sigma > S # [,1] [,2] [1,] [2,] > A <- t(chol(s)) # > A # [,1] [,2] [1,] [2,] > z <- matrix(rnorm(2*1000),2,1000) # N(0,1) 2x1000 > var(t(z)) # z [,1] [,2] [1,] [2,] > x <- A %*% z + c(5,0) # z 1000 > var(t(x)) # x [,1] [,2] [1,] [2,] > plot(z[1,],z[2,]) # z µ 2 61

62 > plot(x[1,],x[2,]) # x [ 1.16] X 1, X 2,... p(x t+1 x t ) f(x) p(x t+1 x t )f(x t ) = p(x t x t+1 )f(x t+1 ) X t f(x t ) X t+1 p(x t+1 x t )f(x t ) dx t = p(x t x t+1 )f(x t+1 ) dx t = f(x t+1 ) f(x) 1. q(v x) X 1 t = 1 2. U U(0, 1) V q(v X t ) 3. α(v, X t ) = min { 1, q(x } t V )f(v ) q(v X t )f(x t ) 4. U α(v, X t ) X t+1 = V X t+1 = X t 5. t 2 [ ] x x Metrololis-Hastings f(x) 62

63 z[2, ] x[2, ] z[1, ] x[1, ] z 1000 x 63

64 f(v)/f(x) (Markov Chain Monte Carlo MCMC) q(v x) (proposal distribution) (1.1) (1.5) X 0 N N = 100, 000 X 1, X 2,..., X N M < N M = 10, 000 X M+1, X M+2,..., X N N t=m+1 X t/(n M) (1.1) M burn-in ( ) [ 1.11] 1.10 MCMC MCMC 1.15 MCMC q(v x) { (2d) 2 v[1] x[1] d, v[2] x[2] d q(v x) = 0 ( d, d) q(v x) = q(x v) { α(v, x) = min 1, f(v) } f(x) 64

65 v x > N < # > rho <- 0.5 # > S <- matrix(c(1,rho,rho,1),2,2) # Sigma > mu <- c(5,0) # mu > Sinv <- solve(s) # S > myf <- function(x) exp(-0.5*(x-mu)%*%sinv%*%(x-mu)) # f(x) > d <- 0.5 # > x <- c(0,0) # > xs <- matrix(0,2,n) # > fs <- rep(0,n) # f(x) > cnt <- 0 # > for(t in 1:N) { + v <- x + (2*runif(2)-1)*d # v ~ q(v x) + a <- myf(v)/myf(x) # alpha=min(1,a) + u <- runif(1) # u ~ U(0,1) + if(u <= a) { + x <- v # v + cnt <- cnt + 1 # + } + xs[,t] <- x # x + fs[t] <- myf(x) # f(x) + } > cnt/n # [1] > plot(xs[1,],xs[2,]) # x > segments(c(0,xs[1,-n]),c(0,xs[2,-n]),xs[1,],xs[2,],col="pink") # > plot(log(fs),type="l") # log(f(x)) 65

66 xs[2, ] log(fs) xs[1, ] Index X t log(f(x t )) 66

67 [ 1.17] K p k (x t+1 x t ), k = 1,..., K f(x) X 1, X 2,... π k p k (x t+1 x t ) K k=1 π k = 1 f(x) [ ] 1.17 k = 1, 2,..., K, 1, 2,... p k f(x) [ 1.18] m X 1, X 2,... X f(x) m x i x[i] x x[i] m 1 x[ i] x[ i] h i (x[ i]) = f(x) dx[i] x[ i] x[i] f i (x[i] x[ i]) = f(x) h i (x[ i]) 1. X 1 t = 1 67

68 2. i 3. V f i (V X t [ i]) 4. X t+1 [i] = V X t+1 [ i] = X[ i] 5. t 2 [ ] Gibbs sampler) α = 1 [ 1.19] m X N m (µ, Σ) Y X Y m + 1 ([ X E Y ]) = [ µ, V λ] ([ X Y ]) = [ Σ ] b b c x Y Y x N ( λ + b Σ 1 (x µ), c b Σ 1 b ) Σ e = (c b Σ 1 b) 1, d = eσ 1 b, A = eσ 1 bb Σ 1 ] 1 ] = [ Σ b b c [ Σ 1 + A d 68 d e

69 [ 1.12] X[2] x[1] N ( µ 2 + ρ(x[1] µ 1 ), 1 ρ 2) > N < # > rho <- 0.5 # > S <- matrix(c(1,rho,rho,1),2,2) # Sigma > mu <- c(5,0) # mu > Sinv <- solve(s) # S > myf <- function(x) exp(-0.5*(x-mu)%*%sinv%*%(x-mu)) # f(x) > x <- c(0,0) # > xs <- matrix(0,2,n) # > fs <- rep(0,n) # f(x) > for(t in 1:N) { + i <- floor(runif(1)*2)+1 # + v <- mu[i]+rho*(x[-i]-mu[-i]) + rnorm(1)*sqrt(1-rho^2) # X[i] x[-i] + x[i] <- v # + xs[,t] <- x # x + fs[t] <- myf(x) # f(x) + } > plot(xs[1,],xs[2,]) # x > segments(c(0,xs[1,-n]),c(0,xs[2,-n]),xs[1,],xs[2,],col="pink") # > plot(log(fs),type="l") # log(f(x)) 69

70 xs[2, ] log(fs) xs[1, ] Index X t log(f(x t )) 70

71 [ 1.13] x[i], i = (1, 1),..., (m, m) {+1, 1} H(x) = γ (i,j) x[i]x[j] (Ising model) γ > 0 f(x) exp( H(x)) h[i] = γ j:(i,j) x[j] i j f i(v x[ i]) exp(vh[i]) f i (v x[ i]) = exp(vh[i]) exp(h[i]) + exp( h[i]) i +1, 1 i i +1 cf. > N < # > m <- 50 # m*m > x0 <- matrix(runif(m*m)>0.5,m,m)*2-1 # +1,-1 > gamma <- 1.0 # gamma > 1/2.269 > x <- x0 # x > xi <- function(i) x[(i[1]-1)%%m+1,(i[2]-1)%%m+1] # x[(i[1],i[2])] > d1 <- c(1,0); d2 <- c(-1,0); d3 <- c(0,1); d4 <- c(0,-1) # > for(t in 1:N) { # + i <- trunc(runif(2)*m)+1 # i + a <- exp((xi(i+d1)+xi(i+d2)+xi(i+d3)+xi(i+d4))*gamma) # a=exp(h[i]) + p <- a/(a+(1/a)) # p=f_{i}(+1 x[-i]) + if(runif(1)<=p) v <- 1 else v <- -1 # p v=+1, (1-p) v=-1 71

72 + x[i[1],i[2]] <- v # x[i] := v + } > bw <- rev(gray((0:64)/64)) # 64 > image(x0,axes=f,col=bw) # > image(x,axes=f,col=bw) # N 1.4 [ 1.20] Ω i=1 D i = Ω A D i P (D i A) = P (A D i )P (D i ) j=1 P (A D j)p (D j ) P (A) > 0 [ 1.21] X, Y f(x, y) f f(x) = f(x, y) dy, f(y) = f(x, y) dx f(x y) f(y x) f(x y) = f(y x)f(x) f(y x)f(x) dx f(y) > 0 72

73 14 100,000 73

74 [ ] a X b Y f(x, y) f(y) > 0 f(x y) = f(y x)f(x) R a f(y x)f(x) dx f(x) f(x y) [ 1.22] X, Y Y x x Y X N(0, a 2 ), Y x N(x, b 2 ) X y N ( cy, cb 2), c = b2 a 2 c (i) a = 3, b = 1 (ii) a = 1, b = 3 [ 1.23] µ N(0, τ 2 ) µ n µ σ 2 x 1,..., x n N(µ, σ 2 ) x = (x x n )/n µ τ 2, σ 2 [ 1.24] x[i] y[i], i = (1, 1),..., (m, m) {+1, 1} x y x ɛ (0 < ɛ < 1) 74

75 ( ) { 1 ɛ y[i] = x[i] f(y[i] x[i]) = ɛ y[i] x[i] (i) x y ( f(y x) exp λ i ) x[i]y[i] λ = ɛ log( ɛ ) i x, y (ii) x ( f(x) exp γ ) x[i]x[j] (i,j) γ > 0 (i,j) x x y x ( f(x y) exp λ x[i]y[i] + γ ) x[i]x[j] (1.8) i (i,j) x y [ 1.14] 1.24 (ii) f(x y) X 1, X 2,... 75

76 y x 1.13 H(x) = γ (i,j) x[i]x[j] λ i x[i]y[i] y i λy[i] f(x y) exp( H(x)) h[i] = γ j:(i,j) x[j] + λy[i] x burn-in E(X y) MAP > ## > N < # > M < # burn-in > ns < # > m <- 50 # m*m > gamma <- 1.0 # gamma > 1/2.269 > eps < # 35% (50% ) > lambda <- 0.5*log((1-eps)/eps) # > d1 <- c(1,0); d2 <- c(-1,0); d3 <- c(0,1); d4 <- c(0,-1) # > xi <- function(i) if(i[1]>=1 && i[2]>=1 && i[1]<=m && i[2]<=m) + x[i[1],i[2]] else -1 # x[(i[1],i[2])] > ## ( ) > y0 <- matrix(-1,m,m) # y[i]=-1 > for(i1 in 1:m) for(i2 in 1:m) + if(!(abs(i1-0.5*m)<0.15*m && abs(i2-0.5*m)<0.15*m) && + sqrt((i1-0.5*m)^2+(i2-0.5*m)^2) < 0.4*m ) y0[i1,i2] <- +1 > ## 76

77 > y <- y0 # > i <- runif(m*m) < eps # > y[i] <- - y[i] # > ## > myabserr <- function(x) mean(abs(x-y0)/2) # y0 > mylogf <- function(x) { # log(f(x y)) y0 + s <- 0 + for(i1 in 1:m) for(i2 in 1:m) { + i <- c(i1,i2) + s <- s + xi(i)*(xi(i+d1)+xi(i+d2)+xi(i+d3)+xi(i+d4)) + } + lambda*sum(x*y)+gamma*s/2 + } > ## > x <- y # x > myabserr(x) # [1] 0.35 > abserrs <- rep(0,1+n/ns); logfs <- rep(0,1+n/ns) # abserr logf > abserrs[1] <- myabserr(x); logfs[1] <- mylogf(x) # > xp <- matrix(0,m,m) # > ## > for(t in 1:N) { + i <- trunc(runif(2)*m)+1 # i + a <- exp((xi(i+d1)+xi(i+d2)+xi(i+d3)+xi(i+d4))*gamma + +y[i[1],i[2]]*lambda) # a=exp(h[i]) + p <- a/(a+(1/a)) # p=f_{i}(+1 x[-i],y) + if(runif(1)<=p) v <- 1 else v <- -1 # p v=+1, (1-p) v=-1 77

78 + x[i[1],i[2]] <- v # x[i] := v + if(t>m) xp <- xp + x # x(t) + if((t %% ns) == 0) { # + j <- 1 + t %/% ns + abserrs[j] <- myabserr(x); logfs[j] <- mylogf(x) # + } + } > xn <- x # > myabserr(xn) # [1] > xp <- xp / (N-M) # > myabserr(xp) # [1] > xp2 <- sign(xp) # 2 > myabserr(xp2) # [1] > bw <- rev(gray((0:64)/64)) # 64 > image(y0,axes=f,col=bw) # > image(y,axes=f,col=bw) # > image(xp,axes=f,col=bw) # N-M > image(xp2,axes=f,col=bw) # 2 > plot(0:(n/ns),abserrs,xlab=paste("t/",ns,sep="")) # > abline(v=m/ns,lty=2) # burn-in > plot(0:(n/ns),logfs,xlab=paste("t/",ns,sep="")) # log(f(y x)) > abline(v=m/ns,lty=2) # burn-in [ 1.25] 1.14 ɛ =

79 15 79

80 16 burn-in 2 80

81 abserrs logfs t/2500 t/ P i X t[i] x[i] /(2m 2 ) log f(x t y) (1.8) exp 81

82 > mean(xp2!= y) [1] [ 1.26] 1.14 (i) (ii) (iii) γ λ N M [ 1.15] > nf <- 10 # > d0 <- c(0,0) # > d1 <- c(1,0); d2 <- c(-1,0); d3 <- c(0,1); d4 <- c(0,-1) # > d5 <- c(1,1); d6 <- c(-1,1); d7 <- c(-1,-1); d8 <- c(1,-1) # > xi9 <- function(i) c(xi(i+d0),xi(i+d1),xi(i+d2),xi(i+d3),xi(i+d4), + xi(i+d5),xi(i+d6),xi(i+d7),xi(i+d8)) # 9 > xf <- y # > abserrs2 <- rep(0,nf+1); abserrs2[1] <- myabserr(y) # > ## > for(t in 1:nf) { + x <- xf + for(i1 in 1:m) for(i2 in 1:m) xf[i1,i2] <- median(xi9(c(i1,i2))) + abserrs2[t+1] <- myabserr(xf) + } 82

83 > abserrs2[nf+1] # [1] > plot(0:nf,abserrs2,xlab="iteration") # > image(xf,axes=f,col=bw) # 1.5 [ 1.6] t R X (moment generating function) M X (t) = E(e tx ) [ 1.27] k M X (t) E(X k ) = dm M X (t) t=0 dt k [ 1.28] X N(0, 1) M X (t) E(X k ) k = 2, 4, 6, 8 [ 1.29] [ 1.30] X Y M X+Y (t) = M X (t)m Y (t) X 1,..., X n N(0, 1) X

84 abserrs iteration 18 84

85 Y n = X X 2 n M Y n (t) n Y χ 2 n f(y) = 1 2 n e y 2 y n 2 1, y 0 2 Γ( n 2 ) Y M Y (t) M Yn (t) [ 1.7] t R X (characteristic function) ϕ X (t) = E(e itx ) [ ] ϕ X (t) = M X (it) i = 1 ϕ X (t) f f(x) ϕ X (t) f(t) = Ff(t) = e itx f(x) dx, f(x) = F 1 f(x) = 1 2π e itx f(t) dt ϕ X (t) = f( t) M X (t) ϕ X (t) [ 1.8] X 1 z 1 X G X (z) = E(z X ) 85

86 [ ] G X (e t ) = M X (t) X G X (z) M X (t) [ 1.31] X n > 0 0 < p < 1 Y λ > 0 0 x n, y 0 p X (x) = n! x!(n x)! px (1 p) n x, p Y (y) = e λ λ y y! M X (t) = ( n, ( ) pe t + (1 p)) MY (t) = exp λ(e t 1) [ 1.9] m X = (X 1,..., X m ) m t = (t 1,..., t m ) M X (t) = E(e t X ) 86

87 [ 1.32] m X E(X) = µ V (X) = Σ Y = AX + b A k m b k [ 1.33] X X N m (µ, Σ) Y [ 1.10] Y n Y g lim E(g(Y n)) = E(g(Y )) (1.9) n Y n Y Y n d Y n F d Y Y F [ ] Y n f n (y) Y f(y) Y n d Y f n f g(x) = I(x A) (1.9) lim n f n (A) = f(a) A = (, y) lim n F n (y) = F (y) F (y) y 87

88 Y n d Y [ 1.2] X n E(X 1 ) = 0 V (X 1 ) = 1 n n Y n = 1 n n X i 1.1 Y n = n X n n i=1 Y n d N(0, 1) (1.10) [ ] E(X 1 ) = µ V (X 1 ) = σ 2 Z n = (X n µ)/σ Z n n Z n d N(0, 1) n n Zn N(0, 1) (asymptotic normality) n Z a a n N(0, 1) Zn N(0, 1/n) a X n N(µ, σ 2 /n) E( X n ) = µ V ( X n ) = σ 2 /n X n [ ] E(X1 k ), k = 1, 2,... X 1 ϕ X1 (t) = E(e itx 1 ) = 88

89 k=0 (it) k k! E(X k 1 ) = 1 t2 2 it3 6 E(X3 1 ) + t4 24 E(X4 1 ) + ϕ Yn (t) = [ ϕ X1 (t/ n) ] n = [ 1 t2 2n lim n ϕ Yn (t) = e t2 /2 it3 6n 3/2 E(X3 1 ) + t4 24n 2 E(X4 1 ) + [ 1.16] > ## myclt: n b > ## : n= b= > ## : b > myclt <- function(n,b) { + X <- matrix(runif(n*b),n,b) # n*b U(0,1) + a <- rep(1/n,n) # 1/n n + y <- a %*% X # b + x0 <- seq(min(y),max(y),length=1000) # y truehist(y) # + lines(x0,dnorm(x0,mean=0.5,sd=sqrt(1/120)),col=2) # + } > myclt(5,10000) > myclt(10,10000) 19 [ 1.17] > ## myclt2: n ( p 1 1-p 0) b > ## : n= p=+1 b= > ## : b 89 ] n = ] n [1 t2 (1 + o(1)) 2n

90 y y 19 n n = 5 n = 10 90

91 > myclt2 <- function(n,p,b) { + X <- matrix(as.numeric(runif(n*b)<p),n,b) # n*b a <- rep(1,n) # 1 n + y <- a %*% X # b + x0 <- seq(from=min(y)-0.5,to=max(y)+0.5) # truehist(y,breaks=x0) # + x0 <- seq(from=min(y),to=max(y)) # + points(x0,dbinom(x0,size=n, prob=p)) # + points(x0,dpois(x0,lambda=n*p),col=3) # + x0 <- seq(from=min(y)-0.5,to=max(y)+0.5,length=10000) # lines(x0,dnorm(x0,mean=n*p,sd=sqrt(n*p*(1-p))),col=2) # + } > myclt2(10,0.05,10000) > myclt2(10,0.5,10000) 20 [ 1.34] n p X (i) λ > 0 p = λ/n n X λ 1.31 M X (t) M Y (t) (ii) 0 < p < 1 n (X np)/ n 0 p(1 p) [ 1.35] X n m E(X) = 0 V (X) = Σ = (σ ij ) 91

92 y y 20 n n = 10, p = 0.05 n = 10, p =

93 X n = n i=1 X i/n n X n d Nm (0, Σ) (1.11) N m E(X k 1 1 Xk m m ), k 1,..., k m 0 ϕ n X(t) = E(e i nt X ) 93

94 2 EM Fisher 2.1 m X q(x) X q(x) q(x) X q n m X 1,..., X n q(x) (i.i.d.=independently identically distributed) X 1,..., X n q(x) (i.i.d.) n X = (x 1,..., x n ) [ 2.1] q(x) p θ Θ R p f(x; θ) 94

95 [ 2.2] f(x; θ) θ 0 Θ q = f( ; θ 0 ) x q(x) = f(x; θ 0 ) [ ] [ 2.1] X N(µ, σ 2 ) θ = (µ, σ 2 ), p = 2 m X N m (µ, Σ) θ = (µ, Σ), p = m(m + 3)/2 [ 2.2] k i X f i (x; θ i ) k i p(i) = π i > 0 k i=1 π i = 1 X i f(x, i; θ) = f i (x; θ i )π i θ = (π 1,..., π k 1, θ 1,..., θ k ) f(x, i; θ) = f(x i; θ)p(i; θ) i X N(µ i, σi 2) θ i = (µ i, σi 2) f i(x; θ i ) = 1 exp( (x µ i) 2 ) 2πσ 2 2σ 2 i i [ 2.3] 2.2 X i X 95

96 X f(x; θ) = k f i (x; θ i )π i i=1 [ 2.3] 2.3 (mixture distribution) (normal mixture model) [ ] i x i x i x f(x i) i p(i) i, j f(x i, j) i j x ( X) 96

97 [ 2.4] 2.3 n 1 = 100, k = 3, π 1 = 0.5, π 2 = 0.3, π 3 = 0.2, µ 1 = 0, µ 2 = 4, µ 3 = 3, σ 1 = 1, σ 2 = 2, σ 3 = 1. > ### > ## : n1 > ## : k > ## : pr[1],...,pr[k-1] > ## : mu[1],...,mu[k] > ## : ss[1],...,ss[k] > ## xx1 > n1 <- 100; k <- 3 > pr <- c(0.5,0.3); mu <- c(0,4,-3); ss <- c(1,2,1)^2 > pr[k] <- 1 - sum(pr) > ## > ## x0= x (k,pr,mu,ss)= col= > drawnormmix <- function(x0,k,pr,mu,ss,col) { + ## i f(x0[t] i)*p(i),t=1,2,.. + fi <- matrix(0,length(x0),k) # fi(xt), i=1,...,k, t=1,...,n + for(i in 1:k) fi[,i] <- pr[i]*dnorm(x0,mean=mu[i],sd=sqrt(ss[i])) + ## f(x0[t]),t=1,2,... + f <- apply(fi,1,sum) # f(xt), t=1,...,n + ## + if(col!= FALSE) { + matlines(x0,fi,col=col,lty=1) # + lines(x0,f,col=col,lty=2,lwd=2) # 97

98 + } + invisible(list(f=f,fi=fi)) # + } > ## > ii1 <- sample(k,n1,replace=true,prob=pr) # > ii1[1:30] # [1] > zz1 <- rnorm(n1) # N(0,1) > xx1 <- zz1*sqrt(ss)[ii1] + mu[ii1] # > x0 <- seq(min(xx1),max(xx1),length=400) # 400 > hist(xx1,nclas=10,prob=t) # > rug(xx1) # > f0 <- drawnormmix(x0,k,pr,mu,ss,"darkgreen")$f # (f0 ) (x, i) i = 1,..., k > ## pr1,mu1,ss1 > pr1 <- mu1 <- ss1 <- rep(0,k) # > for(i in 1:k) { + pr1[i] <- sum(ii1==i)/length(ii1) # i + x <- xx1[ii1==i] # i + mu1[i] <- mean(x) # + ss1[i] <- mean((x-mu1[i])^2) # + } > pr1 # pr [1] > mu1 # mu [1] > sqrt(ss1) # sqrt(ss) 98

99 [1] > drawnormmix(x0,k,pr1,mu1,ss1,"red") # n 1 X 1 = {(x t, i t ), t = 1,..., n 1 } log L(θ X 1 ) = n 1 t=1 log (f it (x t ; θ it )π it ) 2.2 [ 2.5] 2.3 x i p(i x) π i (x; θ) = f i(x; θ i )π i f(x; θ) i î(x; θ) = arg max i=1,...,k π i(x; θ) x î(x; θ) θ ˆθ 99

100 x i i.i.d. [ 2.6] 2.4 n 2 = 300 i x î(x; ˆθ) 2.4 (x t, i t ) i > ## : n2 > ## xx2 > n2 <- 300 > ii2 <- sample(k,n2,replace=true,prob=pr) # > ii2[1:30] # [1] > zz2 <- rnorm(n2) # N(0,1) > xx2 <- zz2*sqrt(ss)[ii2] + mu[ii2] # > hist(xx2,nclass=15,prob=t) # > rug(xx2) # > drawnormmix(x0,k,pr,mu,ss,"darkgreen") # > ## i ii2 xx2 ( xx1 ii1 > a <- drawnormmix(xx2,k,pr1,mu1,ss1,false) # n1 > round(t(a$fi[1:10,]),3) # f(xx2[t] i)*p(i) 10 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] [2,] [3,] > round(a$f[1:10],3) # f(xx2[t])

101 [1] > pr2x1 <- a$fi / a$f # p(i x) > round(t(pr2x1[1:10,]),3) # 10 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] [2,] [3,] > ## ( ) : ii2x1 pr2x1 a$fi > ii2x1 <- apply(pr2x1,1,function(p) order(-p)[1]) > ii2x1[1:30] # [1] > sum(ii2x1 == ii2)/length(ii2) # [1] [ 2.7] k = 2 i = 0 (ham) i = 1 (spam) m j = 1,..., m x[j] = 1 x[j] = 0 i x p(x i) f i (x i; θ i ) = p(x[1] i)p(x[2] i) p(x[m] i) i x[j] = 1 p(x[j] i) = θ i [j] x[j] = 0 p(x[j] i) = 1 θ i [j] θ = (π 1, θ 0 [1],..., θ 0 [m], θ 1 [1],..., θ 1 [m]) 101

102 Histogram of xx1 Histogram of xx2 Density Density xx xx2 21 (x t, i t ), t = 1,..., n 1 102

103 Naive Bayes Spambase UCI Repository of machine learning databases ~mlearn/mlrepository.html > ### > load("spam1.rda") # dat1.train,spam.train,dat1.test,spam.test > dim(dat1.train) # n=3601, m=54 [1] > t(dat1.train[1:20,1:10]) # Iword_freq_make Iword_freq_address Iword_freq_all Iword_freq_3d Iword_freq_our Iword_freq_over Iword_freq_remove Iword_freq_internet Iword_freq_order Iword_freq_mail > spam.train[1:20] # spam=1, ham=0 [1] > dim(dat1.test) # [1] > ### theta > ## x : 0,1 ( dat) > ## y : 0,1 ( spam) 103

104 > mymle <- function(x,y) { + py1 <- mean(y) # p(y=1) + px0 <- apply(x[y==0,],2,mean) # p(x[j]=1 y=0) + px1 <- apply(x[y==1,],2,mean) # p(x[j]=1 y=1) + list(py1=py1,px0=px0,px1=px1) # theta + } > ### p(y=1 x) > ## th : list(py1,px0,px1) > ## x : 0,1 ( dat) > mypp <- function(th,x) { + x <- as.matrix(x) # x + x0 <- x1 <- x + for(j in seq(ncol(x))) { # j=1,...,m + a <- x[,j] + 1 # t=1,...,n j =2, =1 + x0[,j] <- c(1-th$px0[j],th$px0[j])[a] # p(x[j] y=0) + x1[,j] <- c(1-th$px1[j],th$px1[j])[a] # p(x[j] y=1) + } + p0 <- apply(x0,1,prod) # p(x y=0) + p1 <- apply(x1,1,prod) # p(x y=1) + th$py1*p1/(th$py1*p1 + (1-th$py1)*p0) # p(y=1 x) + } > ### > myppplot <- function(spam,pp,pth) { + ## + def.par <- par(no.readonly = TRUE); on.exit(par(def.par)) + layout(matrix(1:2,2,1)) + ## spam/ham + sp <- spam==1 # spam=true, ham=false + p0 <- mean(pp[!sp]>pth) # ham spam 104

105 + p1 <- mean(pp[sp]>pth) # spam spam + ## ham + hist(pp[!sp],col="blue",nclass=50,prob=t,main="ham mails", + sub=paste("p(say spam ham mail)=",round(p0,5))) + abline(v=pth,col="green") + ## spam + hist(pp[sp],col="red",nclass=50,prob=t,main="spam mails", + sub=paste("p(say spam spam mail)=",round(p1,5))) + abline(v=pth,col="green") + ## + ret <- c(pth,p0,p1) + names(ret) <- c("pth","p0","p1") + ret + } > ### > ## > th <- mymle(dat1.train,spam.train) # th > names(th) # [1] "py1" "px0" "px1" > th$py1 # p(y=1) [1] > round(th$px0[1:10],3) # p(x[j]=1 y=0) Iword_freq_make Iword_freq_address Iword_freq_all Iword_freq_3d Iword_freq_our Iword_freq_over Iword_freq_remove Iword_freq_internet Iword_freq_order Iword_freq_mail > round(th$px1[1:10],3) # p(x[j]=1 y=1) 105

106 Iword_freq_make Iword_freq_address Iword_freq_all Iword_freq_3d Iword_freq_our Iword_freq_over Iword_freq_remove Iword_freq_internet Iword_freq_order Iword_freq_mail > ## > pp.train <- mypp(th,dat1.train) # > round(pp.train[1:20],3) # p(y=1 x) > myppplot(spam.train,pp.train,0.5) # 0.5 pth p0 p > ## > pp.test <- mypp(th,dat1.test) # > round(pp.test[1:20],3) # p(y=1 x) > spam.test[1:20] # [1] > myppplot(spam.test,pp.test,0.5) # 0.5 pth p0 p > ## ham spam

107 > pth <- quantile(pp.train[spam.train==0],p=0.99) > pth # p(y=1 x)>pth spam 99% > myppplot(spam.train,pp.train,pth) # pth pth p0 p > myppplot(spam.test,pp.test,pth) # pth pth p0 p [ 2.1] 2.5 x m f i (x; θ i ) x N m (µ i, Σ i ) θ i = (µ i, Σ i ) i j π i f i (x; θ i ) > π j f i (x; θ j ) i j x S i (x) = log(π i f i (x; θ i )) [ 2.2] [ 2.3] 2.1 i Σ i = Σ x 2.5 θ = (π 1,..., π k 1, θ 1,..., θ k ) i x 107

108 ham mails ham mails Density pp[!sp] P(say spam ham mail)= spam mails Density Density pp[!sp] P(say spam ham mail)= spam mails Density pp[sp] P(say spam spam mail)= pp[sp] P(say spam spam mail)=

109 R 1 R 2 R k x R i i R = {R 1,..., R k } P (R) = k i=1 R i π i f i (x; θ i ) dx R = {R 1,..., R k } x R i π i f i (x; θ i ) π j f j (x; θ j ), j = 1,..., k P (R ) P (R) 2.3 X = (x 1,..., x n ) f(x; θ) θ ˆθ ˆθ X f(x; θ) ˆθ [ 2.4] (X 1,..., X n ) θ L(θ X ) = f(x 1 ; θ) f(x n ; θ) 109

110 [ 2.5] θ π(θ) π(θ X ) L(θ X )π(θ) θ Θ π(θ X ) MAP Maximum A Posteriori estimator) ˆθ MAP = arg max θ Θ π(θ X ) [ 2.4] X 1,..., X n µ N(µ, 1) (i.i.d.), µ N(0, τ 2 ) µ (i) L(µ X ) exp( n 2 (µ x)2 ) x = n t=1 x t/n (ii) µ MAP ˆµ MAP = x x 1+nτ 2 [ 2.6] estimator) θ Θ L(θ X ) Maximum Likelihood ˆθ ML = arg max θ Θ L(θ X ) MAP π(θ) = 110

111 [ 2.5] X 1,..., X n N(µ, σ 2 ) (i.i.d.) θ = (µ, σ 2 ) ML ˆµ ML = x, ˆσ 2 ML = n t=1 (x t x) 2 /n [ 2.6] X 1,..., X n N m (µ, Σ) (i.i.d.) θ = (θ, Σ) ˆµ ML = x, ˆΣ ML = n t=1 (x t x)(x t x) /n [ 2.7] k 2 X {1, 2,..., k} θ i = P (X = i), i = 1,..., k 1 θ = (θ 1,..., θ k 1 ) P (X = k) = 1 k 1 i=1 θ i x 1,..., x n i z i ˆθ i = z i /n (Z 1,..., Z k ) [ 2.8] 2.2 (x t, i t ), t = 1,..., n X i θ θ i i = 1,..., k π i i ˆθ i = arg max θ i n I(i t = i) log f i (x t ; θ i ), t=1 ˆπ i = z i n z i = n t=1 I(i t = i) ˆµ i = n t=1 I(i t = i)x t /z i ˆσ 2 i = n t=1 I(i t = i)(x t ˆµ i ) 2 /z i

112 X 2.5 [ 2.8] i X 2.3 f(x; θ) 2.6 n 2 x t (i t θ X 2 = {x t, t = n 1 + 1,..., n 1 + n 2 } log L(θ X 2 ) = n 1 +n 2 log f(x t ; θ) = n 1 +n 2 ( k ) log f i (x t ; θ i )π i t=n 1 +1 t=n 1 +1 i=1 > ## > ## : xx2 > ## : k > ## pr2,mu2,ss2 > ## > mylik2 <- function(theta) { + pr <- theta[1:(k-1)]; mu <- theta[k:(2*k-1)]; ss <- theta[(2*k):(3*k-1)] + pr[k] <- 1 - sum(pr) + f <- drawnormmix(xx2,k,pr,mu,ss,false)$f # xx2 + -sum(log(f)) # *(-1) + } > ## (optim > th0 <- c(1/3,1/3,0,2,-2,1,1,1) # > opt2 <- optim(th0,mylik2,method="bfgs",control=list(trace=1,reltol=1e-14),hessian=true) initial value iter 10 value

113 iter 20 value iter 30 value iter 40 value final value converged > th2 <- opt2$par > pr2 <- th2[1:(k-1)]; mu2 <- th2[k:(2*k-1)]; ss2 <- th2[(2*k):(3*k-1)] > pr2[k] <- 1 - sum(pr2) > a <- order(-pr2) # pr2 > pr2 <- pr2[a]; mu2 <- mu2[a]; ss2 <- ss2[a] # > pr2 # pr [1] > mu2 # mu [1] > sqrt(ss2) # sqrt(ss) [1] > hist(xx2,nclass=15,prob=t) # > rug(xx2) # > drawnormmix(x0,k,pr2,mu2,ss2,"blue") # > lines(x0,f0,col="darkgreen",lty=2,lwd=2) ˆθ = (ˆπ 1,..., ˆπ k 1, ˆµ 1,..., ˆµ k, ˆσ 2 1,..., σ 2 k ) i > ## i ii2 xx2 > a <- drawnormmix(xx2,k,pr2,mu2,ss2,false) # n2 x > pr2x2 <- a$fi / a$f > round(t(pr2x2[1:10,]),3) #

114 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] [2,] [3,] > ## ( ) MAP : ii2x2 > ii2x2 <- apply(pr2x2,1,function(p) order(-p)[1]) > ii2x2[1:30] # [1] > sum(ii2x2 == ii2)/length(ii2) # [1] k x 1,..., x n x t x t i [ 2.9] n 1 X 1 = {(x t, i t ), t = 1,..., n 1 } n 2 X 2 = {x t, t = n 1 + 1,..., n 1 + n 2 } log L(θ X 1, X 2 ) = log L(θ X 1 ) + log L(θ X 2 ) = n 1 t=1 log (f it (x t ; θ it )π it ) + n 1 +n 2 t=n 1 +1 ( k ) log f i (x t ; θ i )π i i=1 114

115 [ 2.9] 2.9 X 1 X 2 ˆθ > ## > ## : xx1, ii1, xx2 > ## : k > ## pr3,mu3,ss3 > ## > mylik3 <- function(theta) { + pr <- theta[1:(k-1)]; mu <- theta[k:(2*k-1)]; ss <- theta[(2*k):(3*k-1)] + pr[k] <- 1 - sum(pr) + ## xx1,ii1 + fi <- drawnormmix(xx1,k,pr,mu,ss,false)$fi + f1 <- fi[seq(n1) + (ii1-1)*n1] + ## xx2 + f2 <- drawnormmix(xx2,k,pr,mu,ss,false)$f # xx2 + ## + -sum(log(c(f1,f2))) # *(-1) + } > ## > th0 <- c(1/3,1/3,0,2,-2,1,1,1) # > opt3 <- optim(th0,mylik3,method="bfgs",control=list(trace=1,reltol=1e-14)) # initial value

116 iter 10 value iter 20 value iter 30 value iter 30 value iter 30 value final value converged > th3 <- opt3$par > pr3 <- th3[1:(k-1)]; mu3 <- th3[k:(2*k-1)]; ss3 <- th3[(2*k):(3*k-1)] > pr3[k] <- 1 - sum(pr3) > pr3 # pr [1] > mu3 # mu [1] > sqrt(ss3) # sqrt(ss) [1] > hist(xx2,nclass=15,prob=t) # > rug(xx2) # > drawnormmix(x0,k,pr3,mu3,ss3,"orange") # > lines(x0,f0,col="darkgreen",lty=2,lwd=2) > ## i ii2 xx2 > a <- drawnormmix(xx2,k,pr3,mu3,ss3,false) # n1 (x,i) n2 x > pr2x3 <- a$fi / a$f > round(t(pr2x3[1:10,]),3) # 10 [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] [2,] [3,] > ## ( ) MAP : ii2x3 116

117 > ii2x3 <- apply(pr2x3,1,function(p) order(-p)[1]) > ii2x3[1:30] # [1] > sum(ii2x3 == ii2)/length(ii2) # [1] > ### > sum((c(pr1,mu1,ss1)-c(pr,mu,ss))^2) # n1 (x,i) [1] > sum((c(pr2,mu2,ss2)-c(pr,mu,ss))^2) # n2 x [1] > sum((c(pr3,mu3,ss3)-c(pr,mu,ss))^2) # n1 (x,i) n2 x [1] [ ] R optim EM 2.4 EM [ 2.10] 2.9 r = 0, 1, 2,... r = 0 θ (r) θ (r+1) 117

118 Histogram of xx2 Histogram of xx2 Density Density xx xx2 23 n 2 x t n 1 (x t, i t ) n 2 x t 118

119 n = n 1 + n 2 ˆ E (Expectation step) θ (r) t = 1,..., n 1 + n 2 i p(i t x t ) = π i (x t ; θ (r) ) π i (x t ; θ (r) ) = I(i = i t ), t = 1,..., n 1 X 1 i t π i (x t ; θ (r) ) = f i(x t ; θ (r) i )π (r) i, t = n f(x t ; θ (r) 1 + 1,..., n ) ˆ M (Maximization step) π 1,..., π k ( i=1 π i = 1) π (r+1) i = 1 n n π i (x t ; θ (r) ) t=1 x t w t = π i (x t ; θ (r) ) n t =1 π i(x t ; θ (r) ), t = 1,..., n 119

120 i = 1,..., k µ (r+1) i = σ 2 (r+1) i = n w t x t, t=1 n t=1 t = 1,..., n w t (x t µ (r+1) i ) 2, t = 1,..., n E M optim log L(θ (r) X 1, X 2 ) > ## EM > ## : xx1, ii1, xx2 > ## : k > ## pr4,mu4,ss4 > pr4 <- c(1/3,1/3,1/3); mu4 <- c(0,2,-2); ss4 <- c(1,1,1) # > nr <- 30 # > mystat <- function(pr,mu,ss) { # + lik <- mylik3(c(pr[-k],mu,ss)) # + cat(format(lik,digits=10),round(c(pr,mu,ss),3),"\n") # + c(lik,pr,mu,ss) + } > stat3 <- mystat(pr3,mu3,ss3) # optim mylik > stat4 <- matrix(0,1+nr,length(stat3)) # > stat4[1,] <- mystat(pr4,mu4,ss4) # 120

121 > xx <- c(xx1,xx2) # > pp1 <- matrix(0,n1,k); pp1[seq(n1)+(ii1-1)*n1] <- 1 # xx1 > t(pp1[1:10,]) # [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] [2,] [3,] > ## EM > for(r in 1:nr) { # : break + a <- drawnormmix(xx2,k,pr4,mu4,ss4,false); pp2 <- a$fi/a$f # xx2 + pp <- rbind(pp1,pp2) # xx1 xx2 + pr4 <- apply(pp,2,sum)/(n1+n2) # pr + # wt <- pp/rep(apply(pp,2,sum),rep(n1+n2,k)) # + wt <- sweep(pp,2,apply(pp,2,sum),"/") # sweep + mu4 <- apply(xx*wt,2,sum) # mu + # ss4 <- apply((xx-rep(mu4,rep(n1+n2,k)))^2*wt,2,sum) # ss + ss4 <- apply(sweep(matrix(xx,n1+n2,k),2,mu4,"-")^2*wt,2,sum) # + cat(r,": ") + stat4[r+1,] <- mystat(pr4,mu4,ss4) # + } 1 : : : : : : : :

122 9 : : : : : : : : : : : : : : : : : : : : : : > pr4 # pr [1] > mu4 # mu [1] > sqrt(ss4) # sqrt(ss) [1]

123 > plot(0:nr,stat4[,1],type="b",xlab="iteration",ylab="lik") # > abline(h=stat3[1],lty=2,col="pink") # optim > matplot(0:nr,stat4[,-1],type="b",xlab="iteration",ylab="parameters") # > abline(h=stat3[-1],lty=2,col="pink") # optim [ 2.7] X Y 2.9 X = {X 1, X 2 } = {x 1,..., x n1 +n 2, i 1,..., n n1 } Y = {i n1 +1,..., i n1 +n 2 } (X, Y) L(θ X, Y) X L(θ X ) = L(θ X, Y) dy Y 2.9 L(θ X ) EM r = 0 θ (r) θ (r+1) ˆ E (Expectation step) θ (r) X Y f(y X ; θ (r) ) f(y X ; θ (r) ) = L(θ(r) X, Y) L(θ (r) X ) Y log(l(θ X, Y)) Q(θ, θ (r) ) = log(l(θ X, Y))f(Y X ; θ (r) ) dy 123

124 lik parameters iteration iteration 24 log L(θ (r) X 1, X 2 ), r = 0,..., optim θ (r) = (π (r) 1, π(r) 2, µ(r) 1, µ(r) 2, µ(r) (r) 3, σ2 1, σ 2 (r) 2, σ 2 (r) 3 ) 124

125 ˆ M (Maximization step) Q(θ, θ (r) ) θ θ (r+1) = arg max θ L(θ (r+1) X ) L(θ (r) X ) Q(θ, θ (r) ) [ ] L(θ X ) = L(θ X, Y)/f(Y X ; θ) log L(θ X ) = log L(θ X, Y) log f(y X ; θ) f(y X ; θ (r) ) log L(θ X ) = Q(θ, θ (r) ) H(θ, θ (r) ) H(θ, θ (r) ) = log(f(y X ; θ))f(y X ; θ (r) ) dy θ H(θ, θ (r) ) H(θ (r), θ (r) ) Q(θ, θ (r) ) Q(θ (r), θ (r) ) θ log L(θ X ) log L(θ (r) X ) [ ] (i) EM log L(θ X ) EM (ii) µ 1 125

126 x t log L(θ X ) [ 2.10] f(x) g(x) < x < f(x) > 0, g(x) > 0 log(g(x))f(x) dx log(f(x))f(x) dx [ 2.11] 2.7 EM X = {x 1, x 2,..., x n } X 1, X 2,..., X n g(x) (i.i.d.) g(x) f(x; θ), θ Θ m θ = (θ 1,..., θ m ) f(x; θ) g(x) f(x; θ) X 1,..., X n X E θ ( ) V θ ( ) C θ ( ) L(θ X ) = f(x 1,..., x n ; θ) = f(x 1 ; θ)f(x 2 ; θ) f(x n ; θ) 126

127 ˆθ ML ˆθ ML (X ) ˆθ ML (x 1,..., x n ) ˆθ [ 2.8] f(x; θ) ˆθ (unbiased) E θ (ˆθ(X 1,..., X n )) = θ [ 2.1] ˆθ V θ (ˆθ(X 1,..., X n )) 1 n G(θ) 1 (2.1) G(θ) A, B A B A B (non-negative definite) m m G(θ) G(θ) Fisher G ij (θ) = E θ { log f(x; θ) θ i log f(x; θ) θ j } 127

128 [ ] G(θ) X t n ng(θ) Fisher 2.14 ( ) ng(θ) = E θ 2 log L θ θ [ ] E θ (ˆθ(X 1,..., X n )) = θ ˆθ i (x 1,..., x n )f(x 1,..., x n ; θ) dx 1 dx n = θ i, i = 1,..., m θ j ˆθ i (x 1,..., x n ) log f(x 1,..., x n ; θ) f(x 1,..., x n ; θ) dx 1 dx n = θ i θ j θ j m S(x 1,..., x n ; θ) X = (X 1,..., X n ) S j (x 1,..., x n ; θ) = log f(x 1,..., x n ; θ) θ j, j = 1,..., m E θ {ˆθ(X )S(X ; θ) } = I m (2.2) 128

129 ˆθ i = θ i = 1 E θ {S(X ; θ)} = 0 (2.3) (2.2) (2.3) } C θ {ˆθ(X ), S(X ; θ) = I m (2.4) {[ ]} ˆθ(X ) V θ S(X ; θ) = [ ] Vθ {ˆθ(X )} I m V θ {S(X ; θ)} I m A = V θ {ˆθ(X )}, B = V θ {S(X ; θ)} m a, b [ a b ] [ A Im B I m b = B 1 a ] [ ] a b = a Aa + 2a b + b Bb 0 a Aa 2a B 1 a + a B 1 a = a (A B 1 )a 0 129

130 A B 1 i.i.d. B = ng (2.1) [ 2.12] { } log f(x; θ) E θ θ i = 0, i = 1,..., m [ 2.13] [ 2.2] n θ ML θ 1 n G(θ) 1 n(ˆθ(x 1,..., X n ) θ) d N(0, G(θ) 1 ) (2.5) [ ] [ ] ˆθ ML Θ log f(x 1,..., x n ; θ) = 0 θ ˆθ ML 130