借助分布可加性验证大样本样本均值的渐进分布

一、大数定律

#函数四：大数定律的验证
###########大数定律,生成的数叫y2，要加和，直接生成叫y1
LLN <- function(m,n,fun='Normal',...){
  #m是画图的抽样样本数(列向量)，n是加和样本数-大样本(按行生成)
  #Y1根据可加性规律直接生成m个样本，并标准化
  #rf按照定义，生成m行n列个子样本
  #sumrf按行求和得到m个所需要的样本，
  #Y2 sumrf的标准化
  
  
  args <- list(...)
  ####################################
  if(fun=='Binomial'){
    ######参数说明dbinom(x, size, prob, log = FALSE)  
    #args[[1]]=p 伯努利分布EX=p,DX=pq
    
    #Y1 <- (rbinom(m,n,args[[1]])-n*args[[1]])/sqrt(n*args[[1]]*(1-args[[1]]))  
    X <- matrix(rbinom(m*n, size=1, prob=args[[2]]), nrow = m,ncol = n) #n行m列矩阵，服从b(1,p)
    
    g1 <- function(x){args[[2]]+3*sqrt(args[[2]]*(1-args[[2]])/x)}
    g2 <- function(x){args[[2]]-3*sqrt(args[[2]]*(1-args[[2]])/x)}
    miu <- args[[2]]
    
    
    
  }else if(fun=='Poisson'){
    ######参数说明dpois(x, lambda, log = FALSE)
    #args[[1]]=lamda 泊松分布 EX=lamda,DX=lamda
    
    #Y1 <- (rpois(m,n*args[[1]])-n*args[[1]])/sqrt(n*args[[1]])
    X <- matrix(rpois(m*n, lambda = args[[1]]), nrow = m,ncol = n) #m行n列矩阵，服从b(1,p)
    
    g1 <- function(x){args[[1]]+3*sqrt(args[[1]]/x)}
    g2 <- function(x){args[[1]]-3*sqrt(args[[1]]/x)}
    miu <- args[[1]]
    
  }else if(fun=='Normal'){
    #####参数说明dnorm(x, mean = 0, sd = 1, log = FALSE)
    #args[[1]]=miu,ard[[2]]=sigma^2 正态分布 EX=miu,DX=sigma^2
    
    #Y1 <- (rnorm(m,n*args[[1]],sqrt(n)*args[[2]])-n*args[[1]])/(sqrt(n)*args[[2]])
    X <- matrix(rnorm(m*n, mean = args[[1]], sd = args[[2]]), nrow = m,ncol = n) #n行m列矩阵，服从b(1,p)
    
    g1 <- function(x){args[[1]]+3*args[[2]]/sqrt(x)}
    g2 <- function(x){args[[1]]-3*args[[2]]/sqrt(x)}
    miu <- args[[1]]
    
    
  }else if(fun=='Gamma'){
    #####参数说明dgamma(x, shape, rate = 1, scale = 1/rate, log = FALSE)
    #args[[1]]=alpha,args[[2]]=lamda gamma分布 EX=alpha/lamda,DX=alpha^2/lamda
    
    #Y1 <- (rgamma(m,n*args[[1]],args[[2]])-n*args[[1]]/args[[2]])/sqrt(n*args[[1]]/(args[[2]]^2))
    X <- matrix(rgamma(m*n, shape = args[[1]], rate = args[[2]]), nrow = m,ncol = n) #n行m列矩阵，服从b(1,p)
    
    g1 <- function(x){args[[1]]/args[[2]]+3*sqrt(args[[1]]/x)/args[[2]]}
    g2 <- function(x){args[[1]]/args[[2]]-3*sqrt(args[[1]]/x)/args[[2]]}
    miu <- args[[1]]/args[[2]]
    
  }else if(fun == 'Chi-Aquared'){
    #####参数说明dchisq(x, df, ncp = 0, log = FALSE)
    #args[[1]]='n',卡方分布 EX='n',DX = '2n'
    
    #Y1 <- (rchisq(m,n*args[[1]])-n*args[[1]])/sqrt(n*2*args[[1]])
    
    X <- matrix(rchisq(m*n, df = args[[1]]), nrow = m,ncol = n) #m行n列矩阵，服从b(1,p)
    
    g1 <- function(x){args[[1]]+3*sqrt(2*args[[1]]/x)}
    g2 <- function(x){args[[1]]-3*sqrt(2*args[[1]]/x)}
    miu <- args[[1]]
  }
  
  Y3 <- apply(X, 1, function(z)(cumsum(z)/seq_along(z)))
  Y3 <- matrix(Y3,nrow = m*n,ncol = 1)
  Y <- as.data.frame(cbind(rep(1:n,m),rep(1:m,each = n),Y3))
  colnames(Y) <- c('num','sig','y3')
  
  #画图
  ggplot(Y) + 
    geom_line(
      aes(x = num,
          y = y3,
          colour=sig,
          group=sig
      )
    ) +
    stat_function(
      data = data.frame(x = 1:n, Legend = factor(1)),
      fun = g1,col = 'red',size = 1,linetype = 2
    )+
    stat_function(
      data = data.frame(x = 1:n, Legend = factor(2)),
      fun = g2,col = 'red',size = 1,linetype = 2
    )+
    geom_hline(
      yintercept = miu,
      #Legend = factor(1),
      color = 'black',size = 1.5,linetype = 3
    )+
    labs(
      title = 'Large Numbers Theorem ',
      subtitle = fun,
      x = 'n', y = 'rf'
    ) +
    theme_light()+
    theme(legend.position = 'none')
    
  
  
}

