Adv_hw_2

1

n=25;theta=2*qt(df=3,p=0.75)/1.34;set.seed(10);
alpha=0.05;z=-qnorm(alpha/2);B=1000;
count.1=count.2=count.3=length.1=length.2=length.3=rep(0,1000);

for( k in 1: 1000){
X=rt(n=n,df=3);

#Jacknife
T=T.tilde=NULL;
for( i in 1:n){ 
T[i]=IQR(X[-i])/1.34 
T.tilde[i] = n*IQR(X)/1.34 - (n-1)*T[i]
}
U=mean(T.tilde)+z*sqrt(var(T.tilde)/n)
L=mean(T.tilde)-z*sqrt(var(T.tilde)/n)
if((U-theta)*(L-theta)<0){count.1[k] = 1}
length.1[k] = 2*z*sqrt(var(T.tilde)/n);

#Normal interval
nonpara=NULL;
for( i in 1:B){
nonpara[i] = IQR(sample(X, size=n, replace=TRUE))/1.34
}
U=mean(nonpara)+z*sqrt(var(nonpara));
L=mean(nonpara)-z*sqrt(var(nonpara));
if ((U-theta)*(L-theta)<0){count.2[k] = 1}
length.2[k] = 2*z*sqrt(var(nonpara));

#Percentile interval
U=quantile(x=nonpara, probs=c(alpha/2,(1-alpha/2)))[2];
L=quantile(x=nonpara, probs=c(alpha/2,(1-alpha/2)))[1];
if ((U-theta)*(L-theta)<0){count.3[k] = 1}
length.3[k] = U-L;

}
mean(count.1);mean(count.2);mean(count.3);

## [1] 0.377

## [1] 0.968

## [1] 0.978

mean(length.1);mean(length.2);mean(length.3);

## [1] 1.249

## [1] 1.351

## [1] 1.314

1번 결과입니다. theta가 IQR에 비례하는데 이 경우는 Jacknife가 안좋다.

2-(a)

theta=1;n=50;B=1000;brk=20;set.seed(10);
X=runif(n=n,min=0,max=theta);theta.hat=max(X)
para=nonpara=NULL;
for( i in 1:B){
nonpara[i] = max(sample(X, size=n, replace=TRUE))
para[i] = max(runif(n=n,min=0,max=theta.hat))
}
t=seq(1:1000)/1000;pdf=(B/brk)*n*(t^(n-1))
par(mfrow=c(1,2))
hist(nonpara,xlim=c(0,1),breaks=brk)
lines(t, pdf,lwd=1.5)
hist(para,xlim=c(0,1),breaks=brk)
lines(t, pdf,lwd=1.5)

실제 pdf를 실선으로 비교했다. nonparametric bootstrap 방법으로는 전혀 비슷해보이지 않으나 parametric은 매우 비슷해보인다.

6

setwd("C:\\Users\\snu\\Desktop\\YC\\SNU\\2014-2\\Advanced Statistics\\")
data.raw=read.table(file="data.txt",header=TRUE,sep=",")
attach(data.raw);

## The following objects are masked _by_ .GlobalEnv:
## 
##     AGE, SAL

AGE=as.vector(AGE);SAL=as.vector(SAL);AGE=AGE[-30];SAL=as.numeric(SAL[-30])

source("ch6.r")
require(MASS)

n=length(AGE);
hist.h = cv.hist.fun(SAL)$mbest;
sigma.hat = min(sd(SAL),IQR(SAL)/1.34)
h.normal = 1.06*sigma.hat/n^(0.2)
h.cv = ucv(SAL)

f1=density(SAL,width=h.normal)
f2=density(SAL,width=h.cv)

par(mfrow=c(1,3))
truehist(SAL,nbins=hist.h)
plot(f1)
plot(f2)

Confidence band는 수업시간에 다루지 않았으므로 실행을 못했다. 그 외는 그래프에서 볼 수 있다.