Newer
Older
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
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
# Define a server for the Shiny app
function(input, output) {
# Fill in the spot we created for a plot
output$Plot <- renderPlot({
myalpha1=input$alpha
mu0=input$mu0
sigma=input$sigma
n=input$n
mu1=input$mu1
criterion=qnorm(1-myalpha1)*sigma/sqrt(n)+mu0 # criterion c
betamu=pnorm((sqrt(n)/sigma)*(mu0-mu1)+qnorm(1-myalpha1)) #beta
power=1-pnorm((sqrt(n)/sigma)*(mu0-mu1)+qnorm(1-myalpha1)) # power
par(mar=c(4,0, 4, 0))
# distribut under null \bar{X}~N(mu0,sig/\sqrt{n})
plot(seq(-3,3,length=1000),dnorm(seq(-3,3,length=1000),mu0,sigma/sqrt(n)),type="l",
lwd=1,ylim=c(0,max(dnorm(seq(-3,3,length=1000),mu0,sigma/sqrt(n)))+0.2),cex.lab=2.2,xlab="",
main="",ylab="",cex=1.4, cex.axis=2,bty="n",
xlim=c(-5,6),yaxt="n",xaxt="n",col="black")
lines(seq(-3,5,length=1000),dnorm(seq(-3,5,length=1000),mu1,sigma/sqrt(n)),type="l",
lwd=1,ylim=c(0,1.5),cex.lab=2.2,xlab="",
main="",ylab="",cex=1.4, cex.axis=2,bty="n",
xlim=c(-5,7),yaxt="n",xaxt="n",col="black")
abline(v=criterion, lty=1, lwd=3, col="red")
axis(1, c(-5:5),cex.axis=2, tck=-.02)
mtext(expression(bar(X)[n]),at=c(0),side=1,line = 3,cex=2)
left <- seq(-3, criterion, length.out=100)
right <- seq(criterion, 3, length.out=100)
yH0r <- dnorm(right, mu0, sigma/sqrt(n))
yH1l <- dnorm(left, mu1, sigma/sqrt(n))
yH1r <- dnorm(right, mu1, sigma/sqrt(n))
polygon(c(right, rev(right)), c(yH0r, numeric(length(right))), border=NA,
density=5, lty=3, lwd=4, angle=-45,col="black")
polygon(c(left, rev(left)), c(yH1l, numeric(length(left))), border=NA,
density=5, lty=3, lwd=4, angle=45,col="black")
polygon(c(right, rev(right)), c(yH0r, numeric(length(right))), border=NA,
col=rgb(1, 0.3, 0.3, 0.6))
polygon(c(left, rev(left)), c(yH1l, numeric(length(left))), border=NA,
col=rgb(0.3, 0.3, 1, 0.6))
text(criterion-1.5, max(dnorm(seq(-3,3,length=1000),mu0,sigma/sqrt(n)))+0.2, adj=0, label=
as.expression(bquote(c~"="~.(format(criterion, digits = 3)))),cex=1.5)
text(-2.5, 0.3*max(dnorm(seq(-3,3,length=1000),mu0,sigma/sqrt(n)))
, adj=1, label= expression("Distrib. sous "*H[0]),cex=1.5)
text(mu1+4, 0.3*max(dnorm(seq(-3,13,length=1000),mu0,sigma/sqrt(n))), adj=1, label= expression("Distrib. sous "*H[a]),cex=1.5)
text(-2.5, 0.2*max(dnorm(seq(-3,3,length=1000),mu0,sigma/sqrt(n))), adj=1, label= expression("N(0,"*sigma/sqrt(n)*")"),cex=1.5)
text(mu1+4, 0.2*max(dnorm(seq(-3,3,length=1000),mu0,sigma/sqrt(n))), adj=1, label= expression("N(1,"*sigma/sqrt(n)*")"),cex=1.5)
source("http://www.math.mcmaster.ca/bolker/R/misc/legendx.R")
legend("topright",c(
as.expression(bquote(alpha[n](mu)~"="~.(format(myalpha1, digits = 3)))),
as.expression(bquote(beta[n](mu)~"="~.(format(betamu,digits = 3))))),
text.col="white",
col = c(NA,NA,NA),lty=c(NA,NA,1),
border = c("#FF4D4D99","#4D4DFF99",NA),
fill=c(rgb(1, 0.3, 0.3, 0.6),rgb(0.3, 0.3, 1, 0.6),NA),
cex=1.5,bty = "n",angle=c(-45,45,90),
box.cex=c(2,2),y.intersp=3)
legend("topright",c(
as.expression(bquote(alpha[n](mu)~"="~.(format(myalpha1, digits = 3)))),
as.expression(bquote(beta[n](mu)~"="~.(format(betamu,digits = 3))))),
text.col="white",
col = c(NA, NA,NA),lty=c(2,2,NA),
border = c("#FF4D4D99","#4D4DFF99",NA),
fill=c(rgb(1, 0.3, 0.3, 0.6),rgb(0.3, 0.3, 1, 0.6),NA),
cex=1.5,density=c(10,10,0),bty = "n",angle=c(-45,45,90),
box.cex=c(2,2),y.intersp=3)
legend("topright",c(
as.expression(bquote(alpha[n](mu)~"="~.(format(myalpha1, digits = 3)))),
as.expression(bquote(beta[n](mu)~"="~.(format(betamu,digits = 3))))),
text.col="black",
col = c(NA, NA,NA),lty=c(3,3,NA),
border = c("#FF4D4D99","#4D4DFF99",NA),
cex=1.5,density=c(10,10,0),bty = "n",angle=c(-45,45,90),
box.cex=c(2,2),y.intersp=3)
})
}