• 结果

LLN(100,1000,fun = 'Binomial',size = 1,prob = 0.3)

LLN(100,1000,fun = 'Poisson',lambda = 2)

LLN(100,1000,fun = 'Normal',mean = 5, sd = 1.5)#正态分布标准化，n=1就可

LLN(100,1000,fun = 'Gamma',shape = 2, rate = 3)

LLN(100,1000,fun = 'Chi-Aquared',df = 3.5)

二、概率分布的可加性

#函数一：可加性验证
Add <- function(n = 2,m = 20000,fun='Normal',...){
  #m是画图的抽样样本数，n是加和样本数-大样本
  args <- list(...)
  if(fun=='Binomial'){
    ######参数说明dbinom(x, size, prob, log = FALSE)
    rf <- matrix(rbinom(m*n,args[[1]],args[[2]]),nrow = m,ncol = n,byrow = TRUE)
    sumrf <- rowSums(rf)
    x <- seq(500,1200,by = 1)
    y <- dbinom(x,sum(args[[1]]),args[[2]][1])

  }else if(fun=='Poisson'){
    ######参数说明dpois(x, lambda, log = FALSE)
    rf <- matrix(rpois(m*n,args[[1]]),nrow = m,ncol = n,byrow = TRUE)
    sumrf <- rowSums(rf)
    x <- seq(0,50,by = 1)#必须是整数
    y <- dpois(x,sum(args[[1]]))

  }else if(fun=='Normal'){
    #####参数说明dnorm(x, mean = 0, sd = 1, log = FALSE)
    rf <- matrix(rnorm(m*n,args[[1]],args[[2]]),nrow = m,ncol = n,byrow = TRUE)
    sumrf <- rowSums(rf)
    x <- seq(-50,50,by = 0.1)
    y <- dnorm(x,sum(args[[1]]),sqrt(sum(args[[2]]^2)))

  }else if(fun=='Gamma'){
    #####参数说明dgamma(x, shape, rate = 1, scale = 1/rate, log = FALSE)
    rf <- matrix(rgamma(m*n,args[[1]],args[[2]]),nrow = m,ncol = n,byrow = TRUE)
    sumrf <- rowSums(rf)
    x <- seq(0,20,by = 0.1)
    y <- dgamma(x,sum(args[[1]]),args[[2]][1])

  }else if(fun == 'Chi-Aquared'){
    #####参数说明dchisq(x, df, ncp = 0, log = FALSE)
    rf <- matrix(rchisq(m*n,args[[1]]),nrow = m,ncol = n,byrow = TRUE)
    sumrf <- rowSums(rf)
    x <- seq(0,50,by = 0.1)
    y <- dchisq(x,sum(args[[1]]))
    
  }
  
  #画图
  ggplot() +
    geom_histogram(#绘制抽得的随机数直方图
      aes(x = sumrf, y = ..density..),
      bins = 30,
      fill = 'cornflowerblue',
      color = 'darkgrey'
    ) +
    geom_line(#绘制抽样密度曲线【按行加和的】
      aes(x = sumrf,color = 'red'),
      stat = 'density',
      size = 1.2
    ) +
    geom_line(#绘制理论密度曲线【直接加和的版本】
      aes(x = x,y = y,color = 'green'),
      size = 1.2
    ) +
    scale_colour_manual(
      name = 'Legend',
      values = c('red', 'green'),
      labels = c( "Theoretical probability density","Sampling probability density")
    ) +
    labs(#加题目,改x,y轴
      title = 'Additivity of the density function',
      subtitle = fun,
      x = 'X-axis',
      y = ' Density'
    ) +
    theme_light() +
    theme(legend.position = c(0.85,0.9))
  

}

