-
Notifications
You must be signed in to change notification settings - Fork 18
/
luR.r
34 lines (34 loc) · 1.18 KB
/
luR.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
luR <-
function(n=5,thetaL=0,thetaU=1,Nsim=1000000)
{
dump("luR","c:\\StatBook\\luR.r")
print(date())
X=matrix(runif(Nsim*n,min=thetaL,max=thetaU),nrow=Nsim)
xb=apply(X,1,mean)
sde=apply(X,1,sd)
mm.L=xb-sqrt(3)*sde;mm.U=xb+sqrt(3)*sde
MSE1=MSE2=MSE3=matrix(ncol=2,nrow=2)
MSE1[1,1]=mean((mm.L-thetaL)^2)
MSE1[1,2]=MSE1[2,1]=mean((mm.L-thetaL)*(mm.U-thetaU))
MSE1[2,2]=mean((mm.U-thetaU)^2)
print("MSE1:");print(MSE1);cat("Total MSE =",sum(diag(MSE1)),"\n")
xmin=apply(X,1,min)
xmax=apply(X,1,max)
MSE2[1,1]=mean((xmin-thetaL)^2)
MSE2[1,2]=MSE2[2,1]=mean((xmin-thetaL)*(xmax-thetaU))
MSE2[2,2]=mean((xmax-thetaU)^2)
print("MSE2:");print(MSE2);cat("Total MSE =",sum(diag(MSE2)),"\n")
xmin.L=(n+1)/(n^2+1)*(n*xmin-xmax)
xmax.U=(n+1)/(n^2+1)*(n*xmax+xmin)
MSE3[1,1]=mean((xmin.L-thetaL)^2)
MSE3[1,2]=MSE3[2,1]=mean((xmin.L-thetaL)*(xmax.U-thetaU))
MSE3[2,2]=mean((xmax.U-thetaU)^2)
print("MSE3:");print(MSE3);cat("Total MSE =",sum(diag(MSE3)),"\n")
print("Eigenvalues of the difference MSE1-MSE3:")
print(eigen(MSE1-MSE3)$values)
print("Eigenvalues of the difference MSE2-MSE3:")
print(eigen(MSE2-MSE3)$values)
d=as.data.frame(cbind(MSE1,MSE2,MSE3))
print(d)
return(date())
}