2 Heteroskedasticity

Suppose the noise variance is itself variable.

Figure 2: Scatter-plot of n = 150 data points from the above model. (Here X isGaussian with mean 0 and variance 9.) Grey: True regression line. Dashed: ordinaryleast squares regression line.


 n=150
      x=rnorm(n,0,3)
      y = 3-2*x + rnorm(n,0,sapply(x,function(x){1+0.5*x^2}))
      
      # Plot the data
      plot(x,y)
      # Plot the true regression line
      abline(a=3,b=-2,col="grey")
      # Fit by ordinary least squares
      fit.ols = lm(y~x)
      # Plot that line
      abline(fit.ols,lty="dashed")

      





par(mfrow=c(1,2))
plot(x,residuals(fit.ols))
plot(x,(residuals(fit.ols))^2)

par(mfrow=c(1,1))




Figure 3: Residuals (left) and squared residuals (right) of the ordinary least squares
regression as a function of x. Note the much greater range of the residuals at large
absolute values of x than towards the center; this changing dispersion is a sign of

heteroskedasticity.



# Generate more random samples from the same model and the same x values,
      # but different y values
      # Inputs: number of samples to generate
      # Presumes: x exists and is defined outside this function
      # Outputs: errors in linear regression estimates
      ols.heterosked.example = function(n) {
        y = 3-2*x + rnorm(n,0,sapply(x,function(x){1+0.5*x^2}))
        fit.ols = lm(y~x)
        # Return the errors
        return(fit.ols$coefficients - c(3,-2))
      }
      # Calculate average-case errors in linear regression estimates (SD of
      # slope and intercept)
      # Inputs: number of samples per replication, number of replications (defaults
      # to 10,000)
      # Calls: ols.heterosked.example
      # Outputs: standard deviation of intercept and slope
      ols.heterosked.error.stats = function(n,m=10000) {
        ols.errors.raw = t(replicate(m,ols.heterosked.example(n)))
        # transpose gives us a matrix with named columns
        intercept.sd = sd(ols.errors.raw[,"(Intercept)"])
        slope.sd = sd(ols.errors.raw[,"x"])
        return(list(intercept.sd=intercept.sd,slope.sd=slope.sd))
      }
      
      

      $intercept.sd
[1] 0.6180515

$intercept.sd

[1] 0.6180515



2.1 Weighted Least Squares as a Solution to Heteroskedasticity


   # Plot the data
      plot(x,y)
      # Plot the true regression line
      abline(a=3,b=-2,col="grey")
      # Fit by ordinary least squares
      fit.ols = lm(y~x)
      # Plot that line
      abline(fit.ols,lty="dashed")
      fit.wls = lm(y~x, weights=1/(1+0.5*x^2))
      abline(fit.wls,lty="dotted")
      

      


Figure 6: Figure 2, plus the weighted least squares regression line (dotted)





### As previous two functions, but with weighted regression
# Generate random sample from model (with fixed x), fit by weighted least
# squares
# Inputs: number of samples
# Presumes: x fixed outside function
# Outputs: errors in parameter estimates
wls.heterosked.example = function(n) {
y = 3-2*x + rnorm(n,0,sapply(x,function(x){1+0.5*x^2}))
fit.wls = lm(y~x,weights=1/(1+0.5*x^2))
# Return the errors
return(fit.wls$coefficients - c(3,-2))
}
# Calculate standard errors in parameter estiamtes over many replications
# Inputs: number of samples per replication, number of replications (defaults
# to 10,000)
# Calls: wls.heterosked.example
# Outputs: standard deviation of estimated intercept and slope
wls.heterosked.error.stats = function(n,m=10000) {
wls.errors.raw = t(replicate(m,wls.heterosked.example(n)))
# transpose gives us a matrix with named columns
intercept.sd = sd(wls.errors.raw[,"(Intercept)"])
slope.sd = sd(wls.errors.raw[,"x"])
return(list(intercept.sd=intercept.sd,slope.sd=slope.sd))

}



2.2 Some Explanations for Weighted Least Squares





  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值