# LV Blatt 8
#
# Aufg. 8.2
#
# ###############################################################
#
# Name:      Wsk.komplette.Sammlung
#
# Aufgabe:   Berechnet die Wahrscheinlichkeit einer kompletten 
#            Sammlung von n.div Bildern beim Kauf von n.buy Bildern
#
# Parameter: n.div  Anzahl der verschiedenen Bilder
#            n.buy  Anzahl der gekauften Bilder        
#
# Rückgabewert: Wsk. einer kompletten Sammlung
#
# ################################################################

Wsk.komplette.Sammlung <- function( n.div, n.buy )
{
  return (1-sum ((-1)^(1:n.div+1)*choose(n.div,1:n.div)*(1 - 1:n.div/n.div)^n.buy))
}

# ###############################################################
#
# Name:      anz.Bilder.kompl.Sammlung   
#
# Aufgabe:   Liefert die Mindestanzahl der Bilder zurück, die man
#            kaufen muss, um mit einer Wahrscheinlichk. von p eine 
#            komplette Sammlung von n.div Bildern zu besitzen            
#
# Parameter: n.div  Anzahl der verschiedenen Bilder
#            p      gewünschte Wsk. für komplette Sammlung
#  
# Rückgabewert: Anzahl der zu kaufenden Bilder für kompl. Samml.
#
# ##############################################################
anz.Bilder.kompl.Sammlung <- function( n.div, p )
{
  n<-n.div # man muss mind. n.div Bilder kaufen
  while (Wsk.komplette.Sammlung(n.div, n) < p)
    n<-n+1
  
  return ( n )
}

# p=0,95
anz.Bilder.kompl.Sammlung(30,0.95)
# p=0,99
anz.Bilder.kompl.Sammlung(30,0.99)
#
# Näherungsweise Lsg. unter Verwendung der (nicht geg.) Unabh.
#
(log(1-0.95^(1/30)))/log(29/30)
#
# bei 660 Aufklebern
(log(1-0.95^(1/660)))/log(659/660)
# zum Vgl.
anz.Bilder.kompl.Sammlung(660,0.95)
# Formel versagt, da die Zahlen zu groß werden (siehe Bemerkung unten)!!!

# Probier's mal mit Simulation

# #####################################################################
#
# Name: first.Compl
#
# Aufgabe: ermittelt zu einer Bildersammlung den ersten Index, bei
#          dem die Sammlung komplett ist 
#
# Parameter: x      Vektor, der die Sammlung enthält 
#            n.div  Anzahl der verschiedenen Bilder
#
# Rückgabewert: Index, bei dem zum ersten Mal eine Sammlung kompl. ist
#
# ##############################################################
first.Compl <- function(x, n.div)
{
  return(max(match(c(1:n.div),x)))
}

# #####################################################################
#
# Name: anz.Bilder.kompl.Sammlung.sim
#
# Aufgabe: simuliert den n.sim-fachen Kauf von jeweils (hoffentlich) genuegend
#          vielen Bildern und berechnet damit einen Schätzwert für die
#          Anzahl der Bilder, bei denen man mit einer Wsk. von p eine
#          komplette Sammlung besitzt
#
# Parameter: n.div  Anzahl der verschiedenen Bilder
#            p      gewünschte Wsk. für komplette Sammlung
#            n.sim  Anzahl der simulierten Sammlungen
#
# Rückgabewert: Anzahl der zu kaufenden Bilder für kompl. Samml.
#
# ##############################################################
anz.Bilder.kompl.Sammlung.sim <- function( n.div, p, n.sim )
{
  bilder  <- matrix(sample(c(1:n.div),n.sim*2*round(log(1-p^(1/n.div))/log((n.div-1)/n.div)),
                           replace=TRUE),ncol=n.sim)
  sammlKompl  <- apply(bilder, 2, first.Compl,n.div = n.div) # je Spalte nachschaun, wann Sammlung komplett 
  
  # ersten Index in der empirischen Vtlgs.fkt. der Bilderanzhal für kompl. Sammlung,
  # bei dem p überschritten wurde, zurück geben
  return (quantile(ecdf(sammlKompl),p))
}

anz.Bilder.kompl.Sammlung.sim( 30, 0.95, 100000) # Hanutas
anz.Bilder.kompl.Sammlung.sim( 30, 0.99, 100000) # Hanutas
anz.Bilder.kompl.Sammlung.sim( 660, 0.95, 10000) # Panini
anz.Bilder.kompl.Sammlung.sim( 660, 0.99, 10000) # Panini

