如何用 R::animation 作動畫
December 22, 2010
上面這個圖是用下面的R
程式碼畫的
require(ggExtra) require(animation) x <- seq(-4.2, 4.2, length=10000) saveMovie( for (i in 1:10) { end.pin <- 1000*i dat1 <- data.frame(x=x[1:end.pin], y=pnorm(x[1:end.pin]),z=dnorm(x[1:end.pin])) h1 <- ggplot(dat1, aes(x=x, y=y)) h1 <- h1 + xlim(min(x), max(x)) + ylim(0, 1) h1 = dnorm(x[end.pin])) { h1 <- h1 + geom_ribbon(aes(y=y, ymin=0, ymax=y), fill='pink', alpha=0.75) h1 <- h1 + geom_ribbon(aes(y=z, ymin=0, ymax=z), fill='cyan', alpha=0.75) } else { h1 <- h1 + geom_ribbon(aes(y=z, ymin=0, ymax=z), fill='cyan', alpha=0.75) h1 <- h1 + geom_ribbon(aes(y=y, ymin=0, ymax=y), fill='pink', alpha=0.75) } h1 <- h1 + geom_text(aes(x=-2.4, y=0.7, label=paste('x=', round(x[end.pin],3))), size=8, colour='darkgreen') h1 <- h1 + geom_text(aes(x=0.2, y=0.2, label="φ(x)="), size=8, colour='blue') h1 <- h1 + geom_text(aes(x=2.0, y=0.2, label=round(dnorm(x[end.pin]), 3)), size=8, colour='blue') h1 <- h1 + geom_text(aes(x=-1.4, y=0.5, label="Φ(x)="), size=8, colour='red') h1 <- h1 + geom_text(aes(x=0.4, y=0.5, label=round(pnorm(x[end.pin]), 3)), size=8, colour='red') h1 <- h1 + geom_line(colour='red') h1 <- h1 + geom_line(aes(y=z), colour='blue') print(h1) }, moviename="normal_dist_ggplot2.gif")
R 的 integrate 問題
December 21, 2010
R::GSL
December 17, 2010
C++ SVD with Armabillo: 雖不滿意, 但可接受
December 17, 2010
用以下的程式碼:
#include <iostream> #include <armadillo> using namespace std; using namespace arma; int main(int argc, char* argv[]) { mat A = randn<mat>(1000,1000); mat U; vec S; mat V; bool status = svd(U,S,V,A); return 0; }
然後編譯執行:
$ g++ test2.cpp -o test2 -O2 -larmadillo $ time ./test2 13.284 secs
還可以接受啦(雖然還是很慢)
P.S.: 但如果要算很多次大矩陣, 恐怕就很難讓人接受了, 以下的程式碼:
#include <iostream> #include <armadillo> using namespace std; using namespace arma; int main(int argc, char* argv[]) { mat A = randn<mat>(1000,1000); mat U; vec S; mat V; int iter = 0; for (iter = 0; iter < 10; iter ++) { A = randn<mat>(1000,1000); bool status = svd(U,S,V,A); } return 0; }
跑了足足125.773秒, 而R
的快上很多:
> system.time(for(i in 1:10) y <- svd(matrix(rnorm(1e6), 1e3))) user system elapsed 35.360 0.510 19.596
這代表Armadillo
的延遲性是線性的, 這真的就不太好了
Calling C functions, R v.s. Python
December 16, 2010
Take a look on the following simple C
function:
double triplesum(double a, double b, double c) { return a + b + c; }
And we want to use it in Python
:
>>> from ctypes import * >>> lib1 = cdll.LoadLibrary('triplesum.so') >>> lib1.triplesum.argtypes = [c_double, c_double, c_double] >>> lib1.triplesum.restype = c_double >>> lib1.triplesum(1.2, 3.4, 2.3) 6.8999999999999995
In R
:
C
code:
void triplesum(double* a, double* b, double* c, double* d) { *d = *a + *b + *c; }
Compile it by
$ R CMD SHLIB forR.c
R
code:
> dyn.load("forR.so") > a <- 1.2; b <- 4.3; c <- 2.3; d .C("triplesum", a, b, c, d)[[4]] [1] 6.9
It’s hard to say which one is more convenient, but Python
seems better because you don’t need to modify the code.
R
‘s mechanism makes me feel like writing Fortran
subroutine but in C
syntax.
Google 到底想說什麼?
December 9, 2010
一門英烈
December 9, 2010
孔子和平獎
袓孫皆俊傑
磕頭相輝映
Bayesian posterior distribution 20101125-1
November 25, 2010
In Bayesian hierarchical setting:
\[
\left.\boldsymbol x_i\right|\boldsymbol\mu, \mathbf\Sigma\sim\mathscr N(\boldsymbol\mu, \mathbf\Sigma)
\]
$i=1,2,\ldots,n$, where
\[
\left.\boldsymbol\mu\right|\mathbf\Sigma\sim\mathscr N\left(\boldsymbol\mu_0, \frac{1}{\kappa_0}\mathbf\Sigma\right)
\]
and
\[
\mathbf\Sigma\sim\mathscr{IW}(\nu_0, \mathbf\Sigma_0)
\]
and, most statisticians can tell me that
\[
\left.\boldsymbol\mu\right|\boldsymbol x_1,\boldsymbol x_2,\ldots, \boldsymbol x_n,\mathbf\Sigma\sim\mathscr N\left(\frac{\kappa_0}{n+\kappa_0}\boldsymbol\mu_0+\frac{n}{n+\kappa_0}\overline{\boldsymbol x}, \frac{1}{\kappa_0}\mathbf\Sigma\right)
\]
where $\frac{1}{n}\sum_{i=1}^n\boldsymbol x_i$ and
\[
\left.\mathbf\Sigma\right|\boldsymbol x_1,\boldsymbol x_2,\ldots, \boldsymbol x_n\sim\mathscr{IW}\left(n+\nu_0, \mathbf\Sigma_0 + \sum_{i=1}^n(\boldsymbol x_i-\overline{\boldsymbol x})(\boldsymbol x_i-\overline{\boldsymbol x})^\top + \frac{\kappa_0n}{\kappa_0+n}(\boldsymbol x_i-\boldsymbol\mu_0)(\boldsymbol x_i-\boldsymbol\mu_0)^\top\right)
\]
right away, but, if we set
\[
\left.\boldsymbol\mu\right|\mathbf\Sigma\sim\mathscr N\left(\boldsymbol\mu_0,\mathbf{A\Sigma A}^\top\right)
\]
can we get such a beautiful posterior form?
BBC 新聞部表示
November 25, 2010
直排測試
November 24, 2010
朱雀橋邊野草花
烏衣巷口夕陽斜
舊時王謝堂前燕
飛入尋常百姓家
《般若波羅密多心經》
觀自在菩薩,行深般若波羅蜜多時,照見五蘊皆空,度一切苦厄。
舍利子,色不異空,空不異色;色即是空,空即是色。
受、想、行、識,亦復如是。
舍利子,是諸法空相,不生不滅,不垢不淨,不增不減,
是故空中無色,無受、想、行、識;
無眼、耳、鼻、舌、身、意;
無色、聲、香、味、觸、法;
無眼界,乃至無意識界;
無無明,亦無無明盡;
乃至無老死,亦無老死盡。
無苦、集、滅、道,無智亦無得,以無所得故。
菩提薩埵,依般若波羅蜜多故,心無罣礙,
無罣礙故,無有恐怖,遠離顛倒夢想,究竟涅槃。
三世諸佛,依般若波羅蜜多故,得阿耨多羅三藐三菩提。
故知般若波羅蜜多,是大神咒,是大明咒,是無上咒,是無等等咒,能除一切苦,真實不虛。
故說般若波羅蜜多咒,即說咒曰:揭諦揭諦,波羅揭諦,波羅僧揭諦,菩提薩婆訶。