Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

My problem is related to statistical testing with weighted data. I have a data sample and I assume that a part of the data follows a specific distribution (Pareto). So I fit my data to a distribution by using EstimatedDistribution and in turn I obtain a p-value related to this fit according to different statistical tests (Cramer-Von Mises, Kolmogorov-Smirnov and possibly others).

Everything goes well as long as I do not include weights (I really know it works, since I have replicated my results in another software package...).

When including weights by using the function WeightedData I run into problems. EstimatedDistribution still works well on WeightedData (again results are replicated in another package), but the DistributionFitTest does not work any longer.

When trying to apply DistributionFitTest on WeightedData I get an error message of the following form:

"rctnln: The argument [my weighted data object] at position 1 should be a rectangular array of real numbers with length greater than the dimension of the array."

I got the intuition behind these error message and intensely searched for an appropriate work-around or alternative Function, but could not find any suitable tools.

Could anybody help me out on this?

share|improve this question
1  
A little bit of code to reproduce this: data = {8, 3, 5, 4, 9, 0, 4, 2, 2, 3}; w = {0.15, 0.09, 0.12, 0.10, 0.16, 0., 0.11, 0.08, 0.08, 0.09}; wd = WeightedData[data, w]; DistributionFitTest[wd]. Seems like DistributionFitTest wasn't updated when WeightedData was introduced. However, DistributionFitTest[EmpiricalDistribution[w -> data]] doesn't work either and EmpiricalDistribution was introduced a version ago. – Sjoerd C. de Vries Jul 11 '13 at 20:51
    
I also tried the variant you suggest (i.e. using DistributionFitTest[EmpiricalDistribution[w -> data]]) when I was still using version 8. Basically it was the same issue as far as I remember. Another variant I tried was to integrate a RandomVariate to draw form the empirical distribution and test this artificial data against the distribution. This worked out well, but I was unsure whether this really represents what I am trying to do... – Joko Jul 11 '13 at 22:40
    
You may consider reporting this at support@wolfram.com. – Sjoerd C. de Vries Jul 11 '13 at 23:03
    
I got the following comment from Wolfram: "Unfortunately, DistributionFitTest does not currently support WeightedData objects. Our developers are aware of user interest in this feature, and are working to include it in a future version of Mathematica. On your behalf, I've filed your comments with our development team for their review as they move forward. We are always interested in improving Mathematica, and I want to thankyou ...." So there will be no solution to my problem on this front - has anyone a suggestion for an appropriate work-around? – Joko Jul 20 '13 at 10:50
1  
How did the weights come about? Did you sample some parts of more intensely than others? Are some observations estimated more precisely than others? How one might appropriately use weights should start with a definition of what the weights are. – Jim Baldwin Dec 16 '15 at 4:40

DistributionFitTest[data] test the hypothesis that the data came from normal distribution.

This is typically done by testing if the empirical CDF is close enough to the CDF, expected under the null hypothesis.

The first visual test is to inspect QQ-plot, implemented in Wolfram Language as QuantilePlot.

Consider synthetic data:

refD = NormalDistribution[1, 2];
baseD = LogisticDistribution[1, 2];
data = RandomVariate[baseD, 10^3];
weights = PDF[refD][data]/PDF[baseD][data];
wd = WeightedData[data, weights];

Visualize QQ-plot against normal:

qd = QuantilePlot[wd, NormalDistribution[1, 2]]

parametric plot of empirical quantile vs.

The fit is pretty good here, and many tests assess the goodness of fit by computing some distance between these two lines and assessing whether the so obtained number is typical under the null.

For instance, we can compute Anderson-Darling test statistic as follows:

ed = EmpiricalDistribution[wd];

ad$test = 
 NExpectation[(CDF[ed, x] - CDF[refD, x])^2/(CDF[refD, 
      x] SurvivalFunction[refD, x]), x~Distributed~refD]

(* Out[13]= 0.000747316 *)

or Cramer-von-Mises test statistic:

cvm$test = 
 NExpectation[(CDF[ed, x] - CDF[refD, x])^2, x~Distributed~refD]

(* Out[14]= 0.000128502 *)

or Kolmogorov-Smirnov test statistic:

ks$test = MaxValue[Abs[CDF[ed, x] - CDF[refD, x]], x]

(* Out[15]= 0.0203543 *)

Since the sample size is rather large for this dataset, we might be OK with using the asymptotic distributions to compute the p-value. It is less clear what to do with small datasets. I would expect the test statistics' distribution to depends on the actual weight vector.

The safest route would be to use Monte-Carlo method, which computes test statistics on samples from the null to determine if the test statistics computed on data is typical. The problem is that sampling does not give you WeightedData, hence we might be comparing apples to oranges.

share|improve this answer
    
@whuber If you have comments, or references on the matter, I would greatly appreciate it if you shared them. Thanks! – Sasha Dec 18 '15 at 15:08

The problem is that EmpiricalDistribution[] does not generate a PDF function, only a CDF. A work around could be to create an un-weighted data set out of your data. It's sub-optimal to say the least, and you may run out of memory on some data sets before getting your answer, but it works.

Flatten[ConstantArray @@ # & /@ Transpose[{data, Ceiling[100*w]}]]

Note that I used 100 as the multiplier since your weights were all in hundredths. Also, weight 0. for InputData 0 makes no sense. You have an element which has no weight, e.g. is not possible in the data. Just remove it. Or set frequency to something very small.

share|improve this answer
1  
I am really not sure why you believe PDF is of any relevance here. Moreover, PDF is defined: PDF[EmpiricalDistribution[{0.8, 0.2} -> {1, 2}], {0, 1, 2, 3}] outputs {0, 0.8, 0.2, 0}. – Sasha Dec 15 '15 at 20:30

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.