R help - Oct 2004 - strange results with dmvnorm

I am experiencing strange results using dmvnorm. I define a scaled distance
matrix from the coordinates bellow and then calculate a covariance matrix using
a spherical correlation function.  Then with certain combinations of
range and sill parameters dmvnorm is returning values greater than 1.  Surely
the results of dmvnorm should be in the interval 0:1 (or do I just nead a
holiday?).  In addition cov.spatial and sphercov are giving different results. 
Any clues as to what is happening here might just save my sanity.

Many Thanks
David

x <-
c(22038,50284,56009,46979,21417,10149,31305,34473,42883,55084,51097,67990,55710,20256,31972,46423,20321,70886
,18758,43727,23834,13612,70217,38389,7980,41711,19462,35191,13924,14460,52657,68895,57616,48606,53874,42444
,17307,4514,38631,20517,42091,51844,42137,57977,22168,19906,31072,56385,51056,26229,20438,34406,28462,17218
,47406,84182,29189,46444,21067,37383,24395,38658,38617,70943,65331,63586,50953,39076,51230,64931,13615,68210
,23633,2455,56004,58891,59927,60778,70023,2075,7118,47964,40564,37949,3033,28401,41224,49202,50643,57889
,48008,11062,73256,38052,64000,60843,39599,45633,70129,25122,54528,42006,54445,38568,28561,37529,43752,67653
,11958,54801,52046,52773,42120,18872,63421,72778,57048,22565,52248,12549,21301,14375,961,34307,40936,44523
,45671,2651,52247,60016,50844,24076,55843,19700,12395,30860,21322,33055,43677,33405,66787,36370,28667,31794
,51464,44090,26633,25823,34199,52506,67256,72894,44651,29108,69867,37312,24407,59456,63329,61886,30339,69144
,13877,47091,13980,24024,68432,1007,6432,13923,31904,41093,18176,22861,25122,45251,35276,8508,19127,73074
,29659,41765,27467,49725,39012,26689,8794,11480,33809,37180,71153,44432,18397,66592,49294,31554,8914,40984
,37200,43233,47493,33728,58487,50731,42288,52529,33187,28559,21471,45940,40803,24063,35777,31187,22374,42185
,37019,33101,17435,21451,64085,48624,0,72571,18484,51457,44206,49327,37643,71719,45009,47327,58397,38843
,27057,28251,34786,32818,26227,16233,51900,42102,70651,49805,33860,19510,6536,67032,27022,20490,74443,56503
,23201,8504)

y<-c(19656,24181,21863,51776,59164,40430,28033,59460,17720,28689,39794,44737,33744,60868,58488,41163,55696,28305
,39258,12051,48751,51806,31246,9147,24850,43030,62512,39028,37460,60933,56162,15139,36460,19472,62064,16044
,52145,42288,15933,52472,4508,25879,28955,39947,54684,48753,35232,43390,33528,19161,15281,8860,37159,41399
,29888,43382,31347,19873,42867,38766,46147,12663,26343,18649,15733,41152,31031,14679,49137,48205,46452,31114
,52894,22521,29526,29937,21176,12439,47782,34605,42022,21109,14249,56541,31478,15300,45095,45090,26764,34286
,14803,24848,13272,29076,34605,31612,7278,28536,454,38202,21544,13010,25735,48565,45640,5288,32648,33012
,35712,35420,21666,65959,21632,28182,28024,19770,23977,39348,19231,27179,45174,14399,30525,25703,31037,19954
,23285,47695,21027,13917,43449,34856,38771,14150,47183,41164,21451,34398,47306,15420,17779,51974,42830,51551
,61578,8379,41219,48378,7699,40958,25624,23373,7979,53196,14330,56537,58067,14237,26000,33423,39967,28781
,27549,33739,33114,41793,20995,42079,34549,41766,39550,1595,35390,42489,22246,7580,20631,29990,23001,21306
,20456,39752,43611,42768,0,35443,48950,45480,48325,10746,42252,29786,45882,45156,48526,19182,51870,37749
,18564,52335,59215,38424,27250,12685,10016,37640,10945,53196,49621,37746,15965,15008,7207,14364,27606,35459
,36399,44440,31733,25446,21065,65535,37627,26158,64824,10465,58176,36662,19769,9080,45041,16080,49648,59638
,23339,40697,47107,8536,57035,55090,44414,1321,12861,21108,32654,27068,38365,22255,31550,11789,45404,53969
,13509,36350)