###########################################
#                                         #  
# Achtung                                 #
#                                         #
###########################################
#
# Die Funktion "Wsk.komplette.Sammlung" liefert seltsame Ergebnisse, z.B.
#
# Wsk.komplette.Sammlung(30,2)=1.142855e-08 (sollte nat. 0 sein)
#
# und 
#
# Wsk.komplette.Sammlung(30,30)=6.00231e-12 (hier sollte 30!/30^30=1.288e-12 raus kommen)
#
# und 
#
# Wsk.komplette.Sammlung(660,660)=5.018835e+64 (sollte 660!/660^660=1.4947e-285 sein)
#
# Das Problem ist die interne Genauigkeit, mit der R rechnet, naemlich 16 gueltige
# Ziffern. Da in der alternierenden Summe Summanden vorkommen, die viel
# groesser als das Endergebnis sind, wie z.B. fuer n=2 und k=13:
# (-1)^(k+1)*choose(30,k)*(1 - k/30)^n=3.845622e+07
# schlaegt sich eine Ungenauigkeit an der 16. Stelle im Endergebnis bereits an
# der neunten Nachkommastelle nieder (weil sieben Stellen vor dem Komma stehen).
# 
# Um diese numerische Ungenauigkeit zu vermeiden, 
# muss man die Package gmp verwenden, die mit beliebig grossen, 
# nur durch den Speicherplatz beschraenkten, Integers rechnen kann 
# (fuer die uebungsaufgabe selbst spielt die Ungenauigkeit keine Rolle, da es
#  es in der Uebungsaufgabe um eine Wahrscheinlichkeit von 95% geht)
#
#
# Hinweis: Das Paket Rmpfr ermöglicht es, die interne Rechengenauigkeit 
#          in bits anzugeben. Der Standard betraegt 53 bits, das sind ungefaehr
#          16 gueltige Stellen. Erhoeht man diese Genauigkeit in der
#          Funktion kompletteSammlung durch "k/30" durch "k/mpfr(30,74)"
#          und zwingt R damit auf 74 bits, das entspricht ungefaehr 22 gueltigen
#          Stellen, geanau zu rechnen, so erhaelt man ebenfalls sehr gute
#          Ergebnisse (nicht ganz 0 für n<30, aber wesentlich kleiner als
#          alle auftretenden Wahrscheinlichkeiten)

library("Rmpfr")  # lädt die Pakete Rmpfr UND gmp; bei Bedarf vorher mit 
# install.packages("gmp") und install.packages("Rmpfr") 
# installieren

# Beschreibung wie bei Wsk.komplette.Sammlung

Wsk.komplette.Sammlung.bigInt <- function( n.div, n.buy )
{
  anzMgl <- as.bigz(n.div)^n.buy;  # Gesamtzahl an Moeglichkeiten (intern ohne 
  # Rundung dargestellt)
  
  # Moeglichkeiten fuer Sammlungen, die nicht komplett sind, abziehen
  # (Formel nach Sylvester-Poincare)
  # Beachte: choose funlktioniert nicht korrekt für große Zahlen
  for (k in 1:n.div)
  {
    anzMgl <- anzMgl - (-1)^(k+1)*
      as.bigz(factorial(as.bigz(n.div))*
                as.bigz(n.div - k)^n.buy/(factorial(as.bigz(n.div-k))*factorial(as.bigz(k))))
  }
  
  erg <- mpfr(anzMgl/as.bigz(n.div)^n.buy,precBits = 100); # zum Schluss in Dezimalbruch umwandeln
  # (dabei muss die Package mpfr verwendet werden,
  # da die normale Rechengenauigkeit nicht genuegt)
  
  return (erg)
}

# jetzt passt alles
for (n in 1:30) print(Wsk.komplette.Sammlung.bigInt(30,n))

Wsk.komplette.Sammlung.bigInt(660,660)
# das entspricht 660!/660^660
mpfr(factorial(as.bigz(660))/as.bigz(660)^660,precBits = 100)

anz.Bilder.kompl.Sammlung.bigInt <- function( n.div, p )
{
  n<-n.div # man muss mind. n.div Bilder kaufen
  while (Wsk.komplette.Sammlung.bigInt(n.div, n) < p)
    n<-n+20 # 20er Schritte, weil es sonst zu lang dauert
  
  return ( n )
}

# Ab wann hat man zu 50% eine komplette Sammlung?
anz.Bilder.kompl.Sammlung.bigInt(660, 0.5)   # dauert ein bisschen
(log(1-0.5^(1/660)))/log(659/660) # zugehörige Näherung

anz.Bilder.kompl.Sammlung.bigInt(660, 0.95)  # dauert noch länger
anzBilder.bigInt(660, 0.99)  # und noch länger