Decomposition of MSE, an example of Standard Deviation
Decomposition of MSE, an example of Standard Deviation
Embedded Files
Given such a basic setting (a setting that appeared in my Research Module of Econometrics and Statistics):
Here sigma_hat is a biased estimator, and sigma_tild is an unbiased estimator of standard deviation, but actually sigma_tild as standard deviation is also biased, according to https://en.wikipedia.org/wiki/Variance, thus here I also show the Bias^2 of sigma_tild^2, which is the variance.
Here as shown above, the blue lines are the Bias^2 of sigma_hat and sigma_tild, and we can see that just like normally believed, the "unbiased" one (which is called Bias^2 (n - 1) ) really has a lower bias square than the "biased" one. However, as the sample size "n_seq" is small, they are all biased. Nevertheless, the green triangles, which stand for the unbiased variance estimator, keep being really unbiased, which equals zero.
Meanwhile, the MSE of the two estimators are similar, and we can observe that MSE = Var + Bias^2.
R Code (click to watch)
R Code (click to watch)
muf <- function(X,n){
output <- sum(X)/n
return(output)
}
sigmafn <- function(X,n){
dev <- X-mean(X)
dev <- 1/n*sum(dev^2)
output <- sqrt(dev)
return(output)
}
sigmafn1 <- function(X,n){
dev <- X-mean(X)
dev <- 1/(n-1)*sum(dev^2)
output <- sqrt(dev)
return(output)
}
sigma2fn1 <- function(X,n){
dev <- X-mean(X)
dev <- 1/(n-1)*sum(dev^2)
return(dev)
}
n_seq <- seq(2,30)
mu=1
sigma <- 1
Times <- 10000
resultsmu <- matrix(NA,Times,length(n_seq))
resultsigma <- matrix(NA,Times,length(n_seq))
resultsigman <- matrix(NA,Times,length(n_seq))
resultsigman2 <- matrix(NA,Times,length(n_seq))
set.seed(527)
for (i in 1:Times){
for (j in 1:length(n_seq)) {
X <- rnorm(n_seq[j],1,1)
resultsmu[i,j] <- muf(X,n_seq[j])
resultsigma[i,j] <- sigmafn(X,n_seq[j])
resultsigman[i,j] <- sigmafn1(X,n_seq[j])
resultsigman2[i,j] <- sigma2fn1(X,n_seq[j])
}
}
MSEf <- function(si,Ti){
dev2 <- sum((si-1)^2)
output <- dev2/Ti
return(output)
}
Varf <- function(si,Ti){
mu <- mean(si)
dev2 <- sum((si-mu)^2)
out <- dev2/Ti
return(out)
}
Bias2f <- function(si,Ti){
bia <- (mean(si)-1)^2
return(bia)
}
MSE <- c()
Var <- c()
Bias2 <- c()
MSEn <- c()
Varn <- c()
Bias2n <- c()
Bias22 <- c()
for (i in 1:length(n_seq)) {
MSE[i] <- MSEf(resultsigma[,i],Times)
MSEn[i] <- MSEf(resultsigman[,i],Times)
Var[i] <- Varf(resultsigma[,i],Times)
Varn[i] <- Varf(resultsigman[,i],Times)
Bias2[i] <- Bias2f(resultsigma[,i],Times)
Bias2n[i] <- Bias2f(resultsigman[,i],Times)
Bias22[i] <- Bias2f(resultsigman2[,i],Times)
}
plot(n_seq,Bias22)
lines(n_seq,Bias2n)
plot(n_seq,MSE,ylim = c(0,0.5),type="l",
main = "Decomposition of MSE",lwd=2)
lines(n_seq,Var,col="red",lwd=2)
lines(n_seq,Bias2,col="blue",lwd=2)
points(n_seq,MSEn,pch=2,"o")
points(n_seq,Varn,col="red",pch=3,"o")
points(n_seq,Bias2n,col="blue",pch=2,"o")
points(n_seq,Bias22,col="green",pch=6,"o")
abline(0,0,lty=2)
legend("topright",c("MSE","Var","Bias^2",
"MSE (n-1)","Var (n-1)","Bias^2 (n-1)","Bias^2 of Variance"),
col=c("black","red","blue","black","red","blue","green"),
pch=c(19,19,19,2,3,2,6),lwd=c(2,2,2,1,1,1,1))
Page updated
Google Sites
Report abuse