#Blatt 4
#a)
getinhomogenerPP<-function(window_x_lim, window_y_lim,supremum,funktion){
#Realisierung homogener Poissonprozess
PP<-getPoissonPP(window_x_lim, window_y_lim, supremum)
#Anzahl der Punkter der Realisierung
n<-length(PP[,1])
#Ausdünnung
inhomPP<-c()
u<-runif(n,0,1)
for(i in 1:n){
if(u[i]<=funktion(PP[i,])/supremum){
inhomPP<-rbind(inhomPP,PP[i,])
}
}
return(inhomPP)
}
##########################################################################
#b)
#Funktion implementieren
givenfunction<-function(x){
return(x[1]*x[2]*10^-5)
}
supremum<-0.1
liste<-list()
for(i in 1:10){
liste[[i]]<-getinhomogenerPP(c(0,100),c(0,100),supremum,givenfunction)
}
##########################################################################
#c)
#window_x_lim, window_y_lim Beobachtungsfenster
#PP beobachteter Poissonprozess
#x beobachteter Punkt, h Bandbreite
lambdadach<-function(window_x_lim, window_y_lim, PP, x, h){
#Anzahl der Punkte der Realisierung
n<-length(PP[,1])
#Zählvariable
count<-0
#Anzahl der Punkte im Ausschnittsfenster
for(i in 1:n){
#Überprüfen ob Punkt im Ausschnittsfenster
if(abs(PP[i,1]-x[1])<=h/2&&abs(PP[i,2]-x[2]<=h/2)){
count<-count+1
}
}
#Volumen des Ausschnittsfensters
x_start<-x[1]-h/2
x_end<-x[1]+h/2
y_start<-x[2]-h/2
y_end<-x[2]+h/2
#Fallunterscheidungen zur Berechnung des Volumen
if(x_start<window_x_lim[1]){
x_start<-window_x_lim[1]
}
if(x_end>window_x_lim[2]){
x_end=window_x_lim[2]
}
if(y_start<window_y_lim[1]){
y_start<-window_y_lim[1]
}
if(y_end>window_y_lim[2]){
y_end=window_y_lim[2]
}
#Volumen
v<-(x_end-x_start)*(y_end-y_start)
return(count/v)
}


##############################################################################
#d)
#Ergebnisvektor
result<-c()
#Ablaufen aller Realisierungen
for(i in 1:10){
#Anzahl der Punkte in der jeweiligen Realisierung
n<-length(liste[[i]][,1])
#Initialwerte für Schätzer von h und Loglikelihoodfkt.
hopt<-0
logmax<--Inf
#Auswahl für hopt aus geg. Menge
for(j in 11:30){
#Hilfsvektor zur Speicherung der jeweiligen Werte
lambdaresult<-c()
#Ablaufen aller Punkte
for(k in 1:n){
if(k>1 & k<n){
lambdaresult<-cbind(lambdaresult,lambdadach(c(0,100),c(0,100),rbind(liste[[i]][1:(k-1),],liste[[i]][(k+1):n,]), liste[[i]][k,], j))
}
else{
if(k==1){
lambdaresult<-cbind(lambdaresult,lambdadach(c(0,100),c(0,100),liste[[i]][(k+1):n,], liste[[i]][k,], j))
}
else{
lambdaresult<-cbind(lambdaresult,lambdadach(c(0,100),c(0,100),liste[[i]][1:(k-1),], liste[[i]][k,], j))
}
}
if(sum(log(lambdaresult))>logmax){
hopt<-j
logmax<-sum(log(lambdaresult))
}
}
}
result[i]<-hopt
}
###############################################################################
#e)
#Ablaufen aller Realisierungen
resultmatrix<-matrix(0,10,19)
for(i in 1:10){
#Ablaufen aller h
for(j in 11:30){
#fehler initialisierung
sqerror<-0
#Ablaufen des Fensters
for(k in 0:99){
for(l in 0:99){
sqerror<-sqerror + ( lambdadach(c(0,100),c(0,100), liste[[i]], c((k+1/2),(l+1/2)), j)- givenfunction(c((k+1/2), (l+1/2))) )^2
}
}
#normierung des fehlers +  speichern
resultmatrix[i,(j-11)]<-sqerror/(10000)
}
}