LoginSignup
0
0

More than 3 years have passed since last update.

@kozakai-ryouta(亮太 小酒井)

サンプルサイズの決め方　永田靖先生　オーム社

Last updated at Posted at 2018-10-25

#サンプルサイズの大きさ

#estimate sample size

n=10000

samples=sample(c(0:1),n,replace=T,prob=c(1-1/3,1/3))

data=data.frame(num=1:n,samples=samples) %>% mutate(freq=cumsum(samples)/num)

data=data %>% mutate(cumsum_samples=cumsum(samples))

#l:許容発生誤差回数、N:必要サイズ、k:求める発生回数

error_values=function(l,N,k){

t=(l*sqrt(N))/sqrt(2*k*(N-k))  

value<-(sqrt(N)/sqrt(2*pi*k*(N-k)))*exp(-t^2)

f=function(y){

return(exp(-y^2))

}

return(c(t,value,(2/sqrt(pi))*integrate(f,0,t)$value))  

}

plot_data=data.frame(l=0,k=0,value=0,prob=0)

L=100;K=100;N=1000

for(l in 1:L){
for(k in 1:K){  

val=error_values(l,N,k) 

plot_data_sub=data.frame(l=l,k=k,value=val[2],prob=val[3])   

plot_data=rbind(plot_data,plot_data_sub)

}
}

plot_data=tail(plot_data,nrow(plot_data)-1)

library(rgl)

plot3d(plot_data$l,plot_data$k,plot_data$prob)

plot3d(plot_data$l,plot_data$k,plot_data$value)



#非心t分布での検定(H1:mu!=mu0)

library(dplyr)

n=100000

sample=sample(c(0:1),n,replace=T,prob=c(0.24,0.76))

#データの値から近似値を推定

mu0=0.76

values=function(w,N,samples){

lambda=(mean(samples)-mu0)/sd(samples)  

z<-(w*(1-1/(4*(N-1)))-lambda)/sqrt((1+(w^2)/(2*(N-1))))  

return(c(z,lambda)  )

}

#t値の誤差がp以下である確率(t=(x_bar-mu)/sqrt(V/n))

alpha=0.01;delta=0.005

plot_data=data.frame(z=seq(10,n,10),prob=0,lambda=0,n=0) 

for(j in 1:nrow(plot_data)){

z=plot_data$z[j]  

plot_data$prob[j]=pnorm(values(qt(alpha/2,df=z-1),z,sample[1:z])[1])+1-pnorm(values(qt(1-alpha/2,df=z-1),z,sample[1:z])[1])

plot_data$lambda[j]=values(qt(1-alpha/2,df=z-1),z,sample[1:z])[2]

plot_data$n[j]=((qnorm(alpha/2)-qnorm(plot_data$prob[j]))/delta)^2+(qnorm(alpha/2)^2)/2

}

plot(plot_data$z,plot_data$prob,type="o")

plot(plot_data$z,plot_data$lambda,type="o")



#母分散の収束サイズ

library(dplyr)

n=100000

sample=sample(c(0:1),n,replace=T,prob=c(0.24,0.76))

#データの値から近似値を推定

num=2

sigma0=round(0.76*(0.24)*10^num)/10^num #p*(1-p)

alpha=0.01;delta=0.999

plot_data=data.frame(z=seq(10,n,10),prob=0,delta=0,n=0,var=0) 


for(j in 1:nrow(plot_data)){

z=plot_data$z[j]

plot_data$var[j]=var(sample[1:z])

plot_data$delta[j]=floor((10^num)*plot_data$var[j]/sigma0)/(10^num)

plot_data$prob[j]=pchisq(qchisq(alpha/2,df=z-1)/(plot_data$delta[j]),df=z-1)+1-pchisq(qchisq(1-alpha/2,df=z-1)/(plot_data$delta[j]),df=z-1)

plot_data$n[j]=(1/2)*((qnorm(alpha/2)-delta*qnorm(plot_data$prob[j]))/(delta-1))^2+3/2

}

plot(plot_data$z,plot_data$prob)

#pchisq(qchisq(alpha/2,df=z-1)/(var(sample[1:z])/sigma0),df=z-1)+1-pchisq(qchisq(1-alpha/2,df=z-1)/(var(sample[1:z])/sigma0),df=z-1)



#p.44

library(dplyr)

n=30

X=rnorm(n)

