# Multiple linear regression with weights
**URL:** https://discourse.julialang.org/t/multiple-linear-regression-with-weights/7449
**Category:** General Usage
**Tags:** question, stats
**Created:** [December 1, 2017, 11:38pm UTC](https://discourse.julialang.org/t/multiple-linear-regression-with-weights/7449 "2017-12-01T23:38:34Z")
**Posts on this page:** 11
**Page:** 1
### Author:  [@yakir12](https://discourse.julialang.org/u/yakir12)
#### Post date: [December 1, 2017, 11:38pm UTC](https://discourse.julialang.org/t/multiple-linear-regression-with-weights/7449/1 "2017-12-01T23:38:34Z")
**Background:** I’m trying to fit a 2D symmetric Gaussian to a an image of a point source. So I have an image of a light source, it’s Gaussian in shape, and I want to find it’s exact center _and_ get a measure of how good the fit is. I plan to `log` the intensities and fit a symmetric paraboloid to the data. To compensate for the over represented tails of the Gaussian, I want to weight the data with itself (so low intensities will have low weights and vice versa). Here’s what it looks like (artificially enlarged):
**Question:** In order to fit a linear model to weighted data I assume I need to use `MultivariateStats.llsq`, but I can’t really see how. What is _the_ way to fit multiple linear regression to weighted data in Julia?
I found tons of functions in `StatsBase` that act on the subtypes of `RegressionModel`, but not one function that produces any instances of its subtypes…
Thanks!
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [December 2, 2017, 7:17am UTC](https://discourse.julialang.org/t/multiple-linear-regression-with-weights/7449/2 "2017-12-02T07:17:20Z")
How is this a regression? My understanding is that you are just fitting a multivariate normal.
If you want to deal with non-normal tails, the right approach would be something that allows for that, eg a [multivariate t](https://en.wikipedia.org/wiki/Multivariate_t-distribution). `Distributions.jl` has it implemented, so you have the log likelihood, but AFAICT it has no `fit_mle` method for it; you can do the MLE yourself with a generic optimization package (eg `Optim.jl`). You may get away with fixing the \nu at some value, eg 4 or 5.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [December 2, 2017, 9:17am UTC](https://discourse.julialang.org/t/multiple-linear-regression-with-weights/7449/3 "2017-12-02T09:17:34Z")
I found my copy of the Kotz-Nadarajah book about multivariate t at home, and apparently for a known \nu, there is a closed form for the bivariate case. See
```nohighlight
@article{tiku1992estimation,
title={Estimation and hypothesis testing for a new family of bivariate nonnormal distributions},
author={Tiku, ML and Kambo, NS},
journal={Communications in statistics-theory and methods},
volume={21},
number={6},
pages={1683--1705},
year={1992},
publisher={Taylor \& Francis}
}
```
Apparently, if I am reading it correctly, the sufficient statistics are independent of \nu. So you can very quickly try various degrees of freedom and just pick the best fit.
---
### Author:  [@yakir12](https://discourse.julialang.org/u/yakir12)
#### Post date: [December 2, 2017, 10:52am UTC](https://discourse.julialang.org/t/multiple-linear-regression-with-weights/7449/4 "2017-12-02T10:52:54Z")
Thanks @Tamas_Papp!
I’m sorry to admit that your help was wee bit over my head and lead me down a different path:
I can fit a normal distribution to the coordinates of the pixels and use the intensities as the weights. This works perfectly when there is no noise. It breaks down completely if there is noise and/or if the image is much larger than the scale (say the FWHM or sigma) of the distribution.
Here’s a quick and (really) dirty example code:
```julia
using Plots
gauss(x,y,mux,muy,sig) = exp(-((x-mux)^2 + (y-muy)^2)/(2*sig^2))
sig = sqrt(3)
mux = 4.5
muy = 6.
f(x,y) = gauss(x,y,mux,muy,sig) + rand()/10
x = Float64.(collect(1:10))
y = Float64.(collect(1:10))
X = x .+ 0*y'
Y = 0*x .+ y'
z = f.(X,Y)
h1 = heatmap(x,y,z)
using Distributions
d = fit_mle(IsoNormal, hcat(vec(X), vec(Y))', vec(z))
z2 = [pdf(d, [i, j]) for (i,j) in zip(X, Y)]
h2 = heatmap(x,y,z2)
z *= sum(z2)/sum(z)
plot(h1,h2)
norm(z - z2)
```
This isn’t the best way.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [December 2, 2017, 11:14am UTC](https://discourse.julialang.org/t/multiple-linear-regression-with-weights/7449/5 "2017-12-02T11:14:11Z")
I don’t think you are doing this right.
First, `gauss` is \propto a pdf, why are you adding random numbers to it? This does not make much sense. To generate your random variates, you would draw from a distribution, _then_ bin them like the image. What you are doing is closer to some non-standard overdispersed normal.
For binned data, your log likelihood function can be approximated as
\sum\_i n\_i \ell(x\_i \mid \theta)
where \ell(x \mid \theta) is the log likelihood for the distribution, x\_i is the center of each bin, n\_i the count, and eg \theta = (\mu, \Sigma) are the parameters.
For normal, there is a closed form \theta that maximizes this (deriving this is a good exercise, you just get weighted means and variances, which you can calculate with `mean` and `cov` from `StatsBase`, with n\_i as the weights — see `WeightVec`).
For multivariate-t, you need the method I mentioned. But now I am not sure you actually need it; you just get the wrong estimates because you are not doing the right thing and simulating something other than a multivariate normal.
---
### Author:  [@yakir12](https://discourse.julialang.org/u/yakir12)
#### Post date: [December 2, 2017, 11:37am UTC](https://discourse.julialang.org/t/multiple-linear-regression-with-weights/7449/6 "2017-12-02T11:37:44Z")
You are of course correct 🙂
Let me be clear and try a more correct way of generating some test data:
```julia
using Images, ImageView, Distributions
d = MvNormal([1,2], 2)
img = centered(zeros(13,13))
for rep in 1:1000
i, j = round.(Int, rand(d))
if checkbounds(Bool, img, i, j)
img[i, j] += 1
end
end
img = imadjustintensity(img)
imshow(img)
```
Given an image similar to the one I generated above, is there a way to estimate the mean and variance of the original isonormal (what I called symmetrical Gaussian) distribution?
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [December 2, 2017, 11:50am UTC](https://discourse.julialang.org/t/multiple-linear-regression-with-weights/7449/7 "2017-12-02T11:50:56Z")
```julia
using Images, Distributions, StatsBase
d = MvNormal([1,2], 2)
img = centered(zeros(13,13))
for rep in 1:1000
i, j = round.(Int, rand(d))
if checkbounds(Bool, img, i, j)
img[i, j] += 1
end
end
# your code up to here, except I am not loading ImageView because I don't
# have it installed, and I don't adjust the image, I don't know what it does
mid_x = indices(img, 1)
mid_y = indices(img, 2)
x = repeat(mid_x; outer = length(mid_y))
y = repeat(mid_y; inner = length(mid_x))
μ1 = mean(x, Weights(vec(img)))
μ2 = mean(y, Weights(vec(img)))
Σ = var(hcat(x, y), Weights(vec(img)), 1)
```
You can probably make the construction of `x` and `y` more sophisticated, but this is the gist.
---
### Author:  [@yakir12](https://discourse.julialang.org/u/yakir12)
#### Post date: [December 2, 2017, 11:55am UTC](https://discourse.julialang.org/t/multiple-linear-regression-with-weights/7449/8 "2017-12-02T11:55:46Z")
Now I get. Right awesome! And that’s a really fast solution too! Thaaaaank you @Tamas_Papp!
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [December 2, 2017, 12:02pm UTC](https://discourse.julialang.org/t/multiple-linear-regression-with-weights/7449/9 "2017-12-02T12:02:50Z")
You are welcome. You can make it more elegant by replacing the last 5 lines with
```julia
XY = hcat(vec(((x,y)->x).(mid_x, mid_y')),
vec(((x,y)->y).(mid_x, mid_y')))
μ, Σ = mean_and_var(XY, Weights(vec(img)), 1)
```
Probably someone will come along and show an even better method for constructing that matrix.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [December 2, 2017, 2:59pm UTC](https://discourse.julialang.org/t/multiple-linear-regression-with-weights/7449/10 "2017-12-02T14:59:20Z")
Sorry for the noise, could not resist:
```julia
XY = hcat(vec(Base.vect.(mid_x, mid_y'))...)'
```
is the most compact syntax I could find. Promise I will stop now 😄
Is there any alternative to `Base.vect` for broadcasting? This is what `[...]` is lowered to, but I could not make that broadcast.
---
### Author:  [@yakir12](https://discourse.julialang.org/u/yakir12)
#### Post date: [December 2, 2017, 3:01pm UTC](https://discourse.julialang.org/t/multiple-linear-regression-with-weights/7449/11 "2017-12-02T15:01:48Z")
Awesome!
For me, the iterative portion of all that code is the
```julia
μ, Σ = mean_and_var(XY, Weights(vec(img)), 1)
```
I can luckily create the `XY` matrix once, and then just reuse it.
Thanks again for all the great help!