0

I want to calculate the integral of product of three kernel density functions. In order to do that, after finding the kernel densities, I should find an approximate function for the product of them and then calculate the integral (with parametric upper limit). Please see the code below that gives me an error. I think there should be a problem in calculating the approximate function or the integration part..

library(ks)
library(rgl)

zz <- " longitude  latitude       depth      time       magnitude
       277.728371  139.925845  0.75103658  26.40786514 0.943718276
       426.087586  168.903095  0.2957441   0.241456485 0.759137864
       331.549444  74.168092   0.55140397  66.51363095 0.776176433
       93.078983   78.588053   0.15328453  104.9418546 0.834896464
       492.359229  11.082291   0.08173915  111.3391451 0.874798119
       86.85704    42.34973    0.23081904  152.8098572 0.878111793
       128.038949  73.935782   0.66160123  157.8933315 0.990100773
       295.300125  1.935765    0.49789785  159.9134319 0.842815655
       294.688309  1.024583    0.44789667  165.7092358 0.886545275
       221.246937  151.217171  0.6337224   167.6213491 0.885163617
       111.240376  156.04214   0.55752237  171.2039395 0.885273526
       25.929383   136.975153  0.0271747   172.6574772 0.812214826
       415.726989  158.482975  0.37340509  184.3767148 0.90837174
       73.921877   60.031908   0.15224511  189.6429637 0.791403228
       8.811256    124.676545  0.26806101  193.7498013 0.813638308"

x <- read.table(text=zz, header = TRUE)
y1 <- x[,1:3]
y2<- x[ ,4]
y3<- x[ ,5]
fhat1 <- kde(x=y1)
fhat2 <- kde(x=y2)
fhat3 <- kde(x=y3)

integrand <- function(p){fhat1(p,y1)* fhat2(p,y2)* fhat3(p,y3)}
Vintegrand <- Vectorize(integrand)
integrate(Vintegrand, lower = 0, upper = t)
Ferra
  • 9
  • 4
  • See page 22 of [this](http://research.cs.tamu.edu/prism/lectures/pr/pr_l7.pdf). – Alex A. Apr 07 '15 at 16:32
  • @AlexA. Thank you Alex. That is exactly what I am doing. I have a multivariate kernel and I am trying to find the total density function by using the product of them and then calculate the integral of the density function. Do you know what is the problem with my code? – Ferra Apr 07 '15 at 16:57
  • It looks like you should be dividing the product by something. – Alex A. Apr 07 '15 at 16:59
  • except that part, I think there is a problem in integration part. Either the calculated approximate function is not correct or the integration part. – Ferra Apr 07 '15 at 17:01
  • I mean that the integrand--the product of the three KDEs--should be divided by something as is done in the link I posted. That's the approximate function you're talking about, right? – Alex A. Apr 07 '15 at 17:03
  • You mean divided by the bandwidths? That's right, but that will not cause the error in this code. I mean, there should be a problem in the way or functions that I used... – Ferra Apr 07 '15 at 17:06
  • You said in the post that your problem was finding the right approximate function. But now you're saying that the problem is that you're getting an error? If that's the case, please include that information in the question. – Alex A. Apr 07 '15 at 17:14
  • @AlexA.You are right Alex. I edited the post. Thank you. – Ferra Apr 07 '15 at 17:16
  • 1
    What is the error you're getting? – Alex A. Apr 07 '15 at 17:33
  • @AlexA. When I use the code as it is it gives me this error: Error in is.finite(upper) : default method not implemented for type 'closure'. I thought it's because it can not calculate the integral with parametric bounds but when I change the upper limit t to for example 1000 it give me this error: Error in (function (p) : could not find function "fhat1" – Ferra Apr 07 '15 at 17:35
  • That's because `fhat1` wasn't declared to be a function. Take a look at the docs: `?kde`. – Alex A. Apr 07 '15 at 17:43
  • So, how can I calculate the approximate function for fhat1 ? Also I need the upper limit to be t, can I use integrate function or it doesn't work with parametric limits? – Ferra Apr 07 '15 at 17:46
  • Given the confusion you are exhibiting about fairly basic R, I'm wondering if your problem definition is also confused. There is no sensible interpretation I can think of for the product of densities. I'm wondering if you really need to do the integration first to get cumulative probabilities. – IRTFM Apr 07 '15 at 17:59
  • @Ferra: `integrate()` performs adaptive quadrature, i.e. numerical integration. It doesn't return a function. `t` needs to evaluate to a numeric constant for `integrate()` to work. – Alex A. Apr 07 '15 at 18:01
  • @BondedDust: Agreed about the problem, but there is a sensible interpretation; see my first comment on the post. – Alex A. Apr 07 '15 at 18:03
  • @BondedDust As alex showed in that pdf, one of the ways to calculate the multivariate kernel density is to calculate the kernel density for each of the variables separately and then find the product of them (if they are independent). I need to find the integration of the calculated density function in order to use it in another calculation. – Ferra Apr 07 '15 at 18:09
  • You are offering `fhat1`, `fhat2` and `fhat3` as functions when you wrote: `integrand <- function(p){fhat1(p,y1)* fhat2(p,y2)* fhat3(p,y3)}`. Since they are not functions, but rather lists , you get an error. Furthermore, you offered a three column matrix to fhat1 and only one column matrices to the other two kde derivations. That cannot be correct. – IRTFM Apr 07 '15 at 18:21
  • @BondedDust About fhat1 you are right, I can divide it into 3 separate kernel densities. Regarding the functions that you said they are not function and they are lists, how can I define them as functions? Also I need to calculate the indefinite integration. I searched online and it seems there is no function in R for calculation indefinite integrals. If so, please tell me to switch to Matlab instead. Thanks. – Ferra Apr 07 '15 at 18:30
  • There is `approxfun` and many spline-function-creation tools available What do you mean when you say "need to calculate the indefinite integration"? R doesn't do symbolic calculus. – IRTFM Apr 07 '15 at 18:35
  • @BondedDust, Therefore, I cannot use R in this case (actually the upper bound limit of the integral is the parameter that I want to estimate in my problem) – Ferra Apr 07 '15 at 19:35
  • @BondedDust: Totally an aside, but R actually _can_ do symbolic derivatives in base R (see `?deriv`) but not symbolic integration without an external library, e.g. Ryacas. – Alex A. Apr 07 '15 at 19:44
  • @AlexA. Right. I guess I should not be too categorical on that matter. R also has the `cubature` package that can do multivariate integration, but as for whether it can assist with this problem still requires a better understanding of exactly what the scientific "problem" might be. Haven't figured out whether the problem involves a path integral or a volume integral. Needs some sort of dimensional analysis to see if these disparate measurements actually make sense when multiplied and then summed. – IRTFM Apr 07 '15 at 19:57
  • @BondedDust: I think what's happening here is that it's creating a kind of mixture of univariate empirical distributions and the OP wants to construct a function that will return Pr[0 – Alex A. Apr 07 '15 at 20:09
  • That's the mathematics. But will the numbers actually have a physical interpretation? – IRTFM Apr 07 '15 at 20:15

0 Answers0