Dist <- sqrt(outer(x,x, "-")^2 + outer(y,y, "-")^2)
Dist <- Dist/max(Dist)

library(spatial)
Cov1 <- sphercov(Dist, 0.8, alpha=0, se=sqrt(2))
Cov2 <- sphercov(Dist, 0.6, alpha=0, se=sqrt(0.55))
library(geoR)
Cov1b <- cov.spatial(Dist, cov.model= "spherical", cov.pars=c(2,
0.8))
Cov2b <- cov.spatial(Dist, cov.model= "spherical", cov.pars=c(0.55,
0.6))

library(mvtnorm)
dmvnorm(rep(0, nrow(Cov)), sigma=Cov1)
[1] 37949.22
dmvnorm(rep(0, nrow(Cov)), sigma=Cov1b)
[1] 2.920084e-05
dmvnorm(rep(0, nrow(Cov)), sigma=Cov2)
[1] 4.722488e+60
dmvnorm(rep(0, nrow(Cov)), sigma=Cov2b)
[1] 1.085243e+51


----------------------------------------------------------------
Concerns about content should be sent to abuse at salford.ac.uk

On Sun, 31 Oct 2004 16:04:53 +0000, "David Pleydell"
<D.R.J.Pleydell at pgr.salford.ac.uk> wrote:
>I am experiencing strange results using dmvnorm. I define a scaled distance
>matrix from the coordinates bellow and then calculate a covariance matrix
using
>a spherical correlation function.  Then with certain combinations of
>range and sill parameters dmvnorm is returning values greater than 1. 
Surely
>the results of dmvnorm should be in the interval 0:1 (or do I just nead a
>holiday?).  In addition cov.spatial and sphercov are giving different
results.
>Any clues as to what is happening here might just save my sanity.
Why would you expect a density to be between 0 and 1?  

Duncan Murdoch

"David Pleydell" <D.R.J.Pleydell at pgr.salford.ac.uk>
writes:
> Surely
> the results of dmvnorm should be in the interval 0:1 (or do I just nead a
> holiday?).  
You need a holiday. Consider

  curve(dnorm(x,sd=.25),from=-1,to=1)

The *integral* has to be 1 though.

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907

David Pleydell <D.R.J.Pleydell <at> pgr.salford.ac.uk> writes:

: I am experiencing strange results using dmvnorm. I define a scaled distance
: matrix from the coordinates bellow and then calculate a covariance matrix 
using
: a spherical correlation function.  Then with certain combinations of
: range and sill parameters dmvnorm is returning values greater than 1.  Surely
: the results of dmvnorm should be in the interval 0:1 (or do I just nead a
: holiday?).  In addition cov.spatial and sphercov are giving different 
results. 
: Any clues as to what is happening here might just save my sanity.

Probabilities must not exceed 1 but densities can.  In fact,
the normal density will approach a delta function as the standard
deviation approaches 0.    For example, in the one dimensional
case note how the density of the standard normal at x=0 increases 
as the standard deviation decreases from 1 to .2:

R>     dnorm(0,0,5:1/5)
[1] 0.3989423 0.4986779 0.6649038 0.9973557 1.9947114

or try plotting the densities to get a visual idea:

R> matplot( sapply(5:1/5, dnorm, x = (-50):50/10, mean = 0), type =
"l")