99999久久久久久亚洲,欧美人与禽猛交狂配,高清日韩av在线影院,一个人在线高清免费观看,啦啦啦在线视频免费观看www

熱線電話:13121318867

登錄
首頁(yè)精彩閱讀標(biāo)準(zhǔn)正態(tài)分布函數(shù)的快速計(jì)算方法
標(biāo)準(zhǔn)正態(tài)分布函數(shù)的快速計(jì)算方法
2016-01-31
收藏

標(biāo)準(zhǔn)正態(tài)分布函數(shù)的快速計(jì)算方法

標(biāo)準(zhǔn)正態(tài)分布的分布函數(shù) $\Phi(x)$ 可以說(shuō)是"數(shù)據(jù)分析師"統(tǒng)計(jì)計(jì)算中非常重要的一個(gè)函數(shù),基本上有正態(tài)分布的地方都或多或少會(huì)用上它。在一些特定的問(wèn)題中,我們"數(shù)據(jù)分析師"需要大量多次地計(jì)算這個(gè)函數(shù)的取值,比如我經(jīng)常需要算正態(tài)分布與另一個(gè)隨機(jī)變量之和的分布,這時(shí)候就需要用到數(shù)值積分,而被積函數(shù)就包含 $\Phi(x)$。如果 $Z\sim N(0,1), X\sim f(x)$,$f$ 是 $X$ 的密度函數(shù),那么 $Z+X$ 的分布函數(shù)就是

我們"數(shù)據(jù)分析師"知道,$\Phi(x)$ 沒(méi)有簡(jiǎn)單的顯式表達(dá)式,所以它需要用一定的數(shù)值方法進(jìn)行計(jì)算。在大部分的科學(xué)計(jì)算軟件中,計(jì)算的精度往往是第一位的,因此其算法一般會(huì)比較復(fù)雜。當(dāng)這個(gè)函數(shù)需要被計(jì)算成千上萬(wàn)次的時(shí)候,速度可能就成為了一個(gè)瓶頸。

當(dāng)然有問(wèn)題就會(huì)有對(duì)策,一種常見(jiàn)的做法是略微放棄一些精度,以換取更簡(jiǎn)單的計(jì)算。在大部分實(shí)際應(yīng)用中,一個(gè)合理的誤差大小,例如 $10^{-7}$,一般就足夠了。在這篇文章中,給大家介紹兩種簡(jiǎn)單的方法,它們都比R中自帶的 pnorm() 更快,且誤差都控制在 $10^{-7}$ 的級(jí)別。

第一種辦法來(lái)自于經(jīng)典參考書(shū) Abramowitz and Stegun: Handbook of Mathematical Functions 的 公式 26.2.17 。其基本思想是把 $\Phi(x)$ 表達(dá)成正態(tài)密度函數(shù) $\phi(x)$ 和一個(gè)有理函數(shù)的乘積。這種辦法可以保證誤差小于 $7.5\times 10^{-8}$,一段C++實(shí)現(xiàn)可以在 這里 找到。(代碼中的常數(shù)與書(shū)中的略有區(qū)別,是因?yàn)榇a是針對(duì)誤差函數(shù) $\mathrm{erf}(x)$ 編寫(xiě)的,它與 $\Phi(x)$ 相差一些常數(shù))

我們來(lái)對(duì)比一下這種方法與R中 pnorm() 的速度,并驗(yàn)證其精度。

library(Rcpp)
sourceCpp("test_as26217.cpp") x = seq(-6, 6, by = 1e-6) system.time(y <- pnorm(x)) ## user  system elapsed ## 1.049   0.000   1.048 system.time(asy <- r_as26217ncdf(x)) ## user  system elapsed ## 0.293   0.019   0.311 max(abs(y - asy)) ## [1] 6.968772e-08

可以看出,A&S 26.2.17 的速度大約是 pnorm() 的三倍,且誤差也在預(yù)定的范圍里,是對(duì)計(jì)算效率的一次巨大提升。

那么還有沒(méi)有可能更快呢?答案是肯定的,而且你其實(shí)已經(jīng)多次使用過(guò)這種方法了。怎么,不相信?看看下面這張圖,你就明白了。

沒(méi)錯(cuò),這種更快的方法其實(shí)就是兩個(gè)字:查表。它的基本想法是,我們預(yù)先計(jì)算好一系列的函數(shù)取值 $(x_i,\Phi(x_i))$,然后當(dāng)我們需要計(jì)算某個(gè)點(diǎn) $x_0$ 時(shí),就找到離它最近的兩個(gè)點(diǎn) $x_k$ 和 $x_{k+1}$,再用線性插值的方法得到 $\Phi(x_0)$ 的近似取值:

什么?覺(jué)得這個(gè)方法太簡(jiǎn)單了?先別急,這里面還有不少學(xué)問(wèn)。之前我們"數(shù)據(jù)分析師"說(shuō)了,我們需要保證這種方法的誤差不超過(guò) $\epsilon=10^{-7}$,因此就需要合理地選擇預(yù)先計(jì)算的點(diǎn)。由于 $\Phi(-x)=1-\Phi(x)$,我們暫且只需要考慮 $x$ 為正的情況。如果讓 $x_i = ih,i=0,1,\ldots,N$,那么對(duì)函數(shù) $f$ 進(jìn)行線性插值的誤差將不超過(guò)( 來(lái)源 )

