[问题类型]: 请把以下不需要的部份删除 程式谘询(我想用R 做某件事情，但是我不知道要怎麽用R 写出来) [软体熟悉度]: 新手(没写过程式，R 是我的第一次) 入门(写过其他程式，只是对语法不熟悉) [问题叙述]: 请简略描述你所要做的事情，或是这个程式的目的 大家好 (以下数学符号以R markdown code表示 或许方便版友有兴趣查阅） 已知 1/n \sum_{i=1}^n Xi^2 ---> E[X^2] when n--> \infty Xi 是一个Bernoulli 分配 成功机率为theta 此处我设定p=theta=1/4 题目请我们运用R来对这个Bernoulli 分配做computational demonstration 来验证这个大数法则 大约是如此 如果我的解读的正确的 我的code如下 我把上课资料的Poisson试图改为rbinom 而且原来的code是看1/n \sum_{i=1}^n X_i 与E[X]之间的距离 现在我要求的是 1/n \sum X_i^2 与 E[X_i^2]之间的距离 所以我修改了code 但显然 我修改得有误 我们已被提醒要先求出E[X^2] 我求出 E[X^2] = V[X]+E[X]^2=(p-2)(p-1)/p^2 就把上述这值与code中的mu做调换 [程式范例]: m <- 100 epsilon <- 0.01 #vn <- c(100, 500, 1600, 3200, 6400, 9600, 18000, 25600, 54000, 108000, 256000, 819200, 1000000) vn <- c(100, 500, 1600, 3200, 6400, 9600, 18000, 25600, 54000, 108000, 256000, 819200) nn <- length(vn) p.good <- numeric(nn) p <- 1/4 mu <- (p-2)*(p-1)/p^2 #mu <- lambda <- 3 ?rbinom for(j in 1:nn) { n <- vn[j] XX <- matrix(rbinom(n*m, size=1, p),ncol=n) #XX <- matrix(rnorm(n*m,mean=9,sd=2),ncol=n) Xbar <- apply(XX,1,XX^2) good <- which(abs(Xbar-mu)<epsilon) p.good[j] <- length(good)/m } windows() plot(vn, p.good, type='b', xlab='n', ylab='Prob[|Xbar-mu|<epsilon]') ########################################################################### 以下是原来的程式码 我是初学者 只能先从改别人的程式码开始. # Computational demonstration of the LLN (Law of Large Numbers # This first demonstration uses the normal distribution # However the LLN applies to all possible distrbutions # m <- 100 epsilon <- 0.01 #vn <- c(100, 500, 1600, 3200, 6400, 9600, 18000, 25600, 54000, 108000, 256000, 819200, 1000000) vn <- c(100, 500, 1600, 3200, 6400, 9600, 18000, 25600, 54000, 108000, 256000, 819200) nn <- length(vn) p.good <- numeric(nn) mu <- 9 mu <- lambda <- 3 for(j in 1:nn) { n <- vn[j] XX <- matrix(rpois(n*m,mu),ncol=n) #XX <- matrix(rnorm(n*m,mean=9,sd=2),ncol=n) Xbar <- apply(XX,1,mean) good <- which(abs(Xbar-mu)<epsilon) p.good[j] <- length(good)/m } windows() plot(vn, p.good, type='b', xlab='n', ylab='Prob[|Xbar-mu|<epsilon]') windows() par(mfrow=c(2,2)) [环境叙述]: 请提供 sessionInfo() 的输出结果， 里面含有所有你使用的作业系统、R 的版本和套件版本资讯， 让版友更容易找出错误 [关键字]: Geometric distribution 选择性，也许未来有用 -- --

※ 发信站: 批踢踢实业坊(ptt.cc), 来自: 8.41.66.201 ※ 文章网址: https://webptt.com/cn.aspx?n=bbs/R_Language/M.1538930158.A.B7B.html ※ 编辑: AmigoSafin (8.41.66.201), 10/08/2018 00:44:14 ※ 编辑: AmigoSafin (8.41.66.201), 10/08/2018 05:22:49 哈罗大家好～ 最後我还是没有把这题做出来 想跟大家请益下我的code哪里需要再改进 首先 题目是： (1/n)*sum\ X_i^2 ---> E[X^2] 当 n--> 无限大 请： 对X_i ~ Bernoulli(theta), theta 范围(0,1) 做出empirical domonstration 提醒请我们须先清楚载明E[X^2]是多少 以免题目变得复杂 就我所知 Bernoulli的E[X^2]还是p 所以我就设定mu值是theta 并且在(0,1)之间 我的code如下 结果并没有跑出我想要的plot 还请大家不吝指导 让我可以知道自己错在哪里 # Bernoulli #par(mfrow=c(2,2)) library(Rlab) m <- 100 epsilon <- 0.01 #vn <- c(100, 500, 1600, 3200, 6400, 9600, 18000, 25600, 54000, 108000, 256000, 819200, 1000000) vn <- c(100, 500, 1600, 3200, 6400, 9600, 18000, 25600, 54000, 108000, 256000, 819200) nn <- length(vn) p.good <- numeric(nn) mu <- seq(from=0, to=1) #mu <- lambda <- 3 for(j in 1:nn) { n <- vn[j] f <-function(X){(1/n*sum(X^2))-mu} #p <-seq(0,1) XX <-matrix(rbern(n*m,mu), ncol = n) #XX <- matrix(rpois(n*m,mu),ncol=n) #XX <- matrix(rnorm(n*m,mean=9,sd=2),ncol=n) Xbar <- apply(XX,1,f) # already made a function good <- which(abs(Xbar-mu)<epsilon) p.good[j] <- length(good)/m} #windows() par(mfrow=c(2,2)) plot(vn, p.good, type='b', xlab='n', ylab='Prob[|Xbar-mu|<epsilon]') 这题接续还有normal poisson & 好多分配要改 要学好R真是不容易 ※ 编辑: AmigoSafin (129.21.70.153), 10/11/2018 00:08:56