[問題類型]: 請把以下不需要的部份刪除 程式諮詢(我想用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/m.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