[問題] Bernoulli distribution驗證weak large n
[問題類型]:
請把以下不需要的部份刪除
程式諮詢(我想用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
→
10/08 08:53,
6年前
, 1F
10/08 08:53, 1F
→
10/08 08:54,
6年前
, 2F
10/08 08:54, 2F
→
10/08 08:55,
6年前
, 3F
10/08 08:55, 3F
→
10/08 11:48,
6年前
, 4F
10/08 11:48, 4F
→
10/08 11:52,
6年前
, 5F
10/08 11:52, 5F
推
10/08 12:31,
6年前
, 6F
10/08 12:31, 6F
→
10/08 13:03,
6年前
, 7F
10/08 13:03, 7F
哈囉大家好～
最後我還是沒有把這題做出來
想跟大家請益下我的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
→
10/11 08:09,
6年前
, 8F
10/11 08:09, 8F
→
10/11 08:59,
6年前
, 9F
10/11 08:59, 9F
→
10/11 09:01,
6年前
, 10F
10/11 09:01, 10F
→
10/11 09:11,
6年前
, 11F
10/11 09:11, 11F
→
10/12 08:38,
6年前
, 12F
10/12 08:38, 12F
R_Language 近期熱門文章
PTT數位生活區 即時熱門文章
