#9.2.2011/ Binomiuskottavuus # Kannatusmittaus, 10.1. ########################## # Binomijakauma P(X=10) par(mfrow=c(1,2)) #par(mfrow=c(1,1)) theta<- seq(.0,.3,len=100) like <- dbinom(10,100,theta) nlike<- like/max(like) plot(theta,like,type="l",xlab=expression(theta), ylab='Uskottavuus') #lines(theta,nlike,lty=2) plot(theta,nlike,type="l",xlab=expression(theta), ylab='Uskottavuus/normeerattu') title(expression('Uskottavuusfunktio ja normeerattu uf, n=100')) ############################### # Esimerkki 10.1/ ########################## # (a) Binomijakauma P(X=5) & P(X<11) par(mfrow=c(1,1)) theta<- seq(.0,.3,len=100) like <- pbinom(10,100,theta) like<- like/max(like) plot(theta,like,type="l",xlab=expression(theta), ylab='Uskottavuus') like<-dbinom(5,100,theta) like<- like/max(like) lines(theta,like,lty=2) text(.12,.7,'x<11',cex=.8) text(.050,.4,'x=5',cex=.8) title(expression('Binomiuskottavuus, n=100')) ############################### # (b) P(X=5) & P(X<11) & P(X>15 ############################### par(mfrow=c(1,1)) theta<- seq(.0,.3,len=100) like <- pbinom(10,100,theta) like<- like/max(like) plot(theta,like,type="l",xlab=expression(theta), ylab='Uskottavuus') like<-dbinom(5,100,theta) like<- like/max(like) lines(theta,like,lty=2) text(.12,.7,'x<11',cex=.8) text(.050,.4,'x=5',cex=.8) like <- 1-pbinom(14,100,theta) # P(X>15) like<- like/max(like) lines(theta,like,lty=3) text(.20,.4,'x>15',cex=.8) title(expression('Binomiuskottavuus, n=100')) ########################################## # Esimerkki 10.2. # Myyrien lukumäärän estimointi ########################################## N1<- 25 # Merkittyjä n<- 60 # Otos merkitsemisen jälkeen x<- 5 # Otoksessa havaitaan x merkittyä NN<-seq(150,1100,by=25) #NN<- c(seq(150,200),seq(200,725,by=25), # seq(725,775),seq(775,1100,by=25)) loglik<-NULL for (N in NN){ ll <- lgamma(N-N1+1) - lgamma(N-N1-n+x+1) - lgamma(N+1) + lgamma(N-n+1) loglik <- c(loglik,ll) } like<- exp(loglik-max(loglik)) plot(NN,like,type='l',xlab='N',ylab='Uskottavuus') title(expression('Esimerkki 10.2: Myyrien lukumäärä')) #............................................. NN[which.max(loglik)] # Myyrien lukumäärän SUE ########################################## ########################################## # Esimerkki 10.2.3 # ....................................... Negatiivinen binomijakauma: # # Geneetikot tutkivat erään harvinaisen genotyypin esiintyvyyttä (prevalenssi) # populaatiossa. Tutkitaan niin monta yksilöä, kunnes löydetään kyseinen # genotyyppi. Genotyyppi löytyi, kun analysoitiin 53:s yksilö. Kun eri yksilöt # ovat toisistaan riippumattomia, niin prevalenssitodennäköisyyden uskottavuus # saadaan geometrisesta jakaumasta. ################################## theta<- seq(.0001,.1,len=100) like <- dgeom(52,theta) #like <- dbinom(1,53,theta) like1<- like/max(like) # uskottavuus/geometrinen jakauma ################################### # Toisessa kokeessa Geneetikot tutkivat niin monta, että löysivät 5 etsimäänsä # genotyyppiä. Tähän tarvittiin 552 yksilöä. ################################### like<- dnbinom(552-5,5,theta) # uskottavuus/negatiivinen binomijakaumajakauma #like<- dbinom(4,551,theta) like2<- like/max(like) plot(c(0,.1),range(c(like1,like2)),type='n', xlab=expression(theta), ylab='Likelihood') lines(theta,like1) lines(theta,like2,lty='dotted',lwd=1.5) title(expression('(c) Prevalenssitodennäköisyys')) ################################### ################################### # Esimerkki 10.2.4 # ....................................... Esimerkki Normaalijakaumasta theta <- seq(-1,6,len=40) like <- pnorm(4,theta) - pnorm(0.9,theta) # Phi(4-theta)-Phi(0.9-theta) like <- like/max(like) plot(theta,like,type='n', xlab=expression(theta), ylab='Uskottavuus') lines(theta,like) # Yhtenäinen viiva like <- dnorm(2.45,theta) like<- like/max(like) lines(theta,like,lty='dashed') # Pätkäviiva n<- 5 xn<- 3.5 llike <- log(dnorm(xn,mean=theta)) + (n-1)*log(pnorm(xn,mean=theta)) like <- exp(llike-max(llike)) lines(theta,like,lty='dotted',lwd=1.5) # Pisteviiva title(expression('Normaalijakauman sijaintiparametrin uskottavuus')) #################################### # Esimerkki/ Otos Bernoullin jakaumasta # Torakan myrkyn vaikutus torakan muniin # 1 = munasta sytyy torakka # 0 = munasta ei synny torakkaa # Munien lukumäärä noudattaa Binomijakaumaa munat<-scan() 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 mean(munat) # SUE theta<-seq(0.1,0.9,0.001) negloglike<- function(theta){-(sum(munat)*log(theta)+sum(1-munat)*log(1-theta))} nlm(negloglike,0.2) # nml on minimointiproseduuri, joten SUE saadaan # minimoimalla negatiivinen log-uskottavuus tulos<-nlm(negloglike,0.2) par(pty = "s") # par(mfrow=c(1,1),"s") theta <- seq(0.1, 0.9, 0.001) plot(theta, - negloglike(theta), type = "n", ylab = "l(theta)") abline(v = mean(munat), col = 13, lwd = 3) lines(theta, - negloglike(theta), col = 6, lwd = 3) # curve(-negloglike(x),0,1) # Maksimoidaan log-uskottavuus funktiolla optimize loglike <- function(theta){(sum(munat)*log(theta) + sum(1-munat) * log(1-theta))} optimize(f=loglike,interval=c(0,1),maximum=TRUE) ############################ # 10.8 Pistesuureen jakauma # Esimerkki: Viiden alkion otos Bernoullin jakaumasta. like<- function(theta,x){(theta^x)*((1-theta)^(5-x))} # uskottavuusfunktio #par(pty = "s") # par(mfrow=c(1,1),"s") theta <- seq(0, 1, 0.001) plot(theta, like(theta,0), type = "n", ylab = "L(theta)") lines(theta, like(theta,0), lwd = 2) lines(theta, like(theta,1), lwd = 2) lines(theta, like(theta,2), lwd = 2) lines(theta, like(theta,3), lwd = 2) lines(theta, like(theta,4), lwd = 2) lines(theta, like(theta,5), lwd = 3) # Generoidaan uskottavuusfunktioita plot(theta, like(theta,1), type = "n", ylab = "L(theta)") lines(theta, like(theta,rbinom(1,5,0.5)), lwd = 2) ###################### # Logaritmoitu uf loglike<-function(theta,x){x*log(theta)+(5-x)*log(1-theta)} # logaritmoitu uf theta <- seq(0, 1, 0.001) plot(theta, loglike(theta,0), type = "n", ylab = "l(theta)") lines(theta, loglike(theta,0), lwd = 2) lines(theta, loglike(theta,1), lwd = 2) lines(theta, loglike(theta,2), lwd = 2) lines(theta, loglike(theta,3), lwd = 2) lines(theta, loglike(theta,4), lwd = 2) lines(theta, loglike(theta,5), lwd = 3) # Generoidaan logaritmoituja uskottavuusfunktiota plot(theta, loglike(theta,0), type = "n", ylab = "l(theta)") lines(theta, loglike(theta,rbinom(1,5,0.5)), lwd = 2)