实验报告六答案

##
## Two Sample t-test
##
## data: x and y
## t = 4.684, df = 17, p-value = 0.0002131
## alternative hypothesis: true difference in means is not
equal to 0
n <- 50 m <- 500 mu0 <- 20 sigma <- 5 mu <- c(seq(15, 25, 0.1)) #alternatives M <- length(mu) power2 <- numeric(M) for (i in 1:M) {
mu2 <- mu[i] pvalues2 <- replicate(m, expr = {
# simulate under alternative mu1 x <- rnorm(n, mean = mu2, sd = sigma) ttest <- t.test(x, alternative = "two.sided", mu = mu0) ttest$p.value }) power2[i] <- mean(pvalues2 <= 0.05)
n <- length(x) xb <- mean(x) if (sigma >= 0)
tmp <- sigma/sqrt(n) * qnorm(1 - alpha/2) else tmp < - sd(x)/sqrt(n) * qnorm(1 - alpha/2)
data.frame(mean = xb, a = xb - tmp, b = xb + tmp) } x <- c(rep(0, 7), rep(1, 10), rep(2, 12), rep(3, 8), rep(4, 3), rep(5, 2), rep(6,
bg = "gray")
file://I:\统计分析软件应用\实验报告及答案\6.html
2014/6/16
w
页码，2/4(W)
教材4.4
objf <- function(x) {
f <- c(-13 + x[1] + ((5 - x[2]) * x[2] - 2) * x[2], -29 + x[1] + ((x[2] +
1) * x[2] - 14) * x[2]) sum(f^2) } x0 <- c(0.5, -2) nlm(objf, x0)
## $minimum ## [1] 48.98 ## ## $estimate ## [1] 11.4128 -0.8968 ## ## $gradient ## [1] 1.415e-08 -1.435e-07 ## ## $code ## [1] 1 ## ## $iterations ## [1] 16
mu1 <- mu[i] pvalues <- replicate(m, expr = {
# simulate under alternative mu1 x <- rnorm(n, mean = mu1, sd = sigma) ttest <- t.test(x, alternative = "two.sided", mu = mu0) ttest$p.value }) power[i] <- mean(pvalues <= 0.05) }
最优解为（11.4128 ，-0.8968）
教材习题4.6
file://I:\统计分析软件应用\实验报告及答案\6.html
2014/6/16
w
页码，3/4(W)
x <- c(140, 137, 136, 140, 145, 148, 140.135, 144, 141) y <- c(135, 118, 115, 140, 128, 131, 130, 115, 131, 125) t.test(x, y, var.equal = TRUE)
} plot(mu, power, type = "l", col = "red") lines(mu, power2, col = "blue") legend(22, 0.8, lty = c(1, 1), c("n=100", "n=50"), col = c ("red", "blue"), text.col = "green",
0)) interval_estimate3(x)
## mean
a
b
## 1 1.905 1.494 2.315
mean(x)
## [1] 1.905
file://I:\统计分析软件应用\实验报告及答案\6.html
2014/6/16
从上面的结果知：总体均值差的95%置信区间为 （7.935 20.94）。教材习题4.9 用教材
190页的命令：实际上就是利用中心极限定理写出的近似区间。（非正态总体的均值的置 信区间）
file://I:\统计分析软件应用\实验报告及答案\6.html
2014/6/16
w
页码，4/4(W)
interval_estimate3 <- function(x, sigma = -1, alpha = 0.05) {
## 95 percent confidence interval:
## 7.935 20.939
## sample estimates:
## mean of x mean o来自 y##141.2
126.8
或者用教材中的命令interval_estimate2
interval_estimate2 <- function(x, y, sigma = c(-1, -1), var.equal = FALSE, alpha = 0.05) {
} } data.frame(mean = xb - yb, df = df, a = xb - yb - tmp, b = xb - yb + tmp) } interval_estimate2(x, y, var.equal = TRUE)
## mean df
a
b
## 1 14.44 17 7.935 20.94
w
页码，1/4(W)
Title
set.seed(100) n <- 100 m <- 500 mu0 <- 20 sigma <- 5 mu <- c(seq(15, 25, 0.1)) #alternatives M <- length(mu) power <- numeric(M) for (i in 1:M) {
} else { S1 <- var(x) S2 <- var(y) nu <- (S1/n1 + S2/n2)^2/(S1^2/n1^2/(n1 - 1) +
S2^2/n2^2/(n2 - 1)) tmp <- qt(1 - alpha/2, nu) * sqrt(S1/n1 + S2/n2) df <- nu
if (var.equal == TRUE) { Sw <- ((n1 - 1) * var(x) + (n2 - 1) * var(y))/
(n1 + n2 - 2) tmp <- sqrt(Sw * (1/n1 + 1/n2)) * qt(1 -
alpha/2, n1 + n2 - 2) df <- n1 + n2 - 2
n1 <- length(x) n2 <- length(y) xb <- mean(x) yb <- mean(y) if (all(sigma >= 0)) {
tmp <- qnorm(1 - alpha/2) * sqrt(sigma[1]^2/n1 + sigma[2]^2/n2)
df <- n1 + n2 } else {