[問題] geometric distribution驗證weak large n

[AmigoSafin ()]

[問題類型]:

請把以下不需要的部份刪除
程式諮詢(我想用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://www.ptt.cc/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
→ andrew43: 你的apply誤用了。例如apply用來取最大值時，
→ andrew43: apply(m, 1, max) 就等於 apply(m,1,function(x)max(x))
→ andrew43: 所以你要適當把你apply中XX^2改成我上面function的寫法
→ AmigoSafin: 謝謝a大～～我有查過?apply但不太了解 好的！我來改改
→ AmigoSafin: a大真的是很厲害～
推 andrew43: 對了，檢查一下你手算的期望值。
→ andrew43: google "r apply anonymous function" 看些例子即可明白
哈囉大家好～

最後我還是沒有把這題做出來
想跟大家請益下我的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
→ andrew43: 先看頭幾行。seq(0,1)就是0和1而已。
→ andrew43: 另外每次rbern()你給的mu不是一個值，這樣可以嗎？
→ andrew43: 我是指你試試rbern(20, c(0,1))看看這是你要的嗎？
→ andrew43: 另外我建議先不要變動mu，先練習把mu=0.5寫成功就好
→ AmigoSafin: 謝謝Ａ大 我在練習看看～～感恩！！