How to calculate the integral of the sum of two kernel density functions in R?

Ferra

New Member
#1
I want to calculate the integral (0 to t) of the sum of two kernel density functions in R. I wrote the code as follows, but it gives me error. I am sure about the kernel density estimation part, but it seems the summation or the integral part is not correct. I appreciate any help in advance.


Code:
library(ks) 
library(rgl)

zz <- " longitude   latitude    depth   time    magnitude
363.218744  52.38412991 3.04532258  9.241866766 0.838097648
363.685488  53.85946071 0.25531946  0.578126952 0.854703175
319.732091  32.23166292 4.2952681   1.06539397  0.781598317
315.812817  36.98421462 3.27938697  2.295483245 0.769473951
119.498255  75.1280934  0.83253364  8.783517513 0.969860223
138.588107  169.629965  1.89874255  9.745653758 0.972609198
139.158872  169.6452912 3.67211991  15.03850063 0.881092856
139.247168  169.3478202 2.85800531  15.26106239 0.809570945
65.1234      147.3603323   2.84478705  16.82280753 0.87957715
65.55518    147.0377814 2.69845292  20.23640417 0.900546144
67.954278   148.6017748 3.41456545  21.21715903 0.968369227
66.968587   146.3484298 2.33379707  22.02860946 0.726200657
359.087074  42.8644051  2.85850316  23.70757068 0.861703105
357.875617  40.08052697 2.90391604  23.84198256 0.828534981
358.483231  41.26875168 1.98773836  26.06830006 0.885146048
56.391874   87.82219658 1.18588976  31.82632342 0.767755206
56.355593   87.68729915 1.37666954  33.07983844 0.731355859
168.510068  6.14462771  3.46656886  96.53093276 0.752218258
167.103837  6.28745699  3.55633497  97.11208455 0.718645583
183.30843   45.57593198 0.24446465  98.29165891 0.821997826
183.036565  51.9529565  0.52456176  98.81054412 0.838918878
276.323019  6.86787982  1.7776151   135.5215818 0.855317054
275.970301  7.91419119  3.21748128  135.6636361 0.718547414
275.727285  8.3243395   1.94806409  136.057032  0.769410166
276.514291  7.01085361  1.43467439  137.343003  0.735545797
330.040302  144.5877302 0.66277381  137.4281315 0.729634011
324.893213  146.4365159 0.88234134  137.6037721 0.967082637
291.11468   49.33367767 2.47072766  138.7895341 0.710016723
253.512945  201.0537249 2.4727616   141.0409587 0.847133697
253.542791  201.2490377 1.98635162  142.6391569 0.879832516
257.795509  216.9359492 0.78064425  145.1118214 0.72706484
90.28441    126.6897782 1.01474716  154.9277926 0.823621983
87.735243   125.8623578 3.00304403  155.2644112 0.815248951
260.385862  96.59499166 3.96772436  156.995524  0.85693796
259.739291  97.04823338 2.61008235  157.5583662 0.96764469
260.311284  96.68968926 4.22429753  157.6724828 0.821776192
260.229276  95.34657991 2.54473831  157.7533106 0.881316669
260.212422  187.7862368 4.31051589  157.9162288 0.728441067
260.231256  188.2226327 3.3502332   158.0074879 0.800794383
261.545309  191.0387749 3.93897914  158.0767802 0.743452552
272.998243  214.0871297 0.58707332  158.0785985 0.773464159
272.69374   214.5296828 1.42084384  159.6240869 0.813691498
274.126863  212.7891744 0.26609394  159.9565838 0.781095884
367.90561   15.41286762 2.95145208  160.5243755 0.76467979
367.864872  15.25877454 2.64007335  161.3173639 0.890119333
367.649781  19.91579374 3.97005021  161.9414598 0.939914671
148.118458  177.8767452 0.7438971   162.0120391 0.762554208
147.655051  176.9144773 1.3703557   162.030937  0.968876858
148.050091  176.9031162 1.30790609  162.3061767 0.923712841
147.649022  175.9846377 0.01212553  162.6304654 0.860915742
356.818492  67.21939123 1.76786772  164.2579604 0.990989558
355.622664  68.75829882 3.19659511  164.3880109 0.869464951
355.844738  68.88550929 4.36517304  164.7900766 0.76530697
108.71847   215.0266242 1.98752526  166.393087  0.957328986
109.014859  214.6488939 3.20229933  167.6177713 0.958599795
109.164433  215.2627269 1.93777216  167.8951075 0.793990658
108.131203  215.0862499 2.35653921  168.2758637 0.798095263
196.231361  70.49924437 1.97858279  168.4511479 0.952004021
197.737244  72.88240664 0.65482871  169.5036831 0.847759836
339.257305  120.888616  1.34736322  169.6356306 0.715570745
340.781081  117.0187847 2.51971051  169.6790878 0.933405511
338.606336  121.5834442 1.80270464  170.3047686 0.828001716
338.041063  120.7029464 1.73078138  170.3913998 0.783988564
338.99619   35.10033383 2.97337925  170.4444717 0.800052404
337.200856  33.51899806 1.09814538  170.613345  0.810081986
336.985755  34.21892728 0.24128014  170.8469088 0.938437723
28.837549   219.7118007 0.88414209  176.6881568 0.882990987
417.469016  60.71793319 3.13009015  178.4274218 0.815114554
420.412299  59.95963986 3.77420788  178.5608842 0.709170932
417.835488  61.34888422 3.412701    179.0698474 0.957269676
266.303184  148.8628134 4.14332284  179.975606  0.881062956
266.264066  149.8292293 4.17012164  180.080254  0.954352122
290.584215  94.18249675 2.09946968  180.3264641 0.851384874
211.043848  164.8354141 0.38617853  180.4780309 0.702373251
394.665914  165.8533229 2.94896243  180.6401804 0.975721953
396.948844  164.3185096 2.95721695  181.189757  0.823211439
394.152446  166.1863908 2.85283338  181.4459117 0.919709724
355.861861  151.482828  1.79907458  181.5735067 0.745331584
355.604564  151.2328554 1.42786694  182.0786425 0.944569983
355.977914  151.7258945 2.60632787  182.2734472 0.97288076
55.193224   70.13457335 1.95421619  182.4681836 0.933382738
55.121245   70.26190069 1.7134954   182.7419885 0.927136941
466.029024  43.60374817 3.02745392  183.048289  0.821559257
466.713681  43.04071215 3.34442703  183.0495735 0.910570324
361.343579  38.84946077 4.09011265  183.1337707 0.869638281
107.494724  224.061154  0.64225969  184.3347164 0.875420906
109.403587  222.7703059 1.02129489  184.5405119 0.80347438
108.536435  221.8397702 1.32737534  185.2302114 0.911269972
108.966261  223.5009692 0.64308213  185.599194  0.748259104
319.258848  64.01371204 4.34421863  185.7281496 0.70045418
321.53673   62.05801435 4.02840039  185.9683567 0.868600256
320.136879  62.7561261  2.02796518  186.1012058 0.863945237
426.321536  79.57156783 1.80279196  186.1655452 0.865248329
425.254692  80.84699203 0.53246863  186.351376  0.908361028
424.301716  81.17660654 0.05602379  186.4132518 0.978248228
424.827296  81.04033515 1.33209986  186.4734765 0.983263093
119.807247  187.958236  3.33382966  186.7652667 0.970978999
120.210092  189.4849366 2.84489764  187.4376306 0.881231886
120.205569  189.1727403 4.18836325  187.4622579 0.977862927
247.808605  146.4244449 0.83110586  188.5004803 0.735083878
249.294833  144.1907226 1.22128356  188.962155  0.802240452
250.204516  146.4749518 1.7944476   189.3051254 0.711019133
334.654137  11.15448488 3.46662664  189.8126837 0.861510955
333.659862  11.64954098 2.1332473   190.1436423 0.92116864
333.427405  10.95652538 2.75561873  191.3076606 0.754366751
9.482643    85.46858555 1.15438219  191.3293578 0.827873998
9.378643    85.45893223 2.23185957  191.4270999 0.729450387
9.333447    84.55839727 1.36554428  192.1174705 0.79096647
439.424485  237.4993259 0.59464444  192.4302016 0.934036784
438.638525  238.1309002 0.47642958  192.835047  0.861088125
439.977322  237.950318  1.51017279  193.0183715 0.93072143
438.636146  237.9862702 1.3178479   193.2070887 0.891683094
330.410302  62.51368701 2.99926551  193.8201965 0.967944474
329.790522  62.17494892 2.86913557  194.8657309 0.718207116
329.460954  62.13150555 3.26717085  195.1562923 0.752731464
474.563651  69.88329745 4.03771076  195.1657446 0.824899558
285.627382  72.98472026 1.95847784  196.1019459 0.921937552
286.820235  72.90395374 1.26502975  196.2494223 0.967886438
236.897161  134.567653  1.15564329  196.4295919 0.707755992
234.806852  133.6690276 0.210115    196.5220457 0.741268909
233.53859   132.9487054 2.28750244  196.5612547 0.82722623
277.005222  47.68301903 4.36155011  196.6331458 0.929370035
277.336998  44.11638043 3.06712344  197.5726448 0.857331156
277.800836  43.03107739 3.70528439  197.5837294 0.926344408
276.793356  44.16753185 3.20462984  197.9290181 0.750939533
165.177531  248.6222625 0.3374573   197.9766126 0.901818102
164.463612  249.5902955 1.32364566  198.1684532 0.885595823
161.499292  251.3536746 1.75217305  198.2615641 0.702053827
440.799878  176.1199082 3.19808685  198.3114499 0.922067657
441.798136  176.5611422 2.98843581  198.9449421 0.997520643
441.027104  176.7791779 3.08648705  199.572324  0.738454407"

y1 <- read.table(text=zz, header = TRUE)

zzz <- " 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"
y2 <- read.table(text=zzz, header = TRUE)

evpts1 <-do.call(expand.grid,  lapply(y1,quantile, prob=c(.1,.15,.2,.25,.3,.35,.4,.45,.5,.55,.6,.65,.7,.75,.8, .85,.9,.95)))
hat1 <- kde(y1, eval.points= evpts1)

evpts2 <-do.call(expand.grid,  lapply(y2,quantile, prob=c(.1,.15,.2,.25,.3,.35,.4,.45,.5,.55,.6,.65,.7,.75,.8, .85,.9,.95)))
hat2 <- kde(y2, eval.points= evpts2)

str(hat1)
str(hat2)

integrand1 <- function(p){hat1(p,y1)+hat2(p,y2)}

Vintegrand <- Vectorize(integrand1)

integrate( Vintegrand, lower = 0, upper = t)
 

Dason

Ambassador to the humans
#4
Just to confirm - you're looking for the integral of the pointwise sum of the kernel density estimates. You're not looking for the kernel density estimate of the sum of the two random variables.

Note that the integral of a sum is the sum of the integrals so you could just integrate each of the densities and then add those together. Assuming this is really what you want.