# Code von Franziska Wandtner, 2011, wandtner@math.uni-frankfurt.de

##########################################################################
# Wir betrachten eine Partition B_1=[0,1/3), B_2=[1/3,2/3), B_3=[2/3,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 gleicher Größe vorstellen, in die rein zufällig Kugeln abgelegt werden, und fragen:
# 1.) In welche Kisten fallen die Punkte? (Protokoll der Reihenfolge für die ersten 10 Kugeln)
# 2.) 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, auch "de Finetti Diagramm" genannt.


##########################################################################
rm(list=ls())

############################
# Anzahl der Wiederholungen
k<-10000

###########################################################################
# Graphikparameter werden gesetzt     
# Für Windows
windows(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,33.33,33.33,0),c(2,2,2.2,2.2),col="yellow")
polygon(c(33.33,66.66,66.66,33.33),c(2,2,2.2,2.2),col="green")
polygon(c(66.66,100,100,66.66),c(2,2,2.2,2.2),col="blue")
text(0,1.9,"0")
text(33.33,1.9,"1/3",cex=0.8)
text(66.66,1.9,"2/3",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]<1/3){n<-1}else{if(s[i]>=1/3 & s[i]<2/3){n<-2}else{if(s[i]>2/3){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]<1/3){count[1]<-count[1]+1}else{if(s[i]>=1/3 & s[i]<2/3){count[2]<-count[2]+1}else{if(s[i]>2/3){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")

readline("Balkendiagramm der ersten Besetzung")
par(font.main=1,mar=c(5,5,5,5)) 
plot(c(0,3.5),c(0,10),type="n",main=paste("Erste Besetzung:","Z = (",count[1],",",count[2],",",count[3],")"),axes=F,xlab="",ylab="Anzahl Kugeln")
axis(1,labels=c("Kiste Nr.1","Kiste Nr.2","Kiste Nr.3"),tick=F,at=c(0.75,1.75,2.75))
axis(2,at=c(1:10))
box()
polygon(c(0.5,0.5,1,1),c(0,count[1],count[1],0),col="yellow")
polygon(c(1.5,1.5,2,2),c(0,count[2],count[2],0),col="green")
polygon(c(2.5,2.5,3,3),c(0,count[3],count[3],0),col="blue")

readline("de Finetti Diagramm")
windows(width=7, height=6)
par(font.main=1,mar=c(2,2,2,2)) 

# Dreiecks-Diagramm vorbereiten
plot(c(-10,100),c(0,90),type="n",axes=F,xlab="",ylab="", main=paste("Erste Besetzung im de Finetti Diagramm"))
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.9)
text(-3,-1,expression(paste(Z[2])),cex=0.9)
text(103,-1,expression(paste(Z[3])),cex=0.9)

readline("Erste Besetzung eintragen")
	points(0.5*10*(count[1]+2*count[3]),10*sqrt(3)*count[1]/2,col="red",cex=2,pch=19)
	
text(7,85,"Erste Besetzung Z:",col="red")

text(-5,80,expression(paste(Z[1], " = ")), col="red")
text(-1,80.1,paste(count[1]),col="red")
text(6,80,expression(paste(Z[2], " = ")), col="red")
text(10,80.1,paste(count[2]),col="red")
text(17,80,expression(paste(Z[3], " = ")), col="red")
text(21,80.1,paste(count[3]),col="red")

vgl<-100*count[1]+10*count[2]+count[3]

readline("Wie typisch war diese Besetzung?")
plot(c(0,100),c(0,100),type="n",xlab="",ylab="",axes=F)
text(50,60,"Wie typisch war diese Besetzung?")

readline("Wiederhole k Mal")
text(50,50,paste("Wiederhole",k,"Mal!"), col="red")


###########################################################################
# Besetzungszahlen für viele Wiederholungen

readline("Simuliere k Mal")

windows(width=7, height=6)
par(font.main=1,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]<1/3){counthelp[j,1]<-counthelp[j,1]+1}else{if(s[i]>=1/3 & s[i]<2/3){counthelp[j,2]<-counthelp[j,2]+1}else{if(s[i]>2/3){counthelp[j,3]<-counthelp[j,3]+1}}}
	}
}

# Häufigkeiten


count<-matrix(rep(0,k),ncol=1)

for(i in 1:k)
{
	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

windows(width=7, height=6)
par(font.main=1,mar=c(2,2,2,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"))
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.9)
text(-3,-1,expression(paste(Z[2])),cex=0.9)
text(103,-1,expression(paste(Z[3])),cex=0.9)

	
# 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])
	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])
}


#Erste Besetzung
readline("Häufigkeit der ersten 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(7,85,"Erste Besetzung Z:",col="red")

text(-5,80,expression(paste(Z[1], " = ")), col="red")
text(-1,80.1,paste(punkt[1]),col="red")
text(6,80,expression(paste(Z[2], " = ")), col="red")
text(10,80.1,paste(punkt[2]),col="red")
text(17,80,expression(paste(Z[3], " = ")), col="red")
text(21,80.1,paste(punkt[3]),col="red")
###########################################################################
#Ende

readline("Ende")
graphics.off()