模擬simulation

?概率函數(shù)

概率函數(shù)通常用來生成特征已知的模擬數(shù)據(jù)，以及在統(tǒng)計函數(shù)中計算概率值。

對于任意分布有四種基本函數(shù)：

前綴作用

d	產(chǎn)生隨機數(shù)
r	估計概率分布的密度
p	估計累積分布
q	估計分位數(shù)

正態(tài)分布：rnorm( )、dnorm( )、pnorm( )、qnorm( )； Gamma分布：rgamma( )、dgamma( )、pgamma( )和 qgamma( )； 泊松分布：rpois( )、dpois( )、ppois( ) 和 qpois( )。 dnorm(x, mean = 0, sd = 1, log = FALSE) pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) rnorm(n, mean = 0, sd = 1)

正態(tài)分布概率函數(shù)默認為標準正態(tài)分布。

?dnorm( ) 的log選項能求密度的對數(shù)值 (log of density) ?pnorm( ) 和 qnorm( ) 中l(wèi)og.p選項也可以對概率求對數(shù) ?lower.tail選項默認計算分布的左尾 (lower tail) ?若Φ是標準正態(tài)分布的累積分布函數(shù)，則pnorm(q)=Φ(q)且qnorm(p)=Φ^-1(p) (qnorm() 等于Φ的反函數(shù))

set.seed( ) 設置隨機數(shù)生成種子 (seed)，可以重復生成相同的隨機數(shù)。

set.seed(1)rnorm(5) #生成5個服從標準正態(tài)分布的隨機數(shù)[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078rnorm(5)[1] -0.8204684 0.4874291 0.7383247 0.5757814 -0.3053884set.seed(1)rnorm(5)[1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078

生成服從泊松分布的隨機數(shù)

x10, x[1] 1 1 1 0 2 1 2 0 1 1summary(x) Min. 1st Qu. Median Mean 3rd Qu. Max. 0 1 1 1 1 2#泊松分布中，均值應該等于平均發(fā)生率

估計泊松分布的累積分布函數(shù)

ppois(2, 2)[1] 0.6766764#在平均發(fā)生率為2的泊松分布中出現(xiàn)小于等于2的隨機變量的概率，Pr(x <= 2)ppois(4, 2)[1] 0.947347# Pr(x <= 4)ppois(6, 2)[1] 0.9954662# Pr(x <= 6)

?從線性模型模擬隨機值

假設線性模型y=β0+β1x+ε，其中β0等于0.5，β1等于2，隨機噪聲ε服從標準差為2的正態(tài)分布。

當x服從正態(tài)分布時：

set.seed(20)x <- rnorm(100) #生成服從正態(tài)分布的隨機變量xe <- rnorm(100, 0, 2) #生成隨機噪聲y<-0.5+2*x+e summary(y) Min. 1st Qu. Median Mean 3rd Qu. Max. -6.4084 -1.5402 0.6789 0.6893 2.9303 6.5052plot(x, y)

用回歸模型模擬x和y，兩者之間存在明顯線性關系，符合回歸模型。

當x服從二項分布時：

set.seed(10)x <- rbinom(100, 1, 0.5) e <- rnorm(100, 0, 2)y<-0.5+2*x+esummary(y) Min. 1st Qu. Median Mean 3rd Qu. Max. -3.4936 -0.1409 1.5767 1.4322 2.8397 6.9410 plot(x, y)

x是二元變量時，y依舊是連續(xù)的且服從正態(tài)分布，x和y之間存在線性趨勢。

?從廣義線性模型模擬隨機值

假設y服從泊松分布，Y ~ Poisson(μ)；log(μ)服從一個截距為β0，斜率為β1的線性函數(shù)，logμ=β0 +β1x；β0為0.5，β1為0.3。

set.seed(1)x <- rnorm(100) #生成服從正態(tài)分布的隨機變量xlog.mu<-0.5+0.3*xy <- rpois(100, exp(log.mu)) #模擬生成100個服從泊松分布的隨機變量summary(y) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 1.00 1.00 1.55 2.00 6.00plot(x, y)

可以看到y(tǒng)不是連續(xù)變量，而是計數(shù)變量，但x和y之間存在線性關系，當x越大，y也增加。?

?隨機抽樣 ?

sample函數(shù)可以從給定的一組對象中隨機抽取樣本。

set.seed(1)sample(1:10, 4)[1] 9 4 7 1 sample(1:10, 4)[1] 2 7 3 6 sample(letters, 5)[1] "r" "s" "a" "u" "w"

默認為無放回抽樣，不會得到重復的元素。

sample(1:10)[1] 10 6 9 2 1 5 8 4 3 7

如果不指定其他參數(shù)，返回這些整數(shù)重新排列后的結(jié)果，向量里1到10的順序被隨機重新排列。

sample(1:10, replace = TRUE)[1] 5 5 2 10 9 1 4 3 6 10

指定replace=TRUE，函數(shù)進行10次有放回抽樣。

分析器profiler