S=sum((X-mean(X))^2)

sigma0=5

pthi=n-1

kai0=S/sigma0

alpha=0.05

upper=qchisq(1-alpha/2,pthi)

under=qchisq(alpha/2,pthi)

#power1:両側、power2:sigma>sigma0、power3:sigma<sigma0

data=data.frame(delta=seq(0,5,0.01),power1=0,power2=0,power3=0)

for(j in 1:nrow(data)){

data$power1[j]=(pchisq(qchisq(alpha/2,pthi)/(data$delta[j]^2),pthi))+1-pchisq(qchisq(1-alpha/2,pthi)/(data$delta[j]^2),pthi)  

data$power2[j]=1-pchisq(qchisq(1-alpha,pthi)/(data$delta[j]^2),pthi)  

data$power3[j]=pchisq(qchisq(alpha,pthi)/(data$delta[j]^2),pthi)

}

#推定されたdelta:alpha,beta,pthiは既知

alpha=0.01;beta=0.1

delta2=qchisq(1-alpha,pthi)/qchisq(beta,pthi)

#必要サンプルサイズ

#両側仮説の場合

#delta>1

delta=2;alpha=0.05;beta=0.01

n=((qnorm(1-alpha/2)-delta*qnorm(beta))/(delta-1))^2/2+3/2

#delta<1

delta=0.1;alpha=0.05;beta=0.01

n=((qnorm(1-alpha/2)-delta*qnorm(beta))/(delta-1))^2/2+3/2

#片側仮説:sigma>sigma0

delta=0.1;alpha=0.05;beta=0.01

n=((qnorm(1-alpha/2)-delta*qnorm(beta))/(delta-1))^2/2+3/2

#片側仮説:sigma<sigma0

delta=0.1;alpha=0.05;beta=0.01

n=((qnorm(1-alpha/2)-delta*qnorm(beta))/(delta-1))^2/2+3/2

#p.59 正規近似　

#(A) fisher

hist(sqrt(2*rchisq(10000,pthi))-sqrt(2*pthi-1))

#(B) wilson-hilferty

hist(sqrt((9*pthi)/2)*((rchisq(10000,pthi)/pthi)^(1/3)-(1-2/(9*pthi))))






#非心t分布での検定(H1:mu!=mu0) p.67~69

library(dplyr)

n=100000

sample=sample(c(0:1),n,replace=T,prob=c(0.24,0.76))

#データの値から近似値を推定

mu0=0.76

values=function(w,N,samples){

lambda=(mean(samples)-mu0)/sd(samples)  

z<-(w*(1-1/(4*(N-1)))-lambda)/sqrt((1+(w^2)/(2*(N-1))))  

return(c(z,lambda)  )

}

#t値の誤差がp以下である確率(t=(x_bar-mu)/sqrt(V/n))

alpha=0.01;delta=0.005

plot_data=data.frame(z=seq(10,n,10),prob=0,lambda=0,n=0) 

for(j in 1:nrow(plot_data)){

z=plot_data$z[j]  

plot_data$prob[j]=pnorm(values(qt(alpha/2,df=z-1),z,sample[1:z])[1])+1-pnorm(values(qt(1-alpha/2,df=z-1),z,sample[1:z])[1])

plot_data$lambda[j]=values(qt(1-alpha/2,df=z-1),z,sample[1:z])[2]

plot_data$n[j]=((qnorm(alpha/2)-qnorm(plot_data$prob[j]))/delta)^2+(qnorm(alpha/2)^2)/2

}

plot(plot_data$z,plot_data$prob,type="o")

plot(plot_data$z,plot_data$lambda,type="o")



library(dplyr)

n=100

X=rnorm(n,mean=1,sd=7)

mu0=3;pthi=n-1

V=sum((X-mean(X))^2)/(n-1)

t=(mean(X)-mu0)/sqrt(V/n)

#power1:両側検定、power2:片側(mu>mu0)、power3:片側(mu<mu0)

data=data.frame(delta=seq(-1,1,0.01),power1=0,power2=0,power3=0)

