. Surajit Ray Ray SAMSI, June 3 2005 - slide #1 Outline Outline - - PowerPoint PPT Presentation

Ray SAMSI, June 3 2005 - slide #1

Alternative Statistical Models

.

Surajit Ray

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #2

Outline

■ Recap of (ordinary) least-squares ■ Violations of statistical assumptions ■ Weighted least-squares ■ Generalized least-squares ■ Concluding comments

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #3

Recap of (ordinary) least-squares

Model for data (t j,y j), j = 1,...,n:

yj = y(t j;q)+ε j

■ y(t j;q) is deterministic model, with parameters q. ■ εj are random errors.

Goal of the inverse problem: estimate q. Standard statistical assumptions for the model:

• 1. y(t j;q) is correct model ⇒ mean of ε j is 0 for all j.
• 2. Variance of ε j is contant for all j, equal to σ2.
• 3. Error terms ε j, εk are independent for j = k.

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #4

OLS estimation

■ Minimize

J(q) =

n

∑

i=1

|yi −y(ti;q)|2 (1)

in q, to give

qols.

■ Estimate σ2 by

• σ2
• ls =

1 n− pJ( qols)

where p = dim(q).

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #5

OLS Estimation

■ converges to q as n increases ■ makes efficient use of the data, i.e., has small

standard error

■ approximate s.e.(

qols,k) = square root of (k,k)

element in Cov(

q) = σ2

• ls
• XTX

−1

where Xr,c =

∂y(tr;q)

∂qc

• evaluated at

qols.

For example, if q = (C,K)T , then

X =      

∂y(t1;q) ∂C ∂y(t1;q) ∂K ∂y(t2;q) ∂C ∂y(t2;q) ∂K

. . . . . .

∂y(tn;q) ∂C ∂y(tn;q) ∂K

      .

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #6

Violations of statistical assumptions

Compute residuals, r j = y j −y(t j;

qols), plot against t j:

0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 Residual vs. Time (Frequency 57 Hz − Damped) Time, tj Residual, rj

Figure 1:

Residual plot for the (damped) spring-mass-dashpot model, fitted using OLS.

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #7

Analysis of Residual plot

• 1. Do we have the correct deterministic model?
• 2. Is variance of ε j constant across time range? No!
• 3. Are errors independent? No!

Implications:

qols is no longer a good estimator for q.

Assuming that answer to #1 is “Yes”, how can we change our statistical model assumptions to better model reality?

■ Transform the data (e.g., log transform)? ■ Explicitly model nonconstant variance, and

correlations between measurements.

■ Incorporate this into the estimation method.

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #8

Weighted least-squares

(a) Deal with nonconstant variance: Assume Var(ε j) = σ2

w j , j = 1,...,n

for known w j.

w j large ⇔ observation (t j,y j) is of high quality.

Instead of OLS, minimize

• J(q) =

n

∑

i=1

wi |yi −y(ti;q)|2

in q.

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #9

Weighted Least Squares

In practice, don’t know w j:

■ Estimate Var(ε j) from repeated measurements at time

t j: w j = σ2

• σ2

j

.

■ If error is larger for larger |y j|, let w−1

j

= y2

j.

■ Alternatively, assume that w−1

j

= y2(t j;q).

■ Assume some other model for w j, e.g.,

w−1

j

= y2θ(t j;q)

where θ is to be estimated from the data.

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #10

Weighted Least Squares

Deal with correlated observations and nonconstant variance: Let ε = (ε1,ε2,...,εn)T and assume Cov(ε) = σ2V, for known matrix V. Let a = (a1,a2,...,an)T , y(q) = (y(t1;q),...,y(tn;q))T ,

W = V −1.

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #11

Weighted Least Squares

The weighted least-squares (WLS) estimator minimizes

Jwls(q) = {a−y(q)}TW{a−y(q)} (2)

in q, to give

qwls.

If above covariance model holds (together with assumption 1) then

qwls has good properties.

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #12

Generalized least-squares (GLS)

Estimate V from the data too! Model Cov(ε) = σ2V in two stages, based on additional parameters θ, α to be estimated from the data. (a) Model Var(ε j). For example, Var(ε j) = σ2y2θ(t j;q), where θ = θ (i.e., scalar). Define diagonal matrices:

G(q,θ) = diag{y2θ(t1;q),y2θ(t2;q),...,y2θ(tn;q)},

and

{G(q,θ)}1/2 = diag{yθ(t1;q),yθ(t2;q),...,yθ(tn;q)}.

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #13

GLS

(b) Model Corr(ε j,εk) for j = k. For example, Corr(ε j,εk) = α|j−k|, where α = α (i.e., scalar) and |α| < 1. Organize into a matrix: Corr(ε) = Γ(α) =

      1 α α2 ··· αn−1 α 1 α ··· αn−2

. . . . . . . . .

···

. . .

αn−1 αn−2 αn−3 ··· 1       .

Put pieces together, write Cov(ε) = σ2{V(q,θ,α)} = σ2{G(q,θ)}1/2 Γ(α){G(q,θ)}1/2.

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #14

GLS

Typical GLS algorithm follows three steps: (i) Estimate q with OLS. Set

qgls = qols.

(ii) Estimate (somehow) (

θ, α). Plug into the model for V

to form estimated weight matrix,

• W = {V(

qgls, θ, α)}−1.

(iii) Minimize (approximate) WLS objective function (2) in

q, namely,

• Jwls(q) = {a−y(q)}T

W{a−y(q)}.

Set the minimizing value equal to

qgls.

Return to step (ii). Iterate until “convergence”.

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #15

GLS

Finally, compute estimate of σ2:

• σ2

gls =

1 n− p{a−y( qgls)}T{V( qgls, θ, α)}−1{a−y( qgls)}.

Methods for estimating (

θ, α) [step (ii)] are beyond the

scope of this talk. Commonly used methods include

■ psuedolikelihood (PL) ■ restricted maximum likelihood (REML).

Outline Recap of (ordinary) least-squares OLS estimation OLS Estimation Violations of statistical assumptions Analysis of Residual plot Weighted least-squares Weighted Least Squares Weighted Least Squares Weighted Least Squares Generalized least-squares (GLS) GLS GLS GLS Concluding comments Ray SAMSI, June 3 2005 - slide #16

Concluding comments

■ Does this approach work better? Are our new

statistical assumptions met?

■ Check! Calculate the GLS weighted residuals,

rgls = {V( qgls, θ, α)}−1/2{a−y( qqls)}

and plot against t j.

■ What about standard errors for

qgls?

■ Approximate s.e.(

qgls,k) = square root of (k,k)

element in

• σ2

gls

• XT{V(

qgls, θ, α)}−1X −1

where Xr,c =

∂y(tr;q)

∂qc

• evaluated at

qgls.