I'm trying to build a bivariate copula-based model of income and wealth in Italy and I'm having trouble handling weighted data. I have access to micro data, a survey of about 10,000 households that includes the corresponding sample weights.

When calculating basic statistics (like mean and median) and even when performing linear regressions it is pretty easy to account for weights, besides there are useful packages for that (e. g. survey). But what do I do when I want to fit a parametric model of the distribution to weighted data? Or to estimate its kernel density?

I have a few ideas, but they seem to be pretty crude. For one, I could inflate my sample to the size of the universe. That is, I could multiply all weights by 100 (which would turn them into integers) and then create a vector that repeats each value of income and wealth a given number of times. But that would lead to a very large sample (which I believe still wouldn't be a perfect representation of the population) and will certainly put some extra strain on my computer.

I could also just round the weights off instead of multiplying them by 100, but this would still make the sample noticeably bigger and will inevitably skew the real proportions.

Another approach I came up with would be to normalize the weights (so that they sum up to one) and then randomly sample with repetitions from my initial sample with the corresponding vector of probability weights. R doesn't allow to draw the samples that are larger in size than the one that they are being drawn from. But I think that drawing the sample of the same size as the initial one will lead to some loss of information about the observed proportions. So I could draw the samples of the initial size as described above several times (how would I know how many is though?) and then combine them into one sample. And again, I will have a larger sample with some of the information lost along the way.

So I was wondering if there is a better way to handle weighted data. In some cases I think I could technically introduce the weights into the formula for computing the maximum likelihood for fitting a particular model, although I certainly wouldn't like to code that from the ground up. I will have to fit a lot of models as part of my project, both univariate (e. g. Singh-Mandala) for income and wealth and bivariate for copulas. I don't think the built in functions in any of the copula-related packages that I'm aware of allow one to account for weights. So any advice would help!


Actually estimating kernel density of weighted data is pretty straightforward. Recall that standard kernel density estimator is

$$ \hat{f_h}(x) = \frac{1}{n} \sum_{i=1}^n K_h(x-x_i) $$

where $K_h(x) = K(x/h)/h$ for kernel $K$ and bandwidth $h$, then it may be thought as weighted kernel density with all weights equal to $w_i = 1/n$,

$$ \hat{f_h}(x) = \sum_{i=1}^n w_i K_h(x-x_i) $$

this easily translates to the unequal weights scenario, where $\sum_{i=1}^n w_i = 1$. Moreover you don't have to look far for software estimating it as it is often available, e.g. R's density has weights parameter for such case.

As about fitting traditional distributions to weighted data, you can introduce weights in the likelihood function. Notice that if you re-sampled your data, so that some points would get duplicated, then you could obtain sample such as $x_1,x_1,x_1,x_2,x_2,\dots,x_k,x_k,x_k,x_k$, then the likelihood function becomes

$$ L(\theta\mid x_1,\dots,x_k) = f(x_1;\theta)\times f(x_1;\theta)\times f(x_1;\theta)\times f(x_2;\theta)\times f(x_2;\theta)\times \dots \\ \times f(x_k;\theta)\times f(x_k;\theta)\times f(x_k;\theta)\times f(x_k;\theta) = \\ f(x_1;\theta)^{3}\times f(x_2;\theta)^{2}\times\dots\times f(x_k;\theta)^{4} = \prod_i f(x_i;\theta)^{n_i} $$

where $f$ are density functions evaluated at $x_i$ points and parametrized by $\theta$. So if weights are non-negative and proportional to numbers of samples, then you can easily, by analogy, introduce weighted likelihood function of the same form

$$ L^w(\theta\mid x_1,\dots,x_k) = \prod_i f(x_i;\theta)^{w_i} $$

that can be maximized to find the optimal value of $\theta$. To lower the risk of numerical precision issues, it'd be better to work with it's logarithm,

$$ \log L^w(\theta\mid x_1,\dots,x_k) = \sum_i w_i \log f(x_i;\theta) $$

So you can use if for fitting the distribution by maximizing the weighted likelihood function.

  • $\begingroup$ As well as I understood the original question is already for weighted data. But how about this case: 1. we have some raw time-series data, 2. we want to fit a normal distribution to it (or any other), 3. we also want to use let's say exponentially declining weights (EWMA). Do we obtain weights for individual observation and then feed raw data and weights to the likelihood function as you provided to obtain parameters for weighted data? Can a likelihood function be designed so that to determine the weights (smoothing parameter in EWMA) and normal distribution parameters at the same time? $\endgroup$ – EmptyHead Jul 21 '17 at 19:27
  • $\begingroup$ Shall I post this as a separate question? $\endgroup$ – EmptyHead Jul 21 '17 at 19:30
  • $\begingroup$ @EmptyHead weights as in example above are something that is known a priori. If you'd try to "estimate" best weights and parameters of normal distribution, then you would need to estimate at least $N+2$ parameters with $N$ data points, so this won't be identifiable. $\endgroup$ – Tim Jul 21 '17 at 19:31

I'm thinking through a similar problem right now - my initial thought is to generate a histogram using the weights (eg. numpy's histogram takes a weights vector) and then fit a parametric distribution to the resulting (normalized) frequency vector. Curve fitting is pretty low headache to implement. I'm happy to go into more detail if it would be helpful.

  • $\begingroup$ The histograms you reference are interesting but it seems they are one-dimensional. Is there an extension of the method to two dimensions so that it will fit in with the OPs situation? $\endgroup$ – Michael R. Chernick Jan 6 '17 at 22:28
  • $\begingroup$ Yup! The folks at numpy have thought of everything - histogram2d and histogramdd are the 2d and multidimensional equivalents, respectively. $\endgroup$ – RedPanda Jan 7 '17 at 5:27
  • $\begingroup$ I feel bad for you. I don.t know how to help you and Tim's nice answer probably won't I hope someone here can. This seems to be important to you and you have done a lot of work on your own already. $\endgroup$ – Michael R. Chernick Jan 7 '17 at 12:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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