for(j in 1:nrow(data)){


val1=qt(alpha/2,pthi);val2=qt(1-alpha/2,pthi)

val1=(val1*(1-1/(4*pthi))-sqrt(n)*data$delta[j])/sqrt(1+val1^2/(2*pthi))

val2=(val2*(1-1/(4*pthi))-sqrt(n)*data$delta[j])/sqrt(1+val2^2/(2*pthi))

data$power1[j]=pnorm(val1)+1-pnorm(val2)

val1=qt(alpha,pthi);val2=qt(1-alpha,pthi)

val1=(val1*(1-1/(4*pthi))-sqrt(n)*data$delta[j])/sqrt(1+val1^2/(2*pthi))

val2=(val2*(1-1/(4*pthi))-sqrt(n)*data$delta[j])/sqrt(1+val2^2/(2*pthi))

data$power3[j]=pnorm(val2);data$power2[j]=1-pnorm(val1)

}

#必要サンプルサイズ　p.75

#両側検定の場合

alpha=0.05;beta=0.01;delta=0.1

n=((qnorm(1-alpha/2)-qnorm(beta))/(delta))+qnorm(1-alpha/2)^2/2

#片側検定(mu>mu0)

alpha=0.05;beta=0.01;delta=0.1

n=((qnorm(1-alpha/2)-qnorm(beta))/(delta))+qnorm(1-alpha/2)^2/2

#片側検定(mu<mu0)

alpha=0.05;beta=0.01;delta=0.1

n=((qnorm(1-alpha/2)-qnorm(beta))/(delta))+qnorm(1-alpha/2)^2/2


#p.79 非心t分布の定義とその性質

n=1000;delta=0.1

lambda=sqrt(n)*delta

pthi=n-1

y=rnorm(n,lambda,1)

t_dash=y/sqrt(rchisq(n,pthi)/pthi)

#非心ｔ分布プロット


hist(t_dash)

plot(seq(0,6,0.1),dt(seq(0,6,0.1),df=pthi,ncp=lambda))


c_star=1-1/(4*pthi);c_star2=1/(2*pthi)


P=0.05

1-pnorm((w*c_star-lambda)/sqrt(1+c_star2*w^2))

data=data.frame(lambda=seq(0,10,0.01),power1=0,power2=0,power3=0)

P=0.05

for(j in 1:nrow(data)){


w1=qt(P/2,df=pthi,ncp=lambda);w2=qt(1-P/2,df=pthi,ncp=lambda)

data$power1[j]=pnorm((w1*c_star-data$lambda[j])/sqrt(1+c_star2*w1^2))+1-pnorm((w2*c_star-data$lambda[j])/sqrt(1+c_star2*w2^2))

w=qt(P,df=pthi,ncp=lambda)

data$power2[j]=1-pnorm((w*c_star-data$lambda[j])/sqrt(1+c_star2*w^2))

w=qt(P,df=pthi,ncp=lambda)

data$power3[j]=pnorm((w*c_star-data$lambda[j])/sqrt(1+c_star2*w^2))

}

plot(data$lambda,data$power1)

plot(data$lambda,data$power2)

plot(data$lambda,data$power3)

#p.87

n1=100;n2=150

X1=rnorm(n1,mean=5,sd=5);X2=rnorm(n2,mean=8,sd=3)

V1=sum((X1-mean(X1))^2)/(length(X1)-1)

V2=sum((X2-mean(X2))^2)/(length(X2)-1)

pthi1=n1-1;pthi2=n2-1

f0=V1/V2

alpha=0.05

qf(1-alpha,pthi1,pthi2)

#6.2 検出力の計算方法

data=data.frame(delta=seq(0,5,0.01),power1=0,power2=0,power3=0)

alpha=0.05

for(j in 1:nrow(data)){

data$power1[j]=pf(1/(qf(1-alpha/2,pthi2,pthi1)*data$delta[j]^2),pthi1,pthi2)+1-pf(qf(1-alpha/2,pthi1,pthi2)/(data$delta[j]^2),pthi1,pthi2)  

data$power2[j]=1-pf(qf(1-alpha,pthi1,pthi2)/(data$delta[j]^2),pthi1,pthi2)

data$power3[j]=pf(1/(qf(1-alpha,pthi2,pthi1)*data$delta[j]^2),pthi1,pthi2)

}

plot(data$delta,data$power1)

plot(data$delta,data$power2)

plot(data$delta,data$power3)

#p.96

#H0:sigma1=sigma2,sigma1>sigma2,sigma1<sigma2

alpha=0.05;beta=0.01;delta=2

n=1+((qnorm(1-alpha/2)-qnorm(beta))/(log(delta)))^2