• 结果

Add(n=3,m=10000,fun='Binomial',Size = c(1000,2000,1500), Prob = 0.2)#Prob参数一样

Add(n=3,m=10000,fun='Poisson',Lamda = c(1,2,3))

Add(n=3,m=10000,fun='Normal',Mean = c(1,3,5),Sd = c(5,6,2))

Add(n=3,m=10000,fun='Gamma',Shape = c(5,4,6), rate = 7) #Rate参数一样

Add(n=3,m=10000,fun='Chi-Aquared',Df = c(10,15,4))

三、中心极限定理

#函数二：三条比较
Comp <- function(m,n,fun='Normal',...){
  #m是画图的抽样样本数(列向量)，n是加和样本数-大样本(按行生成)
  #Y1根据可加性规律直接生成m个样本，并标准化
  #rf按照定义，生成m行n列个子样本
  #sumrf按行求和得到m个所需要的样本，
  #Y2 sumrf的标准化
  
  timestart<-Sys.time()#比较两种生成方法的
       
  
  args <- list(...)
  ####################################
  if(fun=='Binomial'){
    ######参数说明dbinom(x, size, prob, log = FALSE)  
    #args[[1]]=p 伯努利分布EX=p,DX=pq
    
    timestart1<-Sys.time()
    Y1 <- (rbinom(m,n,args[[2]])-n*args[[2]])/sqrt(n*args[[2]]*(1-args[[2]]))
    timeend1<-Sys.time()
    
    timestart2<-Sys.time()
    rf <- matrix(rbinom(m*n,1,args[[2]]),nrow = m,ncol = n)
    sumrf <- rowSums(rf) 
    Y2 <- (sumrf-n*args[[2]])/sqrt(n*args[[2]]*(1-args[[2]]))
    timeend2<-Sys.time()
    

    
  }else if(fun=='Poisson'){
    ######参数说明dpois(x, lambda, log = FALSE)
    #args[[1]]=lamda 泊松分布 EX=lamda,DX=lamda
    
    timestart1<-Sys.time()    
    Y1 <- (rpois(m,n*args[[1]])-n*args[[1]])/sqrt(n*args[[1]])
    timeend1<-Sys.time()
    
    timestart2<-Sys.time()
    rf <- matrix(rpois(m*n,args[[1]]),nrow = m,ncol = n)
    sumrf <- rowSums(rf)
    Y2 <- (sumrf-n*args[[1]])/sqrt(n*args[[1]])
    timeend2<-Sys.time()
    
  }else if(fun=='Normal'){
    #####参数说明dnorm(x, mean = 0, sd = 1, log = FALSE)
    
    timestart1<-Sys.time()    
    Y1 <- (rnorm(m,n*args[[1]],sqrt(n)*args[[2]])-n*args[[1]])/(sqrt(n)*args[[2]])
    timeend1<-Sys.time()
    
    timestart2<-Sys.time()
    rf <- matrix(rnorm(m*n,args[[1]],args[[2]]),nrow = m,ncol = n)
    sumrf <- rowSums(rf)
    Y2 <- (sumrf-n*args[[1]])/(sqrt(n)*args[[2]])
    timeend2<-Sys.time()
    
  }else if(fun=='Gamma'){
    #####参数说明dgamma(x, shape, rate = 1, scale = 1/rate, log = FALSE)
    
    timestart1<-Sys.time()    
    Y1 <- (rgamma(m,n*args[[1]],args[[2]])-n*args[[1]]/args[[2]])/sqrt(n*args[[1]]/(args[[2]]^2))
    timeend1<-Sys.time()
    
    timestart2<-Sys.time()
    rf <- matrix(rgamma(m*n,args[[1]],args[[2]]),nrow = m,ncol = n)
    sumrf <- rowSums(rf)
    Y2 <- (sumrf-n*args[[1]]/args[[2]])/sqrt(n*args[[1]]/(args[[2]]^2))
    timeend2<-Sys.time()
    
  }else if(fun == 'Chi-Aquared'){
    #####参数说明dchisq(x, df, ncp = 0, log = FALSE)
    
    timestart1<-Sys.time()    
    Y1 <- (rchisq(m,n*args[[1]])-n*args[[1]])/sqrt(n*2*args[[1]])
    timeend1<-Sys.time()
    
    timestart2<-Sys.time()
    rf <- matrix(rchisq(m*n,args[[1]]),nrow = m,ncol = n)
    sumrf <- rowSums(rf)
    Y2 <- (sumrf-n*args[[1]])/sqrt(n*2*args[[1]]) 
    timeend2<-Sys.time()
  }
  
  
  #画图
  x = seq(-7,7,by = 0.1)
  Plot <- ggplot() +
    geom_line(#绘制可加性抽样
      aes(x = Y1,color = 'red'),
      stat = 'density',
      size = 1.2
    ) +
    geom_line(#绘制单独抽样
      aes(x = Y2,color = 'green'),
      stat = 'density',
      size = 1.2
    ) +
    geom_line(#绘制理论密度曲线N(0,1)
      aes(x = x,y = dnorm(x),color = 'blue'),
      size = 1.2
    ) +
    scale_colour_manual(
      name = 'Legend',
      values = c('red', 'green','blue'),
      labels = c( "Additively sampling probability density",
                  "Respectively sampling probability density",
                  "Theoretical probability density")
    ) +
    labs(#加题目,改x,y轴
      title = 'Three probability distribution curves comparation',
      subtitle = fun,
      x = 'X-axis',
      y = ' Density'
    ) +
    theme_light() +
    theme(legend.position = c(0.8,0.9))
  
  #运行计时
  
  runningtime1<-c('generate random numbers additively',timeend1-timestart1,'s')
  runningtime2<-c('generate random numbers respectively',timeend2-timestart2,'s')
  
  #输出
  print(c(runningtime1,runningtime2))
  Plot

  
}

