David L Lorenz wrote:
> We have a situation where we know the variance of observations for a
> linear regression problem and need confidence limits on the parameter
> estimates. The documentation in the NOTES section for lm states "In
> addition, S-PLUS does not currently support weighted regression when
> the absolute precision of the observations is known. This situation
> arises often in physics and engineering, when the uncertainty
> associated with a particular measurement is known in advance due to
> properties of the measuring procedure or device. If you know the
> absolute precision of your observations, it is possible to supply
> them to the weights argument. This computes the correct coefficients
> for your model, but the standard errors and other inference tools
> will be incorrect."
>
> I have searched the web, done bibliographic searches, and even
> searched the documentation for a competing product and I found
> nothing regarding the solution of this problem. I would be surprised
> if this problem has not been solved because it should be a high
> priority in physics. At least I remember it being an issue in physics
> classes, if not in the engineering classes I took.
>
> Is anyone aware of an approach to estimating the confidence limits of
> parameter estimates when the variance of the observations is known? I
> think I can see a solution for the case of equal known variance, but I
> don't know that I have the ability to apply that to unequal
> variances.Thanks.
The solution to the problem is trivial ***in theory***. Numerical
considerations might bite you; I dunno.
Let your vector of observations be Y and assume
Y ~ N(X beta, Delta)
where Delta = diag{sigma_1^2, ... sigma_n^2) where the sigma_i`
in the scenario you envisiage are ***known*** values.
Let W = Delta^{-1/2} Y, then W ~ N(Delta^{-1/2} X beta, I), where
I is the n x n identity. Let M = Delta^{-1/2} X for short, so that
the model becomes W ~ N(M beta, I).
The estimate of beta is
beta-hat = (M'M)^{-1}M'W
(where M' means the transpose of M) and the covariance matrix of
beta-hat, on which you can base any inferential procedure you wish,
is (M'M)^{-1}.
Modulo numerical gotchas --- which are always lurking about --- the
calculations can be done in raw S/R in a few lines. If your design
matrix X is reasonably well behaved and the number of columns
(= length(beta)) is not too large, you should be able to roll.
Elsewise something more numerically sophisticated (Cholesky or Q-R
decompositions?) will have to brought into play.
cheers,
Rolf Turner
rolf@math.unb.ca
|