## Intervaly spolehlivosti # rozsah xs<-seq(0,20,.1) mu<-10 so<-4 # cheme vykreslit 2x2 grafy (ala matlab subplot(2,2,x)) a do 'opar' ulozime současné nastavení opar<-par(mfrow=c(2,2)) ## Příklad: glykémie pacientů s cukrovkou. # 1) Vygeneruje náhodný vzorek z daného teoretického rozdělení a vykreslíme jej spolu s teoretickou hustotou psti. # velikost vzorku n.vyber<-100 # náhodný vzorek hodnot glykemie set.seed(123) # kvůli reprodukovatelnosti x<-rnorm(n.vyber,mu,so) # histogram h<-hist(x,breaks=50,plot=FALSE) # toto je jen pro výpočet max. hodnoty hustoty kvůli nastavení většího rozsahu mezí osy Y hist(x,col='gray',probability=TRUE,breaks=50,ylim=c(0,2*max(h$density)),xlab='glykémie [mmol/l]',ylab='hustota psti',main='Hodnoty glykémie') # proložíme hustotu pravděpodobnosti lines(xs,dnorm(xs,mu,so),ty='l',col='red',lwd=2) legend('topright',c('glykémie','rozdělení glykémií'),col=c('gray','red'),lwd=c(10,2)) # 2) Přidáme průmery náhodných vzorků a teoretické rozdělení průměrů. # vykreslíme nejdřív pro srovnání původní obrázek hist(x,col='gray',probability=TRUE,breaks=50,ylim=c(0,2*max(h$density)),xlab='glykémie [mmol/l]',ylab='hustota psti',main='Průměry glykémie') lines(xs,dnorm(xs,mu,so),ty='l',col='red',lwd=2) legend('topright',c('glykémie','rozdělení glykémií','průměry','rozdělení průměrů'),col=c('gray','red',scales::alpha('blue',.5),'blue'),lwd=c(10,2,10,2)) # počet průměrovaných vzorků n<-10 m<-matrix(rnorm(n.vyber*n,mu,so),nrow=n.vyber) # průměry xm<-rowMeans(m) # histogram výběrových průměrů h.mean<-hist(xm,col=scales::alpha('blue',.5),probability=TRUE,breaks=15,add=TRUE) # proložíme hustotu pravděpodobnosti rozdělení průměrů lines(xs,dnorm(xs,mu,so/sqrt(n)),ty='l',col='blue',lwd=2) # 3) interval spolehlivosti pro průměr prumer<-mean(rnorm(n,mu,so)) alpha<-.05 is.dolni<-prumer+so/sqrt(n)*qnorm(alpha/2) is.horni<-prumer+so/sqrt(n)*qnorm(1-alpha/2) is<-c(is.dolni,is.horni) plot(xs,dnorm(xs,mu,so/sqrt(n)),ty='l',col='blue',lwd=2,main='Interval spolehlivosti pro střední hodnotu',xlab='glykémie [mmol/l]',ylab='hustota psti') lines(xs,dnorm(xs,mu,so),ty='l',col='red',lwd=2) lines(rep(mu,2),par('usr')[3:4]/c(1,10),col='blue',lwd=5) points(prumer,0,col='green',cex=3,pch=19) lines(is,rep(0,2),col='green',lwd=5) legend('topright',c('rozdělení glykémií','rozdělení průměrů','interval spolehlivosti','střední hodnota'),col=c('red','blue','green','blue'),lwd=c(2,2,5,5),pch=c(NA,NA,19,NA)) plot(xs,dnorm(xs,mu,so/sqrt(n)),ty='l',col='blue',lwd=2,main='Interval spolehlivosti + výběrové průměry',xlab='glykémie [mmol/l]',ylab='hustota psti') lines(xs,dnorm(xs,mu,so),ty='l',col='red',lwd=2) ys<-par('usr')[3:4]/c(1,20) for (i in 1:length(xm)) lines(rep(xm[i],2),ys,col='black') points(prumer,0,col='green',cex=3,pch=19) lines(is,rep(0,2),col='green',lwd=5) # kolik průměrů padne do intervalu spolehlivosti? n.prumer.v.is<-sum(xm>=is[1]&xm<=is[2]) legend('topright',c('rozdělení glykémií','rozdělení průměrů','interval spolehlivosti', 'výběrové průměry',paste0('(',n.prumer.v.is,' ze ',n.vyber,' v IS)')), col=c('red','blue','green','black','black'),lwd=c(2,2,5,1,NA),pch=c(NA,NA,19,NA,NA)) sum(xm>=is[1]&xm<=is[2]) savePlot('is_mu.png') # obnovíme původní grafické rozložení par(opar)