在開發(fā)規(guī)模較大的程序或者做大型數(shù)據(jù)分析時，運行代碼可能耗費很長時間，分析器可以分析代碼的運行時間并提供一些建議解決問題，優(yōu)化代碼。

?system.time( )

system.time( ) 取出并分析任意一條R語句，返回運行語句所需時間。如果執(zhí)行分析時代碼有錯誤，返回代碼運行直到錯誤發(fā)生所用的時間。?

●用戶時間 (user time)：運行這行語句需要占用中央處理器的時間。?

●運行時間 (elapsed time)：所經(jīng)歷的時間，也稱時鐘時間 (wall clock time)。 ?

?通常情況下，用戶時間和運行時間是比較接近的，計算機執(zhí)行函數(shù)或者語句所需時間大致上等于等待的時間。

?運行時間比用戶時間長，說明中央處理器實際可能耗費很多時間在其它后臺事項上。

?運行時間比用戶時間短，通常是因為機器有多核或多處理器且都處于可利用狀態(tài)。

#運行時間比用戶時間長system.time(readLines("http://www.google.cn")) #readLines() 從遠程服務器讀取網(wǎng)頁 user system elapsed 0.002 0.002 0.305

中央處理器實際上沒有耗費很多時間運行代碼，大部分時間用在等待網(wǎng)絡，等待數(shù)據(jù)經(jīng)過網(wǎng)絡傳遞到計算機里，等待網(wǎng)絡對數(shù)據(jù)的來回處理上。

#運行時間比用戶時間短hilbert function(n) { i<-1:n 1/outer(i-1,i,"+")}x 1000) system.time(svd(x)) #svd()計算這個矩陣的奇異值分解 用戶 系統(tǒng) 流逝 1.605 0.094 0.742

svd函數(shù)使用蘋果機的 Accelerate 框架，一個多線程線性代數(shù)庫，把計算任務分給兩個核進行。

#運行時間與用戶時間大致相等system.time({ n r for(i in 1:n){ x r[i] } }) 用戶 系統(tǒng) 流逝 0.099 0.006 0.109

不需要多線程或網(wǎng)絡活動的簡單表達式，用戶時間和運行時間基本相同。

?Rprof( )

Rprof( ) 啟動R中的分析器，可以在規(guī)律的樣本區(qū)間內(nèi)記錄追蹤函數(shù)調(diào)用棧，即每過0.02秒打印輸出函數(shù)的調(diào)用棧。(當函數(shù)運行時間少于0.02秒時無效)

system.time函數(shù)適用于較小型、較簡單、很明確問題發(fā)生位置的程序。如果不明確問題在哪里，可以用Rprof函數(shù)。(注意，不可以同時使用Rprof( )和sysyem.time( )。)

summaryRprof( )總結(jié)分析器的原始輸出結(jié)果，計算每個函數(shù)各自耗費的時間。

library(ISLR)Wage Rprof(tmp lm(wage~age + jobclass + education,data=Wage) #運行l(wèi)m函數(shù)Call:lm(formula = wage ~ age + jobclass + education, data = Wage)Coefficients: (Intercept) age 59.5800 0.5554 jobclass2. Information education2. HS Grad 4.5107 11.2009 education3. Some College education4. College Grad 23.3304 38.3862 education5. Advanced Degree 62.9115 Rprof()summaryRprof(tmp)$by.self #真正耗時的函數(shù) self.time self.pct total.time total.pct"model.matrix.default" 0.02 25 0.04 50"" 0.02 25 0.02 25"conditionCall" 0.02 25 0.02 25"stop" 0.02 25 0.02 25$by.total #頂層函數(shù)為lm，耗費0.06s total.time total.pct self.time self.pct"lm" 0.06 75 0.00 0"model.matrix.default" 0.04 50 0.02 25"model.matrix" 0.04 50 0.00 0"" 0.02 25 0.02 25"conditionCall" 0.02 25 0.02 25"stop" 0.02 25 0.02 25".External2" 0.02 25 0.00 0"[[.data.frame" 0.02 25 0.00 0"[[" 0.02 25 0.00 0"base::try" 0.02 25 0.00 0"contrasts 0.02 25 0.00 0"eval" 0.02 25 0.00 0"model.frame.default" 0.02 25 0.00 0……$sample.interval[1] 0.02$sampling.time[1] 0.08unlink(tmp)

sysyem.time( )將從R分析器中得到的數(shù)據(jù)規(guī)范化的方式：

?by total：把每個函數(shù)依照總運行時間排序，第一個通常是頂層函數(shù)

因為lm位于頂層，耗費75%的時間，但它只是調(diào)用真正進行運算的輔助函數(shù)，因此知道頂層函數(shù)使用的時間用處不大

?by self：先減去調(diào)用棧中低級函數(shù)花的時間再進行排序，即每個函數(shù)真正的耗費時間

編輯：李雪純 馮文清 校審：張健 羅鵬

總結(jié)

以上是生活随笔為你收集整理的R plot图片背景设置为透明_数据科学06 | R语言程序设计模拟和R分析器的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。