plot(seq(-10,10,0.1),df(seq(-10,10,0.1),pthi1,pthi2))

E_F=pthi2/(pthi2-2);V_F=(2*(pthi1+pthi2-2)*pthi2^2)/(pthi1*(pthi2-4)*(pthi2-2)^2)

#p.105

#H0:mu1=mu2

V=(sum((X1-mean(X1))^2)+sum((X2-mean(X2))^2))/(pthi1+pthi2)

t=(mean(X1)-mean(X2))/sqrt(V*(1/n1+1/n2))

alpha=0.05

qt(1-alpha,pthi1+pthi2)

#H1:mu1!=mu2

delta=0.1;lambda=sqrt((n1*n2)/(n1+n2))*delta

qt(1-alpha,pthi1+pthi2,lambda)


#welch test

V1=sum((X1-mean(X1))^2)/(n1-1);V2=sum((X2-mean(X2))^2)/(n2-1)

t0=(mean(X1)-mean(X2))/sqrt(V1/n1+V2/n2)

pthi_star=(V1/n1+V2/n2)^2/((V1/n1)^2/pthi1+(V2/n2)^2/pthi2)

qt(1-0.05,pthi_star)


alpha=0.05

data=data.frame(delta=seq(-1,1,0.01),power1=0,power2=0,power3=0)

pthi=pthi1+pthi2

for(j in 1:nrow(data)){

lambda=sqrt((n1*n2)/(n1+n2))*data$delta[j]  

val1=qt(alpha,pthi1+pthi2);val2=qt(1-alpha,pthi1+pthi2)  

val1=(val1*(1-1/(4*pthi))-lambda)/sqrt(1+val1^2/(2*pthi))

val2=(val2*(1-1/(4*pthi))-lambda)/sqrt(1+val2^2/(2*pthi))

data$power1[j]=pnorm(val1)+1-pnorm(val2) 

val1=qt(2*alpha,pthi1+pthi2);val2=qt(1-2*alpha,pthi1+pthi2)  

val1=(val1*(1-1/(4*pthi))-lambda)/sqrt(1+val1^2/(2*pthi))

val2=(val2*(1-1/(4*pthi))-lambda)/sqrt(1+val2^2/(2*pthi))

data$power2[j]=pnorm(val1)

data$power3[j]=1-pnorm(val2) 

}

plot(data$delta,data$power1)

plot(data$delta,data$power2)

plot(data$delta,data$power3)

#サンプルサイズの推定　H0:mu1=mu2,mu1>mu2,mu1<mu2

alpha=0.05;beta=0.01;delta=1

n=2*((qnorm(1-alpha/2)-qnorm(beta))/(delta))^2+qnorm(alpha/2)^2/4

#p.122

d=X1-X2

ave_d=mean(d)

V_d=sum((d-ave_d)^2)/(n1+n2-1)

pthi=(n1+n2-1)

t0=ave_d/sqrt(V_d/(n1+n2))

data=data.frame(delta=seq(-1,1,0.01),power=0)

pthi=n1+n2-1

for(j in 1:nrow(data)){

lambda=sqrt(n)*data$delta[j]  

val1=qt(alpha,pthi1+pthi2);val2=qt(1-alpha,pthi)  

val1=(val1*(1-1/(4*pthi))-lambda)/sqrt(1+val1^2/(2*pthi))

val2=(val2*(1-1/(4*pthi))-lambda)/sqrt(1+val2^2/(2*pthi))

data$power[j]=pnorm(val1)+1-pnorm(val2) 


}

plot(data$delta,data$power)

#H0:mu1=mu2

alpha=0.05;beta=0.01;delta=0.1

n=((qnorm(1-alpha/2)-qnorm(beta))/delta)^2+qnorm(1-alpha/2)^2/2

#p.130 一元配置分散分析 分散が既知の場合:sigma0

sigma0=10;n=20

X=t(matrix(c(rnorm(n,2,3),rnorm(n,3,4),rnorm(n,4,5),rnorm(n,8,9)),nrow=n))

S_A=sum((apply(X,1,mean)-mean(X))^2)*n

pthi=nrow(X)-1

a=apply(X,1,mean)-mean(X);mu=mean(X)

kai0=S_A/sigma0

qchisq(1-0.05,pthi)

#検出力

alpha=0.05

data=data.frame(delta=seq(0.01,3,0.01),power=0)

