Hej. Chcę zrobić estymator Parzena-Rosenblatta dla rozkładu dwuwymiarowego z jądrem normalnym w
\(\displaystyle{ \RR}\). Wzór nań to
\(\displaystyle{ \hat f = \frac{1}{nh^d} \sum_{i=1}^n K\left( \frac{x-x_i}{h}\right)}\), gdzie
\(\displaystyle{ d}\) to liczba wymiarów (dwa),
\(\displaystyle{ K}\) to jądro (normalne),
\(\displaystyle{ n}\) to ilość obserwacji, a
\(\displaystyle{ x_i}\) to te obserwacje. Napisałem funkcję kilka dni temu i od tych kilku dni temu szukam w niej błędu, którego znaleźć nie mogę. Funkcja bierze
\(\displaystyle{ x}\),
\(\displaystyle{ y}\) (lub wektor iksów i igreków) oraz ramkę danych i zwraca wartość estymatora w punkcie
\(\displaystyle{ (x,y)}\) (lub wektor wartości dla kolejnych punktów, które zostały wprowadzone jako argumenty funkcji w postaci wektora). Czy ktoś widzi błąd? :/ (Sposób wyboru
\(\displaystyle{ h}\) jest nieważny, to po prostu ma być jakaś mała liczba. Dvnorm to gęstość normalna w 2d).
Kod: Zaznacz cały
library(mvtnorm)
estymatorpr<-function(x,y,dane) {
n<-nrow(dane)
liczbaargumentow<-length(x)
h1<-bw.SJ(dane[,1])
h2<-bw.SJ(dane[,2])
h<-mean(c(h1,h2))
wektorwynikowy<-0
for(j in 1:liczbaargumentow) {
suma<-0
for(i in 1:n) {
suma<-suma+dmvnorm((c(x[j],y[j])-c(dane[i,][[1]],dane[i,][[2]]))/h)
}
wektorwynikowy[j]<-suma/(n*h^2)
}
wektorwynikowy
}
