R語言由Monte Carlo方法計算積分

2021-08-13 16:07:49 字數 1332 閱讀 9614

蒙特·卡羅方法(monte carlo method),也稱統計模擬方法,是二十世紀四十年代中期由於科學技術的發展和電子計算機的發明,而被提出的一種以概率統計理論為指導的一類非常重要的數值計算方法。是指使用隨機數(或更常見的偽隨機數)來解決很多計算問題的方法。

f = function(o) sum(rexp(o, rate = 1/1000))

g = vectorize(f)

p = function(n))

c(mean = mean(b), sd = sd(b))

}

# 簡單 monte carlo 方法

p1 = function(n))

c(mean = mean(b), sd = sd(b))}

# importance sampling 方法

p2 = function(n))

c(mean = mean(smp), sd = sd(smp))}

f <- function(x) exp(- x^(1/2))*(sin(x))^2

g <- function(x) (1/(2*pi))*(1/(1 + x^2/4))*2

ratio <- function(x) f(x)/g(x)

# 首先看看如何產生服從柯西分布的正半軸的密度函式的隨機數

h <- function(n) abs(rcauchy(n, scale = 2))

a = h(100000)

hist(a[a<=40], freq = false, breaks = seq(0, 40, 0.05), border = "darkblue")

curve(g, from = 0, to = 40, add = true, col = "red")

# 下面用 importance sampling 方法生成服從 f(x) 的隨機數

)c(mean = mean(d), sd = sd(d))}

smp(100, 100)

smp(100, 1000)

smp(100, 100000)

# 檢查上面結果是否正確

w = function(x) x*f(x)

integrate(w, lower = 0, upper = 1000)

# 注意到均值的標準差是和根號 n 成反比。

a = null

for(i in seq(1000, 10000, 1000)) a = c(a, smp(100, i)["sd"])

# 得 n = 7000

中子穿牆問題的MonteCarlo求解方法

如下圖所示,代表乙個中子穿過用於遮蔽中子的鉛牆的示意圖。鉛牆的高度遠大於左右的厚度。設中子垂直由左端進入鉛牆,鉛牆中前行乙個單位距離後與乙個鉛原子碰撞。此時會改變方向,但是此時改變的方向是任意的,並且能繼續執行乙個單位後與另乙個鉛原子碰撞。這樣下去,如果中子在鉛牆裡消耗掉所有的能量或者從左端溢位即被...

R語言學習 卡方檢驗

本文是個人學習筆記 卡方檢驗用來檢驗類別變數。gender f table data gender,data group print gender f 12 31616 9277 10 行是性別,列是組 prop.table gender f,2 注意,這裡2代表 gender f 中的第二個變數,...

r語言electricity資料集 R語言 資料集

第二章 建立資料集 1.r語言的資料型別 數值型 字元型 邏輯型 複數型 虛數 和原生型 位元組 2.資料結構 a.向量 儲存數值型 字元型和邏輯型資料的一維陣列 a c 1,2,3,4,5 建立 組合功能的函式c a 1 1 2 3 4 5 a c 1,4 訪問 方括號 向量中指定的元素 1 1 ...