for(j in 1:nrow(data)){

lambda=n*data$delta[j]  

C=(pthi+2*lambda)/(pthi+lambda)

pthi_star=(pthi+lambda)^2/(pthi+2*lambda)

W=qchisq(1-alpha,df=pthi-1)

data$power[j]=1-pnorm(sqrt((2*W)/C)-sqrt(2*pthi_star-1))  

}

plot(data$delta,data$power)

#サンプルサイズ推定

alpha=0.05

data=data.frame(delta=0.01,n=seq(0,1000,1),power=0)

for(j in 1:nrow(data)){

lambda=data$n[j]*data$delta[j]  

C=(pthi+2*lambda)/(pthi+lambda)

pthi_star=(pthi+lambda)^2/(pthi+2*lambda)

W=qchisq(1-alpha,df=pthi-1)

data$power[j]=1-pnorm(sqrt((2*W)/C)-sqrt(2*pthi_star-1))  

}

plot(data$n,data$power)

#非心カイ二乗分布

n=100;lambda=seq(0.1,6,0.1);pthi=n-1

E_X=pthi+lambda

V_X=2*(pthi+2*lambda)

mean(rchisq(100000,df=n-1,ncp=4.5))

var(rchisq(100000,df=n-1,ncp=4.5))

params=data.frame(lambda=lambda,E_X=E_X,V_X=V_X)

#p.148 一元配置分散分析 分散が未知の場合:sigma0

n=20

X=t(matrix(c(rnorm(n,2,3),rnorm(n,3,4),rnorm(n,4,5),rnorm(n,8,9)),nrow=n))

S_A=sum((apply(X,1,mean)-mean(X))^2)*n

CT=sum(X)^2/prod(dim(X))

S_T=sum(X^2)-CT

S_E=S_T-S_A

pthi_A=nrow(X)-1;pthi_E=nrow(X)*(ncol(X)-1)

V_A=S_A/pthi_A;V_E=S_E/pthi_E

f0=V_A/V_E

alpha=0.05

qf(1-alpha,pthi_A,pthi_E)

a=apply(X,1,mean)-mean(X)

mu=mean(X)

#検出力曲線

data=data.frame(delta=seq(0.01,5,0.01),power=0)

for(j in 1:nrow(data)){

lambda=n*data$delta[j]  

C_A=(pthi_A+2*lambda)/(pthi_A+lambda)

pthi_star_A=(pthi_A+lambda)^2/(pthi_A+2*lambda)

W=qf(1-alpha,pthi_A,pthi_E)

data$power[j]=1-pnorm((sqrt(W/pthi_E)*sqrt(2*pthi_E-1)-sqrt(C_A/pthi_A)*sqrt(2*pthi_star_A-1))/sqrt(C_A/pthi_A+W/pthi_E))

}

plot(data$delta,data$power)

#サンプルサイズ推定

delta=1.1

data=data.frame(n=seq(2,100,1),power=0)

for(j in 1:nrow(data)){

lambda=data$n[j]*delta

pthi_E=nrow(X)*(data$n[j]-1)

C_A=(pthi_A+2*lambda)/(pthi_A+lambda)

pthi_star_A=(pthi_A+lambda)^2/(pthi_A+2*lambda)

W=qf(1-alpha,pthi_A,pthi_E)

data$power[j]=1-pnorm((sqrt(W/pthi_E)*sqrt(2*pthi_E-1)-sqrt(C_A/pthi_A)*sqrt(2*pthi_star_A-1))/sqrt(C_A/pthi_A+W/pthi_E))

}

plot(data$n,data$power)

#母不良率の検定

#H0:p=p0

p0=0.6

n=100

samp=sample(0:1,n,replace = T,prob=c(0.2,0.8))

p_hat=mean(samp)

u0=(p_hat-p0)/sqrt(p0*(1-p0)/n)

alpha=0.05

qnorm(1-alpha)

#検出力

P0=0.1;alpha=0.05

data=data.frame(P=seq(0,1,0.01),power=0)

for(j in 1:nrow(data)){

A=sqrt(P0*(1-P0)/(data$P[j]*(1-data$P[j])))

B=(data$P[j]-P0)/sqrt(data$P[j]*(1-data$P[j])/n)

data$power[j]=pnorm(-qnorm(1-alpha/2)*A-B)+1-pnorm(qnorm(1-alpha/2)*A-B)

}

