# Code von Franziska Wandtner, 2011, wandtner@math.uni-frankfurt.de ########################################################################## # Wir betrachten eine Partion B_1=[0,3/10), B_2=[3/10,8/10), B_3=[8/10,1] des Intervalls [0,1] und fragen bei 10 Zufallsvariablen, die uniform auf [0,1] verteilt sind, in welches Teilstück B_i sie gefallen sind. # Wir können uns 3 Kisten unterschiedlicher Größe vorstellen, in die rein zufällig Kugeln abgelegt werden, und fragen: # Wieviele Kugeln fallen in die einzelnen Kisten? ("Besetzungszahlen" Z1, Z2, Z3 für die ersten 10 Kugeln) # Zunächst schauen wir uns die Ausgänge für die ersten 10 Kugeln einzeln an. # Anschließend wiederholen wir den Vorgang k=10000 Mal und veranschaulichen das Ergebnis mit einem Dreiecksdiagramm. # Schließlich vergleichen wir das Ergebnis mit den exakten Gewichten der Multinomialverteilung. ########################################################################## rm(list=ls()) ############################ # Anzahl der Wiederholungen k<-10000 ########################################################################### # Graphikparameter werden gesetzt # Für Mac quartz(width=10, height=6) par(cex.main=1.3,cex.axis=1.0,cex.lab=1.0,pch=18) par(font.main=1,mar=c(2,2,5,2)) ########################################################################### # Kugeln einzeln plot(c(0,100),c(0,2.5),type="n",xlab="",ylab="",axes=F,main=paste("Wohin fallen die ersten 10 Punkte?")) polygon(c(0,100,100,0),c(2,2,2.2,2.2)) polygon(c(0,30,30,0),c(2,2,2.2,2.2),col="yellow") polygon(c(30,80,80,30),c(2,2,2.2,2.2),col="green") polygon(c(80,100,100,80),c(2,2,2.2,2.2),col="lightblue") text(0,1.9,"0") text(30,1.9,"3/10",cex=0.8) text(80,1.9,"8/10",cex=0.8) text(100,1.9,1) text(1,1.5,"Kistennummer:",adj=c(0,0)) readline("Erster Punkt") s<-runif(10) n<-0 for(i in 1:10) { points(s[i]*100,2.1) if(s[i]<3/10){n<-1}else{if(s[i]>=3/10 & s[i]<8/10){n<-2}else{if(s[i]>8/10){n<-3}}} text(1+i*5,1.3,n) readline("Nächster Punkt") } count<-rep(0,3) for(i in 1:10) { if(s[i]<3/10){count[1]<-count[1]+1}else{if(s[i]>=3/10 & s[i]<8/10){count[2]<-count[2]+1}else{if(s[i]>8/10){count[3]<-count[3]+1}}} } readline("Besetzungszahlen") text(1,1,"Besetzungszahlen:",adj=c(0,0),col="red") text(16,0.8,expression(paste(Z[1], " = ")), col="red") text(19,0.81,paste(count[1]),col="red") text(49,0.8,expression(paste(Z[2], " = ")), col="red") text(52,0.81,paste(count[2]),col="red") text(82,0.8,expression(paste(Z[3], " = ")), col="red") text(85,0.81,paste(count[3]),col="red") vgl<-100*count[1]+10*count[2]+count[3] ########################################################################### # Besetzungszahlen für viele Wiederholungen readline("de Finetti für k Wiederholungen") quartz(width=6, height=6) par(mar=c(2,2,2,2)) k<-k # Anzahl Wiederholungen counthelp<-matrix(rep(0,k*3),ncol=3) for(j in 1:k) { s<-runif(10) for(i in 1:10) { if(s[i]<3/10){counthelp[j,1]<-counthelp[j,1]+1}else{if(s[i]>=3/10 & s[i]<8/10){counthelp[j,2]<-counthelp[j,2]+1}else{if(s[i]>8/10){counthelp[j,3]<-counthelp[j,3]+1}}} } } # Häufigkeiten #Achtung:Z2 kann 10 sein!! count<-matrix(rep(0,k),ncol=1) for(i in 1:k) { if(counthelp[i,1]==10){count[i]<-0000}else{if(counthelp[i,2]==10){count[i]<-1111}else{if(counthelp[i,3]==10){count[i]<-2222}else{count[i]<-counthelp[i,1]*100+counthelp[i,2]*10+counthelp[i,3]}}} } x<-as.data.frame(table(count)) y<-as.matrix(x) l<-length(y)/2 # Dreiecks-Diagramm vorbereiten plot(c(-10,100),c(0,90),type="n",axes=F,xlab="",ylab="", main=paste("Häufigkeiten der Besetzungen bei", k,"Wdh"),cex.main=0.8) lines(c(0,100),c(0,0),col="grey") lines(c(0,50),c(0,sqrt(3)*100/2),col="grey") lines(c(50,100),c(sqrt(3)*100/2,0),col="grey") for(i in 1:9) { lines(c(5*i,100-5*i),c(sqrt(3)*10*i/2,sqrt(3)*10*i/2),col="grey") lines(c(5*i,10*i),c(sqrt(3)*10*i/2,0),col="grey") lines(c(10*i,50+5*i),c(0,sqrt(3)*10*(10-i)/2),col="grey") } text(50,89,expression(paste(Z[1])),cex=0.7) text(-3,-1,expression(paste(Z[2])),cex=0.7) text(103,-1,expression(paste(Z[3])),cex=0.7) l<-length(y)/2 # Punkte eintragen punkte<-matrix(rep(0,3*l),ncol=3) g<-grey.colors(max(x[,2]+1),start=0,end=0.75) for (i in 1:l) { a<-as.numeric(y[i,1]) freq<-as.numeric(y[i,2]) if(a!=0000&a!=1111&a!=2222) {punkte[i,1]<-trunc(a/100,1) punkte[i,2]<-trunc((a-100*punkte[i,1])/10,1) punkte[i,3]<-trunc(a-100*punkte[i,1]-10*punkte[i,2],1) text(0.5*10*(punkte[i,1]+2*punkte[i,3]),10*sqrt(3)*punkte[i,1]/2,paste(freq),col=g[max(x[,2])+1-freq]) } else{if(a==0000){text(0.5*10*(10+2*0),10*sqrt(3)*10/2,paste(freq),col=g[max(x[,2])+1-freq]) }else{if(a==1111){text(0.5*10*(0+2*0),10*sqrt(3)*0/2,paste(freq),col=g[max(x[,2])+1-freq])}else(if(a==2222){text(0.5*10*(0+2*10),10*sqrt(3)*0/2,paste(freq),col=g[max(x[,2])+1-freq])})}} } #Erste Besetzung markieren readline("Erste Besetzung") erste<-as.numeric(x[x[1]==vgl]) punkt<-rep(0,3) punkt[1]<-trunc(erste[1]/100,1) punkt[2]<-trunc((erste[1]-100*punkt[1])/10,1) punkt[3]<-trunc(erste[1]-100*punkt[1]-10*punkt[2],1) text(0.5*10*(punkt[1]+2*punkt[3]),10*sqrt(3)*punkt[1]/2,paste(erste[2]),col="red") text(-5,80,expression(paste(p[1], " = 0.3 ")), col="red",cex=0.8) text(10,80,expression(paste(p[2], " = 0.5")), col="red", cex=0.8) text(25,80,expression(paste(p[3], " = 0.2")), col="red", cex=0.8) ########################################################################### #Vergleich mit exakten Multinomialgewichten readline("Vergleich mit erw. Häufigkeiten") quartz(width=11, height=6) par(mfrow=c(1,2)) par(mar=c(1,1,3,1)) # Dreiecks-Diagramm vorbereiten plot(c(-10,100),c(0,90),type="n",axes=F,xlab="",ylab="", main=paste("Häufigkeiten der Besetzungen bei", k,"Wdh"),cex.main=0.8) lines(c(0,100),c(0,0),col="grey") lines(c(0,50),c(0,sqrt(3)*100/2),col="grey") lines(c(50,100),c(sqrt(3)*100/2,0),col="grey") for(i in 1:9) { lines(c(5*i,100-5*i),c(sqrt(3)*10*i/2,sqrt(3)*10*i/2),col="grey") lines(c(5*i,10*i),c(sqrt(3)*10*i/2,0),col="grey") lines(c(10*i,50+5*i),c(0,sqrt(3)*10*(10-i)/2),col="grey") } text(50,89,expression(paste(Z[1])),cex=0.7) text(-3,-1,expression(paste(Z[2])),cex=0.7) text(103,-1,expression(paste(Z[3])),cex=0.7) l<-length(y)/2 # Punkte eintragen punkte<-matrix(rep(0,3*l),ncol=3) g<-grey.colors(max(x[,2]+1),start=0,end=0.75) for (i in 1:l) { a<-as.numeric(y[i,1]) freq<-as.numeric(y[i,2]) if(a!=0000&a!=1111&a!=2222) {punkte[i,1]<-trunc(a/100,1) punkte[i,2]<-trunc((a-100*punkte[i,1])/10,1) punkte[i,3]<-trunc(a-100*punkte[i,1]-10*punkte[i,2],1) text(0.5*10*(punkte[i,1]+2*punkte[i,3]),10*sqrt(3)*punkte[i,1]/2,paste(freq),col=g[max(x[,2])+1-freq]) } else{if(a==0000){text(0.5*10*(10+2*0),10*sqrt(3)*10/2,paste(freq),col=g[max(x[,2])+1-freq]) }else{if(a==1111){text(0.5*10*(0+2*0),10*sqrt(3)*0/2,paste(freq),col=g[max(x[,2])+1-freq])}else(if(a==2222){text(0.5*10*(0+2*10),10*sqrt(3)*0/2,paste(freq),col=g[max(x[,2])+1-freq])})}} } # Dreiecks-Diagramm vorbereiten plot(c(-10,100),c(0,90),type="n",axes=F,xlab="",ylab="", main=paste("Erwartete Häufigkeiten der Besetzungen bei", k,"Wdh (gerundet)"),cex.main=0.8) lines(c(0,100),c(0,0),col="grey") lines(c(0,50),c(0,sqrt(3)*100/2),col="grey") lines(c(50,100),c(sqrt(3)*100/2,0),col="grey") for(i in 1:9) { lines(c(5*i,100-5*i),c(sqrt(3)*10*i/2,sqrt(3)*10*i/2),col="grey") lines(c(5*i,10*i),c(sqrt(3)*10*i/2,0),col="grey") lines(c(10*i,50+5*i),c(0,sqrt(3)*10*(10-i)/2),col="grey") } text(50,89,expression(paste(Z[1])),cex=0.7) text(-3,-1,expression(paste(Z[2])),cex=0.7) text(103,-1,expression(paste(Z[3])),cex=0.7) punktem<-matrix(rep(0,66*3),ncol=3) pk<-1 for(i in 0:10) { j<-0 while(i+j<=10){ punktem[pk,1]<-i punktem[pk,3]<-j j<-j+1 pk<-pk+1 } } for (i in 1:66) { text(0.5*10*(punktem[i,1]+2*punktem[i,3]),10*sqrt(3)*punktem[i,1]/2,paste(round(dmultinom(c(punktem[i,1],10-punktem[i,1]-punktem[i,3],punktem[i,3]),prob=c(0.3,0.5,0.2))*k,0)),col=g[max(x[,2])+1-round((dmultinom(c(punktem[i,1],10-punktem[i,1]-punktem[i,3],punktem[i,3]),prob=c(0.3,0.5,0.2))*k),0)]) } ########################################################################### #Ende readline("Ende") graphics.off()