如何用 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
測試一下:
\[ \Phi(\color{red}x) = \int_{-\infty}^{\color{red}x} \phi(z)dz\]
> integrate(dnorm,0,4384) 0.5 with absolute error < 5.2e-05 > integrate(dnorm,0,4385) 6.124136e-05 with absolute error < 0.00012
理論上應該是要遞增才對﹐ 一下掉到快變零,這是為什麼?
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
朱雀橋邊野草花
烏衣巷口夕陽斜
舊時王謝堂前燕
飛入尋常百姓家
《般若波羅密多心經》
觀自在菩薩，行深般若波羅蜜多時，照見五蘊皆空，度一切苦厄。
舍利子，色不異空，空不異色；色即是空，空即是色。
受、想、行、識，亦復如是。
舍利子，是諸法空相，不生不滅，不垢不淨，不增不減，
是故空中無色，無受、想、行、識；
無眼、耳、鼻、舌、身、意；
無色、聲、香、味、觸、法；
無眼界，乃至無意識界；
無無明，亦無無明盡；
乃至無老死，亦無老死盡。
無苦、集、滅、道，無智亦無得，以無所得故。
菩提薩埵，依般若波羅蜜多故，心無罣礙，
無罣礙故，無有恐怖，遠離顛倒夢想，究竟涅槃。
三世諸佛，依般若波羅蜜多故，得阿耨多羅三藐三菩提。
故知般若波羅蜜多，是大神咒，是大明咒，是無上咒，是無等等咒，能除一切苦，真實不虛。
故說般若波羅蜜多咒，即說咒曰：揭諦揭諦，波羅揭諦，波羅僧揭諦，菩提薩婆訶。
Gong-Yi's SandBox@WordPress 