plot(data$P,data$power,type="o")

#サンプルサイズ推定

P0=0.8;alpha=0.05;power=0.9

data=data.frame(delta=seq(0,P0,0.01),n=0,n1=0,n2=0) 

for(j in 1:nrow(data)){

P1=P0-data$delta[j]  

A=sqrt(P0*(1-P0)/(P1*(1-P1)))

B=(P1-P0)/sqrt(P1*(1-P1))

n1=((qnorm(power)+qnorm(1-alpha/2)*A)/(-B))^2

P2=P0+data$delta[j]  

A=sqrt(P0*(1-P0)/(P2*(1-P2)))

B=(P2-P0)/sqrt(P2*(1-P2))

n2=((qnorm(1-power)+qnorm(1-alpha/2)*A)/(-B))^2

data$n1[j]=n1;data$n2[j]=n2

if(P1<1 & P2<1){

data$n[j]=max(n1,n2)

}else{

if(P1<1){

data$n[j]=n1  

}else{

data$n[j]=n2  

}  

} 

}





#回帰係数の検定　p.170

y=c(2,3,5,7,9);x=c(2,3,6,8,9);n=5

Sxx=sum((x-mean(x))^2)

Syy=sum((y-mean(y))^2)

Sxy=sum((x-mean(x))*(y-mean(y)))

Se=Syy-Sxy^2/Sxx

V_e=S_e/pthi_e

pthi_e=n-2

beta10=1

t0=(cov(x,y)/var(x)-1)/sqrt(V_e/Sxx)

qt(1-alpha,pthi_e)

#検出力

alpha=0.05

data=data.frame(delta=seq(-5,5,0.01),power=0)

for(j in 1:nrow(data)){

lambda=sqrt(Sxx)*data$delta[j]

W=qt(1-alpha,pthi_e)

data$power[j]=pnorm((-W*(1-1/(4*pthi_e))-lambda)/sqrt(1+W^2/(2*pthi_e)))+1-pnorm((W*(1-1/(4*pthi_e))-lambda)/sqrt(1+W^2/(2*pthi_e)))


}

plot(data$delta,data$power)

#Sxxの推定

beta=0.1;alpha=0.05

data=data.frame(delta=seq(0.01,2,0.01),Sxx=0)

for(j in 1:nrow(data)){

data$Sxx[j]=((qnorm(1-alpha/2)-qnorm(beta))/data$delta[j])^2/(1-qnorm(1-alpha/2)^2/(2*pthi_e))  

}

plot(data$delta,data$Sxx)

#母不良率に関する1元配置分散分析の手順　p.174

n=100

P_hat=c(0.2,0.3,0.5,0.4,0.8,0.9)

k=length(P_hat)

sita_hat=asin(sqrt(P_hat))

mu=mean(sita_hat);a=sita_hat-mu

S_A=sum((sita_hat-mean(sita_hat))^2)

kai0=S_A*4*n

pthi_A=k-1

alpha=0.05

qchisq(1-alpha,pthi_A)

d1=max(sita_hat)-min(sita_hat)

d0=max(P_hat)-min(P_hat)

#検出力曲線

data=data.frame(delta=seq(0,5,0.01),power=0)

for(j in 1:nrow(data)){

lambda=n*data$delta[j]  

C=(pthi_A+2*lambda)/(pthi_A+lambda)

pthi_star=(pthi_A+lambda)^2/(pthi_A+2*lambda)

W=qchisq(1-alpha,pthi_A)

data$power[j]=1-pnorm(sqrt((2*W)/C)-sqrt(2*pthi_star-1))

}

plot(data$delta,data$power)

#検出力の最小値

k=4;alpha=0.05;d0=0.1;n=300

pthi_A=k-1

P_max=0.5+d0/2

P_min=0.5-d0/2

d1=asin(sqrt(P_max))-asin(sqrt(P_min))

lambda=2*n*d1^2

C=(pthi_A+2*lambda)/(pthi_A+lambda)

pthi_star=(pthi_A+lambda)^2/(pthi_A+2*lambda)

W=qchisq(1-alpha,pthi_A)

1-pnorm(sqrt((2*W)/C)-sqrt(2*pthi_star-1))

#サンプルサイズ推定

k=4;alpha=0.05;d0=0.1;

pthi_A=k-1

