0

I fit a multivariate distribution in R with below data. I need help in understanding the output.

      library(mvtnorm)
      h2
           a        b       c        d        e
      0.3837  1.11516  1.9256 -0.50875  1.88859
      0.8426 -0.84579 -0.7873  0.02639  0.33517
      0.1542 -0.09256 -0.1094 -1.38610  0.54481
     -1.6815  3.56519  3.1228  1.48618  1.53065
     -0.5342  0.71441  0.6989 -1.55011 -0.41603
      1.0721  4.52062  3.8051 -0.16023 -0.01286
     -0.3047  1.27455  0.7179  0.27030 -0.34911
      2.9078  1.61245  1.2929  0.21195  1.32390
     -0.3047 -0.51061 -0.4536  1.54295  1.08431
     -0.7636  0.96761  1.2947  0.73289 -0.59614

      mu = colMeans(h2)
      C1 = cov(h2)

    # Below function will compute probability for each row
      F2 = function(b)
      {
         pmvnorm(lower=-Inf, upper=b, mean=mu, sigma=C1)
      }
     # Below function will compute density for each row
      F3 = function(b)
      {
        dmvnorm(x=b, mean=mu, sigma=C1, log=FALSE)
      }


      proba=apply(h2,1, F2)
      dens = apply(h2,1, F3)
      cbind(h2,proba,dens)

          a        b       c        d        e    proba      dens
      0.3837  1.11516  1.9256 -0.50875  1.88859 0.072544 0.0007242
      0.8426 -0.84579 -0.7873  0.02639  0.33517 0.019907 0.0065792
      0.1542 -0.09256 -0.1094 -1.38610  0.54481 0.005645 0.0031873
     -1.6815  3.56519  3.1228  1.48618  1.53065 0.047796 0.0006827
     -0.5342  0.71441  0.6989 -1.55011 -0.41603 0.001885 0.0030992
      1.0721  4.52062  3.8051 -0.16023 -0.01286 0.112288 0.0013975
     -0.3047  1.27455  0.7179  0.27030 -0.34911 0.022331 0.0062871
      2.9078  1.61245  1.2929  0.21195  1.32390 0.271875 0.0012214
     -0.3047 -0.51061 -0.4536  1.54295  1.08431 0.032490 0.0026146
      -0.7636  0.96761  1.2947  0.73289 -0.59614 0.013268 0.0007913

My question, why the values of probability and density are very low for each row.

My interpretation is that, each of the rows does not belong to a multivariate normal distribution with mu and C1. Is that correct?

Tim
  • 108,699
  • 20
  • 212
  • 390
Chandra
  • 143
  • 7

2 Answers2

1

The values of cumulative probability and density are as low as they are due to the effects of being in 5 dimensions. Your expectation as to what is reasonable is perhaps based to a large extent on 1 dimensional normal.

I checked the results for probability and density using MATLAB, and got exact agreement, so I think your results are correct.

I will show you just the cumulative probability results using only some of the columns, and you can start to get a feel. These correspond to your proba output column, but only using some of the columns of input data. I am not showing you every possibility.

Just the 1st column (1D):

   0.565637288338181
   0.702816002597073
   0.492667067696748
   0.068452074479719
   0.284582761085686
   0.763058850958907
   0.349888595566043
   0.985563411940202
   0.349888595566043
   0.225766862907170

Just the 2nd column (1D):

   0.472501841180498
   0.110154392554774
   0.217285938950877
   0.915626527024584
   0.380039621684968
   0.973799841208276
   0.509987926467789
   0.588761043921356
   0.151976916659519
   0.438008974406623

Just the 1st and 2nd columns (2D):

   0.264841990150589
   0.076405351627506
   0.105236178888105
   0.062355658370148
   0.106157429165297
   0.742779509346500
   0.176154283011264
   0.580041734712487
   0.051831899178908
   0.097057907212853

