假设检验-p.value-对汤银才此程序的疑问

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

对汤老师的p-value有点小小的疑惑？R语言与统计分析-by汤银才-p171
此题是方差sigma已知的情况下， 对假设进行检验，并对均值区间进行估计
z.test<-function(x,n,sigma,alpha,u0=0,alternative="two.sided"){
                        options(digits=4)
                        result<-list( )
                        mean<-mean(x)
                        z<-(mean-u0)/(sigma/sqrt(n))
                        p<-pnorm(z,lower.tail=FALSE)
                        result$mean<-mean
                        result$z<-z
                        result$p.value<-p
                        if(alternative=="two.sided"){
                         p<-2*p
                         result$p.value<-p
                       }
                         else if (alternative == "greater"|alternative =="less" ){
                         result$p.value<-p
                       }
                         else return("your input is wrong")
                   result$conf.int<- c(mean-sigma*qnorm(1-alpha/2,mean=0, sd=1,lower.tail = TRUE)/sqrt(n),
                                                 mean+sigma*qnorm(1-alpha/2,mean=0, sd=1,lower.tail = TRUE)/sqrt(n))
result
}
红色部分，是疑惑的地方，p.value是取得 lower.tail=F，是上侧的面积，那么 less时即假设H0：u>= u0时，p.value =p是不是错误的呢？另外，双侧检验的时候，z如果是负数，那么2*p(上侧面积）会不会出现 大于 1的情况！！！！

也尝试修改下，大大们帮忙看看有没有统计知识的问题撒？
z.test <- function(x,n,sigma,alpha,u0=0,alternative="two.sided"){
  options(digits=4)
  result <- list()
  mean <- mean(x)
  z <-(mean-u0)/(sigma/sqrt(n))
   result$mean <- mean
  result$z <- z
  ifelse(alternative=="two.sided",result$p.value <- 2*(1-pnorm(abs(z))),      #双侧检测
         ifelse(alternative=="less",result$p.value <- pnorm(z,lower.tail=T),           #下侧检验
                ifelse(alternative=="greater",result$p.value <-pnorm(z,lower.tail=F),   #上侧检验
                       return("your input is wrong"))))
  result$conf.int <- c(mean-qnorm(1-alpha/2)*sigma/sqrt(n),
                       mean+qnorm(1-alpha/2)*sigma/sqrt(n)) #有误
  result
}  

刚学到 假设检验， p-value 又是太重要了
自学统计， 没老师。。。
请大大们帮忙指点统计知识的问题撒 谢谢

谢谢 epoh老师的指导 之前把alternative 给搞反了，p.value 的面积真是有点绕

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏2 回帖

关键词：value 假设检验 alue 汤银才 p-value 程序 检验

Z test for known population standard deviation
z.test(x, mu = 0, stdev, alternative = c("two.sided", "less", "greater"), sd = stdev, conf.level = 0.95, ...)

   http://rss.acs.unt.edu/Rdoc/libr ... os/html/z.test.html

    z.test<-
function(x, mu=0, stdev,
                   alternative = c("two.sided", "less", "greater"),
                   sd=stdev,
                   conf.level = 0.95, ... ){

  if(missing(stdev) && missing(sd)) stop("You must specify a Standard Deviation of the population")

  alternative <- match.arg(alternative)

  n <- length(x)
  z <- (mean(x)-mu)/(sd/sqrt(n))

  out <- list(statistic=c(z=z))
  class(out) <- 'htest'

  out$parameter <- c(n=n,"Std. Dev." = sd,
                     "Std. Dev. of the sample mean" = sd/sqrt(n))

  out$p.value <- switch(alternative,
                    two.sided = 2*pnorm(abs(z),lower.tail=FALSE),
                    less = pnorm(z),
                    greater = pnorm(z, lower.tail=FALSE) )

  out$conf.int <- switch(alternative,
                         two.sided = mean(x) +
                           c(-1,1)*qnorm(1-(1-conf.level)/2)*sd/sqrt(n),
                         less = c(-Inf, mean(x)+qnorm(conf.level)*sd/sqrt(n)),
                         greater = c(mean(x)-qnorm(conf.level)*sd/sqrt(n), Inf)
                         )
  attr(out$conf.int, "conf.level") <- conf.level

  out$estimate <- c("mean of x" = mean(x))
  out$null.value <- c("mean" = mu)
  out$alternative <- alternative
  out$method <- "One Sample z-test"
  out$data.name <- deparse(substitute(x))
  names(out$estimate) <- paste("mean of", out$data.name)

  return(out)
}

epoh 发表于 2012-12-3 18:12
Z test for known population standard deviation
z.test(x, mu = 0, stdev, alternative = c("two.sided" ...

z.test(y$Weight,mu=3,stdev=0.18,"less",conf.level=0.99)        
One Sample z-test
data:  y$Weight
z = -2.6667, n = 36.00, Std. Dev. = 0.18, Std. Dev. of thesample mean = 0.03, p-value = 0.00383
alternative hypothesis: true mean is less than 3
99 percent confidence interval:   
-Inf 2.98979
sample estimates:
mean of y$Weight            
2.92
epoh 老师，像是这个结果输出格式，是怎么控制的？是 class(out) <- "htest"的原因吗，特别是出现
alternative hypothesis: true mean is less than 3
原代码好像是没有这个的？

另外，t.test的代码怎么查看呢？
t.testfunction (x, ...)
UseMethod("t.test")
<bytecode: 0x0ea3488c>
<environment: namespace:stats>
> methods(t.test)
[1] t.test.default* t.test.formula*   
Non-visible functions are asterisked
而在 z.test中detail，说Once the student has learned about the T distribution the t.test function should be used instead of this one (but the syntax is very similar, so this function should be an appropriate introductory step to learning t.test).

谢谢 epoh老师的指导

jmpamao 发表于 2012-12-4 22:50
epoh 发表于 2012-12-3 18:12
Z test for known population standard deviation
z.test(x, mu = 0, stde ...

哈哈,你的水平很好,一说就对.
的确是来自class htest (hypothesis test)

print.htest <- function(x, digits = 4, quote = TRUE, prefix = "", ...)
{
    cat("\n")
    cat(strwrap(x$method, prefix = "\t"), sep="\n")
    cat("\n")
    cat("data: ", x$data.name, "\n")
    out <- character()
    if(!is.null(x$statistic))
        out <- c(out, paste(names(x$statistic), "=",
                            format(round(x$statistic, 4))))
    if(!is.null(x$parameter))
        out <- c(out, paste(names(x$parameter), "=",
                            format(round(x$parameter, 3))))
    if(!is.null(x$p.value)) {
        fp <- format.pval(x$p.value, digits = digits)
        out <- c(out, paste("p-value",
                            if(substr(fp, 1L, 1L) == "<") fp else paste("=",fp)))
    }
    cat(strwrap(paste(out, collapse = ", ")), sep="\n")
    if(!is.null(x$alternative)) {
        cat("alternative hypothesis: ")
        if(!is.null(x$null.value)) {
            if(length(x$null.value) == 1L) {
                alt.char <-
                    switch(x$alternative,
                           two.sided = "not equal to",
                           less = "less than",
                           greater = "greater than")
                cat("true", names(x$null.value), "is", alt.char,
                    x$null.value, "\n")
            }
            else {
                cat(x$alternative, "\nnull values:\n")
                print(x$null.value, ...)
            }
        }
        else cat(x$alternative, "\n")
    }
    if(!is.null(x$conf.int)) {
        cat(format(100 * attr(x$conf.int, "conf.level")),
            "percent confidence interval:\n",
            format(c(x$conf.int[1L], x$conf.int[2L])), "\n")
    }
    if(!is.null(x$estimate)) {
        cat("sample estimates:\n")
        print(x$estimate, ...)
    }
    cat("\n")
    invisible(x)
}


.../ src/library/stats/R/htest.R
.../src/library/stats/R/t.test.R

t.test.rar (2 KB) 本附件包括：

• t.test.R

2012-12-5 11:01:12 上传

epoh 发表于 2012-12-5 10:59
哈哈,你的水平很好,一说就对.
的确是来自class htest (hypothesis test)

epoh老师 再向您请教个问题，

print.myout <- function(x){print("it's very good")}
out <- list()
class(out) <- "myout"
out
>"it's very good"
print.myout 是怎么识别myout的？或者说 out 怎么知道给他一个"myout"字符， 就知道指的是 print.myout这个function。

另外 能否 简单介绍下  mode  class attr 的关系与区别吗？ 谢谢 epoh 老师

jmpamao 发表于 2012-12-8 23:45
epoh老师 再向您请教个问题，

print.myout

看了底下例子,也许更容易明白.
例子来自The Art of R Programming
Chapter 12 Object-Oriented Programming
page 97/193

One can apply print()(Generic Functions),to many types of objects
methods(print)
  [1] print.acf*                                 
  [2] print.anova                                
  [3] print.aov*                                 
  [4] print.aovlist*                             
  [5] print.ar*     
...
...
[168] print.warnings                             
[169] print.xgettext*                           
[170] print.xngettext*                           
[171] print.xtabs*      

#create a object of class "employee"
j <- list(name="Joe", salary=55000, union=T)
class(j) <- "employee"
attributes(j)
print(j)
$name
[1] "Joe"

$salary
[1] 55000

$union
[1] TRUE

attr(,"class")
[1] "employee"

#now write print function for object employee
print.employee <- function(wrkr) {
cat(wrkr$name,"\n")
cat("salary",wrkr$salary,"\n")
cat("union member",wrkr$union,"\n")
}

print(j)
Joe
salary 55000
union member TRUE

######
methods(print)
  [1] print.acf*                                 
  [2] print.anova                                
  [3] print.aov*                                 
  [4] print.aovlist*                             
  [5] print.ar*     
...
...
[65]print.employee  ...
...
[168] print.vignette*                           
[169] print.warnings                             
[170] print.xgettext*                           
[171] print.xngettext*                           
[172] print.xtabs*