# Improve Model Fit with Weights

This example shows how to fit a polynomial model to data using both the linear least-squares method and the weighted least-squares method for comparison.

Generate sample data from different normal distributions by using the `randn`

function.

rng("default") % For reproducibility rnorm = []; idx = []; for k=1:20 r = k*randn([20,1]) + (1/20)*(k^3); rnorm = [rnorm;r]; idx = [idx;ones(20,1).*k]; end

The dependent variable `rnorm`

contains sample data from 20 normal distributions. The independent variable `idx`

contains integers indicating whether two elements in `rnorm`

are sampled from the same normal distribution.

Fit a third-degree polynomial model to `idx`

and `rnorm`

. Return information about the coefficient estimates and the algorithm used to fit the model.

`[p3fit,~,fitinfo] = fit(idx,rnorm,"poly3")`

p3fit = Linear model Poly3: p3fit(x) = p1*x^3 + p2*x^2 + p3*x + p4 Coefficients (with 95% confidence bounds): p1 = 0.05156 (0.0438, 0.05932) p2 = -0.03993 (-0.2875, 0.2076) p3 = 0.1418 (-2.124, 2.408) p4 = 0.0462 (-5.585, 5.678)

`fitinfo = `*struct with fields:*
numobs: 400
numparam: 4
residuals: [400x1 double]
Jacobian: [400x4 double]
exitflag: 1
algorithm: 'QR factorization and solve'
iterations: 1

`p3fit`

contains the estimates for the model coefficients with 95% confidence intervals. The default fitting method for fitting a polynomial model is linear least squares. `fitinfo`

contains information about the fitting algorithm used to fit the coefficients to the data. The error in the data can be estimated by the residuals stored in `fitinfo`

.

Plot the residuals using a stem plot.

stem(idx,fitinfo.residuals) xlabel("idx") ylabel("residuals")

The plot of the residuals shows that their variance increases as the values in `idx`

increase. The nonconstant variances across the different values of `idx`

indicate that the weighted least-squares fitting method is more appropriate for calculating the model coefficients.

Create a vector of zeros for later storage of the weights.

W = zeros(length(rnorm),1);

The weights you supply transform the residual variances so that they are constant for different values of `idx`

. Define the weight for each element in `rnorm`

as the reciprocal of the residual variance for the corresponding value in `idx`

. Then fit the model with the weights.

for k=1:20 rnorm_idx = rnorm(idx==k); recipvar = 1/var(rnorm_idx); w = (idx==k).*recipvar; W = W+w; end [wp3fit,~,wfitinfo] = fit(idx,rnorm,"poly3","Weights",W)

wp3fit = Linear model Poly3: wp3fit(x) = p1*x^3 + p2*x^2 + p3*x + p4 Coefficients (with 95% confidence bounds): p1 = 0.04894 (0.04419, 0.0537) p2 = 0.03601 (-0.08744, 0.1595) p3 = -0.4262 (-1.253, 0.4009) p4 = 0.9836 (-0.1959, 2.163)

`wfitinfo = `*struct with fields:*
numobs: 400
numparam: 4
residuals: [400x1 double]
Jacobian: [400x4 double]
exitflag: 1
algorithm: 'QR factorization and solve'
iterations: 1

`wp3fit`

and `wfitinfo`

contain the results of the weighted least-squares fitting.

Display `p3fit`

, `wp3fit`

, and `rnorm`

in the same plot.

plot(p3fit,idx,rnorm) hold on plot(wp3fit,"g") xlabel("idx") ylabel("rnorm") legend(["rnorm","linear least-squares fit","weighted least-squares fit"]) hold off

The plot shows `wp3fit`

closely tracking `p3fit`

.

You can determine whether `wp3fit`

is a better fit than `p3fit`

by plotting the residuals.

stem(idx,wfitinfo.residuals) xlabel("idx") ylabel("residuals")

The output shows that the `wp3fit`

residuals are smaller than the `p3fit`

residuals. The variances of the `wp3fit`

residuals are also more similar for different values of `idx`

than the variances of the `p3fit`

residuals.