Just the 1st, 2nd, and 3rd columns (3D):

   0.264617423149077
   0.057300913569709
   0.083949731409476
   0.060366682143580
   0.093198795122322
   0.731508503440227
   0.126580702128226
   0.510771711413619
   0.040244089729650
   0.094487914775201

Just the 1st, 2nd, 3rd, and 4th columns (4D):

   0.073559770009587
   0.031623966714283
   0.006727668874189
   0.050106397108077
   0.003951924232790
   0.277948705795869
   0.067743297644953
   0.306330906498960
   0.037353799984810
   0.064902776076491

You can of course look at other combinations of the columns.

Going back to the original problem in 5D, I calculated the Mahalanobis Distances of the input data points relative to the fitted mean and covariance, and converted these to the corresponding probability levels from the Chisquare distribution which they would follow if truly multivariate normal with that mean and covariance. The match is not great, although the data set is small. If these were sorted and the sample large enough, and the match were good, then it should be close to a straight line. This is essentially getting at the distribution of points ellipsoidally about the mean relative to the fitted covariance - that has taken the effect of number of dimensions out of the picture, and is perhaps what you really (should) care about.

    0.7359
    0.1573
    0.3759
    0.7460
    0.3844
    0.6015
    0.1701
    0.6325
    0.4349
    0.7203

Full disclosure: I didn't do these Mahalanobis Distance calculations very carefully, so could have made a mistake.

Mark L. Stone
  • 12,546
  • 1
  • 31
  • 51
0

Any real number is within support of normal distribution, this applies also to multivariate normal distribution. So any real number "belongs" to multivariate normal with any parameters, some values are just observed with smaller probabilities.

The probabilities that you see are low because this is a multivariate distribution. Try to visualize yourself univariate and bivariate normal distributions. In the first case you observe probabilities of points on $x$-axis, while in the second case probabilities of points in two dimensional space. Obviously multivariate distribution relates to greater total space, so individual probabilities for each point (or interval) are proportionally smaller than in univariate case. Imagine that you are interested in marginal distribution of one of the variables in bivariate normal distribution, to obtain it you would sum up (take integral) all the points from one of the dimensions,

$$ P(X = x) = \sum_y P(X = x, Y = y) $$

so from this alone you can imagine that the individual probabilities in bivariate case were smaller than in univariate case.

You can frame your question differently and ask if your sample fits the 95% (or other) highest density interval, i.e. if they fall below some predefined threshold. You can be also interested how likely they are under such model and if another model does fit your data better. Those however are different questions that need different approach than looking at the raw probabilities alone.

Tim
  • 108,699
  • 20
  • 212
  • 390
  • Tim, I am interested in the answers for the questions that you mentioned. But, I was baffled by the results in first step and wanted to understand what is going on in multivariate case. I will post a separate question after doing some study in the direction, Mark pointed above. Thank you for the reply. – Chandra Dec 19 '15 at 16:52
  • @Chandra I understand and so my answer is: nothing is going on. If you are defining "belonging" to distribution in terms of quantiles, than you must ask yourself what is the threshold that marks boundary between fitting and not fitting. You could also compare fit of different models to your data. But the probabilities alone does not provide any such information. So I am saying that either your question is ill-posed or you use inappropriate methods for your problem. – Tim Dec 19 '15 at 17:36
  • Objective of my exercise is to do anomaly detection. Any anomaly should have very low probability of belonging to fitted multivariate distribution. However, I was surprised by the low probability levels in multivariate distribution, during the training. I guess any one working with multivariate distribution would have similar question. I did refer to some books on multivariate analysis and could not find a satisfactory answer. So, I don't think this Q is ill-posed. However, I certainly need to find a better method. Please let me know, if you know a good reference. Thanks. – Chandra Dec 19 '15 at 22:22
  • @Chandra it depends how you define "low" probabilities, that is why I was saying about arbitrary threshold that has to be chosen. Try to imagine multivariate distribution -- probability of observing any single combination of points (in continuous case interval) would be low because it is a combination of circumstances. – Tim Dec 19 '15 at 22:54
  • @Chandra see my edit – Tim Dec 19 '15 at 23:14