其中 $\Vert f’’ \Vert_{\infty}$ 是函數(shù)二階導(dǎo)絕對(duì)值的最大值。對(duì)于正態(tài)分布函數(shù)來(lái)說(shuō),它等于 $\phi(1)\approx 0.242$。于是令 $E(x)=10^{-7}$,我們就可以解出 $h\approx 0.001818$。最后,只要 $x_N>5.199$,即 $N\ge 2860$ 并另所有 $x>x_N$ 的取值等于1,就可以保證整個(gè)實(shí)數(shù)域上 $\Phi(x)$ 的近似誤差都不超過(guò) $10^{-7}$。

這種簡(jiǎn)單方法的實(shí)現(xiàn)我放在了 Github 上 ,源程序和測(cè)試代碼也可以在文章最后找到。下面給出它的表現(xiàn):

library(Rcpp)
sourceCpp("test_fastncdf.cpp") x = seq(-6, 6, by = 1e-6) system.time(fasty <- r_fastncdf(x)) ## user  system elapsed ## 0.043   0.024   0.066 max(abs(y - fasty)) ## [1] 9.99999e-08

與之前的結(jié)果相比,相當(dāng)于速度是 pnorm() 的15倍!

我們似乎一直以為,在計(jì)算機(jī)和統(tǒng)計(jì)軟件普及以后,一些傳統(tǒng)的做法就會(huì)慢慢被淘汰,例如現(xiàn)在除了考試,或許大部分的時(shí)間我們都是在用軟件而不是正態(tài)概率表。從教學(xué)與實(shí)際應(yīng)用的角度來(lái)看,這種做法是 應(yīng)該進(jìn)行推廣和普及的 ,但這也不妨礙我們從一些“舊知識(shí)”中汲取營(yíng)養(yǎng)。關(guān)于這種大巧若拙的做法的故事還有很多,比如廣為流傳的 這一則 。在計(jì)算資源匱乏的年代,數(shù)據(jù)科學(xué)家"數(shù)據(jù)分析師"們想出了各種巧妙的辦法來(lái)解決他們遇到的各種問(wèn)題?,F(xiàn)如今計(jì)算機(jī)的性能已經(jīng)遠(yuǎn)不是當(dāng)年可以媲跡,但前人的很多智慧卻依然穿透了時(shí)間來(lái)為現(xiàn)在的我們提供幫助,不得不說(shuō)這也是一種緣分吧。

數(shù)據(jù)分析咨詢請(qǐng)掃描二維碼

若不方便掃碼,搜微信號(hào):CDAshujufenxi

數(shù)據(jù)分析師資訊
更多

OK
客服在線
立即咨詢
客服在線
立即咨詢
') } function initGt() { var handler = function (captchaObj) { captchaObj.appendTo('#captcha'); captchaObj.onReady(function () { $("#wait").hide(); }).onSuccess(function(){ $('.getcheckcode').removeClass('dis'); $('.getcheckcode').trigger('click'); }); window.captchaObj = captchaObj; }; $('#captcha').show(); $.ajax({ url: "/login/gtstart?t=" + (new Date()).getTime(), // 加隨機(jī)數(shù)防止緩存 type: "get", dataType: "json", success: function (data) { $('#text').hide(); $('#wait').show(); // 調(diào)用 initGeetest 進(jìn)行初始化 // 參數(shù)1:配置參數(shù) // 參數(shù)2:回調(diào),回調(diào)的第一個(gè)參數(shù)驗(yàn)證碼對(duì)象,之后可以使用它調(diào)用相應(yīng)的接口 initGeetest({ // 以下 4 個(gè)配置參數(shù)為必須,不能缺少 gt: data.gt, challenge: data.challenge, offline: !data.success, // 表示用戶后臺(tái)檢測(cè)極驗(yàn)服務(wù)器是否宕機(jī) new_captcha: data.new_captcha, // 用于宕機(jī)時(shí)表示是新驗(yàn)證碼的宕機(jī) product: "float", // 產(chǎn)品形式,包括:float,popup width: "280px", https: true // 更多配置參數(shù)說(shuō)明請(qǐng)參見(jiàn):http://docs.geetest.com/install/client/web-front/ }, handler); } }); } function codeCutdown() { if(_wait == 0){ //倒計(jì)時(shí)完成 $(".getcheckcode").removeClass('dis').html("重新獲取"); }else{ $(".getcheckcode").addClass('dis').html("重新獲取("+_wait+"s)"); _wait--; setTimeout(function () { codeCutdown(); },1000); } } function inputValidate(ele,telInput) { var oInput = ele; var inputVal = oInput.val(); var oType = ele.attr('data-type'); var oEtag = $('#etag').val(); var oErr = oInput.closest('.form_box').next('.err_txt'); var empTxt = '請(qǐng)輸入'+oInput.attr('placeholder')+'!'; var errTxt = '請(qǐng)輸入正確的'+oInput.attr('placeholder')+'!'; var pattern; if(inputVal==""){ if(!telInput){ errFun(oErr,empTxt); } return false; }else { switch (oType){ case 'login_mobile': pattern = /^1[3456789]\d{9}$/; if(inputVal.length==11) { $.ajax({ url: '/login/checkmobile', type: "post", dataType: "json", data: { mobile: inputVal, etag: oEtag, page_ur: window.location.href, page_referer: document.referrer }, success: function (data) { } }); } break; case 'login_yzm': pattern = /^\d{6}$/; break; } if(oType=='login_mobile'){ } if(!!validateFun(pattern,inputVal)){ errFun(oErr,'') if(telInput){ $('.getcheckcode').removeClass('dis'); } }else { if(!telInput) { errFun(oErr, errTxt); }else { $('.getcheckcode').addClass('dis'); } return false; } } return true; } function errFun(obj,msg) { obj.html(msg); if(msg==''){ $('.login_submit').removeClass('dis'); }else { $('.login_submit').addClass('dis'); } } function validateFun(pat,val) { return pat.test(val); }