###################################################################################################### # Ziehen ohne Zurücklegen und die Normalapproximation ###################################################################################################### # Code von Benjamin Straub, modifiziert von AW unter Verwendung eines Hinweises von Igor Scheller # zur Vereinheitlichung der Fenstersetzung für Windows, Linux und MacOS ###################################################################################################### # hier lassen sich die Parameter setzen ---- N<-1000 # Anzahl der Simulationen zum Schätzen vom Stchprobenmittelwert M g<-90 # Anzahl Individuen in der Population n<-30 # Wie oft soll aus der Population ohne Zurücklegen gezogen werden ###################################################################################################### ###################################################################################################### # Teil 1: Die Verteilung des Stichprobenmittelwerts ---- ###################################################################################################### # Graphikparameter werden gesetzt if (Sys.info()["sysname"] == "Windows") { windows(width=8, height=6, pointsize=12) } else if (Sys.info()["sysname"] == "Linux") { x11(width=8, height=6) } else { system('rm plotStofiAufgabe30_*.pdf') pdf('plotStofiAufgabe30_1.pdf', width=16, height=12) } par(mar = c(4, 5, 1.5, 1) + 0.6) par(mfrow = c(1,1)) par(cex.lab=2) par(cex=2) par(mgp=c(3.2,1.2,0)) # Definieren der Population n_10<-45 # Größe 10 n_5<-15 # Größe 5 n_20<-30 # Größe 20 population<-c(rep(10,n_10),rep(5,n_5),rep(20,n_20)) # Unsere Population der Größen ###################################################################################################### # Simulation der Verteilung des Stichprobenmittels M:=1/k*(X1+...+Xk) M<-rep(0,N) for (i in 1:N){ X<-sample(population,n,replace=FALSE) M[i]<-mean(X) } # Erstellen Histogramm hist(M,probability=TRUE,xlim=c(10,15),ylim=c(0,1.4),ylab="Dichte",cex=0.4, main="Histogramm von M") mu_M<-mean(population) abline(v=mu_M,col="blue",lwd=3) text(mu_M,1.3,bquote(mu[M]==.(mu_M)),pos=4,col="blue",cex=1) # Normalapproximation MIT 'endlicher Population' Korrektur x<-seq(10,15,0.01) var_M<-1/n*(var(population)-var(population)/(length(population)-1)*(n-1)) lines(x,dnorm(x,mu_M,sqrt(var_M)),col="red",lwd=3) text(10.0,1.3,"Normalapproximation",pos=4, col="red",cex=0.8) text(10.2,1.1,"mit",pos=4, col="red",cex=0.8) text(9.9,0.9,"Varianz-Korrektur",pos=4, col="red",cex=0.8) # Normalapproximation OHNE 'endlicher Population' Korrektur var_M_ohne_korrektur<-1/n*var(population) lines(x,dnorm(x,mu_M,sqrt(var_M_ohne_korrektur)),col="orange",lwd=3) text(10.65,1.1,"/ohne",pos=4, col="orange",cex=0.8) ###################################################################################################### # Einzeichnen der asymptotischen und des empirischen 95% Konfidenzintervalle für M arrows(qnorm(0.975,mu_M,sqrt(var_M)),0,qnorm(0.975,mu_M,sqrt(var_M)),0.4,length = 0, angle = 30, code = 3, col ='red',lwd=3,lty=3) arrows(qnorm(0.0275,mu_M,sqrt(var_M)),0,qnorm(0.0275,mu_M,sqrt(var_M)),0.4,length = 0, angle = 30, code = 3, col ='red',lwd=3,lty=3) arrows(qnorm(0.975,mu_M,sqrt(var_M_ohne_korrektur)),0,qnorm(0.975,mu_M,sqrt(var_M_ohne_korrektur)),0.4,length = 0, angle = 30, code = 3, col ='orange',lwd=3,lty=3) arrows(qnorm(0.0275,mu_M,sqrt(var_M_ohne_korrektur)),0,qnorm(0.0275,mu_M,sqrt(var_M_ohne_korrektur)),0.4,length = 0, angle = 30, code = 3, col ='orange',lwd=3,lty=3) arrows(quantile(M,0.975),0,quantile(M,0.975),0.4,length = 0, angle = 30, code = 3, col ='black',lwd=3,lty=3) arrows(quantile(M,0.025),0,quantile(M,0.025),0.4,length = 0, angle = 30, code = 3, col ='black',lwd=3,lty=3) text(11.45,0.7,"Empir. Quantile",pos=4, col="black",cex=0.8) text(quantile(M,0.975)-0.1,0.43,"97.5 %",pos=4, col="black",srt=90) text(quantile(M,0.025)-0.1,0.43,"2.5 %",pos=4, col="black",srt=90) ###################################################################################################### ###################################################################################################### ###################################################################################################### # Teil 2 - Überdeckungswkt. des Kofidenzintervalls ----mit "endliche Population"-Korrektur ###################################################################################################### # Wie oft liegt der Erwartungswert mu in dem geschätzten Konfidenzintervall # mit "endliche Population"-Korrektur if (Sys.info()["sysname"] == "Windows") { windows(width=10, height=6, pointsize = 12) } else if (Sys.info()["sysname"] == "Linux") { x11(width=10, height=10) } else { dev.off() pdf('plotStofiAufgabe30_2.pdf', width=16, height=12) } par(mar = c(4, 5, 1.5, 1) + 0.6) par(mfrow = c(1,1)) par(cex.lab=2) par(cex=2) par(mgp=c(3.2,1.2,0)) ###################################################################################################### # Plot der Konfidenzintervalle mit 'endlicher Population' Korrektur plot("",pch=20,ylim=c(7,17),xlim=c(1,N),xlab="Simulation",ylab="Intervall I",main="Mit 'endliche Population'-Korrektur") links_raus<-which(M>mu_M+2*sqrt(var_M)) rechts_raus<-which(Mmu_M+2*sqrt(var_M_ohne_korrektur)) rechts_raus<-which(M0){ drinnen<-c(1:N)[-c(links_raus,rechts_raus)] } for (i in drinnen){ arrows(i,M[i]-2*sqrt(var_M_ohne_korrektur),i,M[i]+2*sqrt(var_M_ohne_korrektur),length = 0, angle = 30, code = 3, col ='black',lwd=1) } for (i in links_raus){ arrows(i,M[i]-2*sqrt(var_M_ohne_korrektur),i,M[i]+2*sqrt(var_M_ohne_korrektur),length = 0, angle = 30, code = 3, col ='red',lwd=3) } for (i in rechts_raus){ arrows(i,M[i]-2*sqrt(var_M_ohne_korrektur),i,M[i]+2*sqrt(var_M_ohne_korrektur),length = 0, angle = 30, code = 3, col ='green',lwd=3) } abline(h=mu_M,col="blue",lwd=3) text(N/20,8,expression(P(mu%notin%I))) text(N/5,8,"=") text(N/3.5,8,round(length(links_raus)/N,3),col="red") text(N/2.5,8,"+") text(N/2,8,round(length(rechts_raus)/N,3),col="green") text(N/1.6,8,"=") text(N/1.3,8,round((length(rechts_raus)+length(links_raus))/N,3)) if (Sys.info()["sysname"] == "Windows") { } else if (Sys.info()["sysname"] == "Linux") { } else { dev.off() system('open plotStofiAufgabe30_*.pdf') } cat("Zum Loeschen der Graphiken Entder druecken"); readline("..."); dummy <- readLines("stdin",n=1); #Befehl bei R-Studio löschen! if (Sys.info()["sysname"] == "Windows") { graphics.off() } else if (Sys.info()["sysname"] == "Linux") { graphics.off() } else { system('rm plotStofiAufgabe30_*.pdf') } cat("Graphiken erfolgreich entfernt. "); ###################################################################################################### ######################################################################################################