P_max=0.5+d0/2

P_min=0.5-d0/2

d1=asin(sqrt(P_max))-asin(sqrt(P_min))

delta=2*d1^2

data=data.frame(n=seq(2,1000,1),power=0)

for(j in 1:nrow(data)){

lambda=2*data$n[j]*d1^2

C=(pthi_A+2*lambda)/(pthi_A+lambda)

pthi_star=(pthi_A+lambda)^2/(pthi_A+2*lambda)

W=qchisq(1-alpha,pthi_A)

data$power[j]=1-pnorm(sqrt((2*W)/C)-sqrt(2*pthi_star-1))

}

plot(data$n,data$power)




#p.182 区間推定法に基づくサンプルサイズの設計法

#1つの母平均に基づくサンプルサイズの設計（分散既知）

n=100

X=rnorm(n,mean=2,sd=3)

alpha=0.05

under=mean(X)-qnorm(1-alpha/2)*sqrt(var(X)/n)

upper=mean(X)+qnorm(1-alpha/2)*sqrt(var(X)/n)

#必要サンプルサイズ

delta=1

n=(2*qnorm(1-alpha/2)*var(X))^2/(delta^2)

#1つの母分散に基づくサンプルサイズの設計

n=100

X=rnorm(n,mean=2,sd=3)

alpha=0.05;S=sum((X-mean(X))^2)

under=S/qchisq(1-alpha/2,n-1)

upper=S/qchisq(alpha/2,n-1)

#必要サンプルサイズ

delta=1.5

n=(((1+sqrt(delta))*qnorm(1-alpha/2))/(sqrt(delta)-1))^2/2+3/2

#1つの母平均の区間推定に基づくサンプルサイズの設計方法（分散未知）

n=100;alpha=0.05

X=rnorm(n,mean=2,sd=3)

mu=mean(X);S=sum((X-mean(X))^2);V=S/(n-1)

under=mu+qt(alpha/2,n-1)*sqrt(V/n)

upper=mu+qt(1-alpha/2,n-1)*sqrt(V/n)

pthi=n-1

C_star=ifelse(n>10,1,(sqrt(2)*gamma((n/2)))/(sqrt(n-1)*gamma((n-1)/2)))

sigma=2^2

n=(2*qt(alpha,n-1)*C_star*sqrt(sigma)/delta)^2

#12.4 2つの母分散の比の区間推定に基づくサンプルサイズの設計方法

n1=100;n2=200

X1=rnorm(n1,mean=2,sd=3);X2=rnorm(n2,mean=3,sd=5)

V1=sum((X1-mean(X1))^2)/(n1-1)

V2=sum((X2-mean(X2))^2)/(n2-1)

under=(V1/V2)/qf(1-alpha/2,n1-1,n2-1)

upper=(V1/V2)*qf(1-alpha/2,n2-1,n1-1)

#サンプルサイズ推定

data=data.frame(n=seq(2,100,1)) %>% mutate(delta=qf(1-alpha/2,n-1,n-1)^2)

# 2つの母平均の差の区間推定に基づくサンプルサイズ設計法

n1=100;n2=200

X1=rnorm(n1,mean=2,sd=3);X2=rnorm(n2,mean=3,sd=3)

X=c(X1,X2)

V=sum((X-mean(X))^2)/(n1+n2-1)

under=mean(X1)-mean(X2)+qt(alpha/2,n1+n2-2)*sqrt(V*(1/n1+1/n2))

upper=mean(X1)-mean(X2)+qt(1-alpha/2,n1+n2-2)*sqrt(V*(1/n1+1/n2))

#サンプルサイズ推定

delta=3;sigma=2

n=8*(sigma^2)*qnorm(1-alpha/2)^2/(delta^2)

#対応がある場合の母平均の差の区間推定に基づくサンプルサイズの設計方法

n=200

X1=rnorm(n,mean=2,sd=3);X2=rnorm(n,mean=3,sd=3)

d=X1-X2

V_d=sum((d-mean(d))^2)/(n-1)

alpha=0.05

under=mean(d)-qt(1-alpha,n-1)*sqrt(V_d/n)

upper=mean(d)+qt(1-alpha,n-1)*sqrt(V_d/n)

#サンプルサイズ

sigma=2;delta=3;alpha=0.05

n=(2*qnorm(1-alpha/2)*sigma)^2/(delta^2)
0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0