2016-12-06 26 views

回答

5

因此,這是你的函數現在(希望你知道怎麼寫的R功能;如果沒有,檢查writing your own function):

f <- function (x) (pi/2) * (1/(1 + 0.25 * x^2)) 

f(-Inf, Inf)定義,因此在這個範圍內的整合給出了一個不定積分。幸運的是,接近Infx^(-2)速度,所以整體是明確界定,並可以計算:

C <- integrate(f, -Inf, Inf) 
# 9.869604 with absolute error < 1e-09 

C <- C$value ## extract integral value 
# [1] 9.869604 

那麼你一定要規範化f,因爲我們知道一個概率密度應該集成到1:

f <- function (x) (pi/2) * (1/(1 + 0.25 * x^2))/C 

您可以通過繪製其密度:

curve(f, from = -10, to = 10) 

density

+0

像往常一樣,好的選擇。 – akrun

+0

謝謝,再次精彩的迴應!現在我有了可能的分佈函數,我想知道如何使用這個新的分佈函數創建一個說n = 1000的隨機樣本? – Chris95

+0

@ZheyuanLi完成。 – Chris95

1

既然我有可能的分佈函數,我想知道如何使用這個新的分佈函數創建一個隨機樣本n = 1000

一個離題的問題,但沒有你做一個新的線程就可以回答。它很有用,因爲它變得微妙。

enter image description here

比較

set.seed(0); range(simf(1000, 1e-2)) 
#[1] -56.37246 63.21080 
set.seed(0); range(simf(1000, 1e-3)) 
#[1] -275.3465 595.3771 
set.seed(0); range(simf(1000, 1e-4)) 
#[1] -450.0979 3758.2528 
set.seed(0); range(simf(1000, 1e-5)) 
#[1] -480.5991 8017.3802 

所以我覺得e = 1e-2是合理的。我們可以借鑑的樣本,做一個(縮放)直方圖和覆蓋密度曲線:

set.seed(0); x <- simf(1000) 
hist(x, prob = TRUE, breaks = 50, ylim = c(0, 0.16)) 
curve(f, add = TRUE, col = 2, lwd = 2, n = 201) 

enter image description here