模擬simulation

?概率函數

概率函數通常用來生成特征已知的模擬數據，以及在統計函數中計算概率值。

對于任意分布有四種基本函數：

前綴作用

d	產生隨機數
r	估計概率分布的密度
p	估計累積分布
q	估計分位數

正態分布：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)

正態分布概率函數默認為標準正態分布。

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

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

set.seed(1)rnorm(5) #生成5個服從標準正態分布的隨機數[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

生成服從泊松分布的隨機數

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#泊松分布中，均值應該等于平均發生率

估計泊松分布的累積分布函數

ppois(2, 2)[1] 0.6766764#在平均發生率為2的泊松分布中出現小于等于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的正態分布。

當x服從正態分布時：

set.seed(20)x <- rnorm(100) #生成服從正態分布的隨機變量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和y之間存在線性趨勢。

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

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

set.seed(1)x <- rnorm(100) #生成服從正態分布的隨機變量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不是連續變量，而是計數變量，但x和y之間存在線性關系，當x越大，y也增加。?

?隨機抽樣 ?

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

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

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

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

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

分析器profiler

在開發規模較大的程序或者做大型數據分析時，運行代碼可能耗費很長時間，分析器可以分析代碼的運行時間并提供一些建議解決問題，優化代碼。

?system.time( )

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

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

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

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

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

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

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

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

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

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

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

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

?Rprof( )

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

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

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

library(ISLR)Wage Rprof(tmp lm(wage~age + jobclass + education,data=Wage) #運行lm函數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 #真正耗時的函數 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 #頂層函數為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分析器中得到的數據規范化的方式：

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

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

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

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

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。