【学习记录-R】自编函数进行单因素方差分析（等容或不等容）

一.每个水平等容 

原理：

 该图有错误，G/F的自由度写成了（dft,dfT)，应该改为（dft,dfe)

步骤：

计算meanx,meanj,sse,ssb,dfe,dft,dfT(样本总均值，每个水平样本均值，误差项平方和，水平项平方和，误差项平方和的自由度，水平项平方和自由度)

计算fx样本得来的观测值，(ssb/dft)/(sse/dfe)

计算p值，pf(fx,df1,dft,dfe)

难点：

 输入值为长数据，怎么将其分开计算每个水平样本均值：二重for循环

使用跑得慢但写的快的苏联代码(bushi

等容：

meanx<-function(x,k,n){
  xk<-c(0);z1=1
  for(i in 1:k){
    m=0#重新计算第i个水平的样本总和m
    z2=z1
    z3=z2+n-1#第i个水平共有n个样本，目前的第z2:z2+n-1个样本为第i行的样本
    for(j in z2:z3){
      m= m+x[j]#计算第i个水平的样本总和m
      z1=z1+1
      #cat("\nz2=",z2,",j=",j,",xj=",x[j],",m=",m)
    }
    xk[i]=m/n#计算第i个水平的样本均值
  }
  return(xk[])
}

不等容：

 

meanx<-function(x,y,k){
  xk<-c(0);z1=1
  for(i in 1:k){
    m=0#重新计算第i个水平的样本总和m
    z2=z1
    z3=z2+y[i]-1
    for(j in z2:z3){#第i个水平共有y[i]个样本，目前的第z2:z2+y[i]-1个样本为第i行的样本
      m= m+x[j]#计算第i个水平的样本总和m
      z1=z1+1
      #cat("\nz2=",z2,",j=",j,",xj=",x[j],",m=",m)
    }
    xk[i]=m/y[i]#计算第i个水平的样本均值
  }
  return(xk[])
}

fcfx<-function(a,x,k,n){
  xk<-meanx(x,k,n)
  meanx=sum(x)/(k*n)
  sse<-sse(x,k,n,xk)
  ssb=n*sum((xk-meanx)^2)
  dft=k-1
  dfT=k*n-1
  dfe=k*(n-1)
  fx=(ssb/dft)/(sse/dfe)
  f=qf(1-a,dft,dfe)
  p<-1-pf(fx,dft,dfe)
  list(dft=dft,dfT=dfT,dfe=dfe,sse=sse,ssb=ssb,fx=fx,f=f,p=p)
}

sse<-function(x,k,n,xk){
  z1=1;sse=0
  for(i in 1:k){
    z2=z1
    z3=z2+n-1
    for(j in z2:z3){
      sse=sse+(x[j]-xk[i])^2
      z1=z1+1
     #cat("\nz2=",z2,",j=",j,",xj=",x[j],",sum=",sum,",xk=",xk[i])
    }
  }
  return(sse)
}

meanx<-function(x,k,n){
  xk<-c(0);z1=1
  for(i in 1:k){
    m=0
    z2=z1
    z3=z2+n-1
    for(j in z2:z3){
      m= m+x[j]
      z1=z1+1
      #cat("\nz2=",z2,",j=",j,",xj=",x[j],",m=",m)
    }
    xk[i]=m/n
  }
  return(xk[])
}

 

例题：

 

 

 

二.不等容 

 原理：

  该图有错误，G/F的自由度写成了（dft,dfT)，应该改为（dft,dfe),

代码：

fcfx2<-function(a,x,y,k){#样本不等容量
  xk<-meanx(x,y,k)
  sse<-sse(x,y,k,xk)
  sumn<-sum(y)
  meanx=sum(x)/sumn
  ssb=sum(y*(xk-meanx)^2)
  dft=k-1
  dfT=sumn-1
  dfe=sumn-k
  fx=(ssb/dft)/(sse/dfe)
  f=qf(1-a,dft,dfe)
  p<-1-pf(fx,dft,dfe)
  list(dft=dft,dfT=dfT,dfe=dfe,sse=sse,ssb=ssb,fx=fx,f=f,p=p)
}

sse<-function(x,y,k,xk){
  z1=1;sse=0
  for(i in 1:k){
    z2=z1
    z3=z2+y[i]-1
    for(j in z2:z3){
      sse=sse+(x[j]-xk[i])^2
      z1=z1+1
      #cat("\nz2=",z2,",j=",j,",xj=",x[j],",sum=",sum,",xk=",xk[i])
    }
  }
  return(sse)
}

meanx<-function(x,y,k){
  xk<-c(0);z1=1
  for(i in 1:k){
    m=0
    z2=z1
    z3=z2+y[i]-1
    for(j in z2:z3){
      m= m+x[j]
      z1=z1+1
      #cat("\nz2=",z2,",j=",j,",xj=",x[j],",m=",m)
    }
    xk[i]=m/y[i]
  }
  return(xk[])
}

 

例题：

1.不等容

 

 

2.等容

题目同上等容例题