• 结果

Comp(m = 10000,n = 100,fun = 'Binomial',size = 1,prob = 0.3)

## [1] "generate random numbers additively"  
## [2] "0.00200009346008301"                 
## [3] "s"                                   
## [4] "generate random numbers respectively"
## [5] "0.0820000171661377"                  
## [6] "s"

Comp(m = 10000,n = 10,fun = 'Poisson',lambda = 2)

## [1] "generate random numbers additively"  
## [2] "0.000999927520751953"                
## [3] "s"                                   
## [4] "generate random numbers respectively"
## [5] "0.00600004196166992"                 
## [6] "s"

Comp(m = 10000,n = 10,fun = 'Normal',mean = 5, sd = 1.5)#正态分布标准化，n=1就可

## [1] "generate random numbers additively"  
## [2] "0.00199985504150391"                 
## [3] "s"                                   
## [4] "generate random numbers respectively"
## [5] "0.00999999046325684"                 
## [6] "s"

Comp(m = 10000,n = 100,fun = 'Gamma',shape = 2, rate = 3)

## [1] "generate random numbers additively"  
## [2] "0.000999927520751953"                
## [3] "s"                                   
## [4] "generate random numbers respectively"
## [5] "0.226000070571899"                   
## [6] "s"

Comp(m = 10000,n = 100,fun = 'Chi-Aquared',df = 3.5)

## [1] "generate random numbers additively"  
## [2] "0.00200009346008301"                 
## [3] "s"                                   
## [4] "generate random numbers respectively"
## [5] "0.226999998092651"                   
## [6] "s"