如何用 R::animation 作動畫

上面這個圖是用下面的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")

Advertisements

R 的 integrate 問題

測試一下:
\[ \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

用了一下子,只有一個感想, R::GSL完成度很高。

幾乎所有的函式都可以用gsl:: + TAB找到

看來R離完整的計算工具愈來愈近了,一個合用的計算工具應該具備以下的功能:

• 基本計算(至少不能比 C 標準函式庫的功能少), 同時要支援 IEEE的擴充規格 (可由 GSL 得到)
• 有應用數學常用函式, 比如說一本工數課本有提到的所有函數 (可由 GSL 得到)
• 統計常用函數 (這點 R 算是目前作得相當完整的)
• 能解偏微分方程, ODE 和 PDE (deSolve 差強人意, 算是堪用)
• 精確而快速的數值積分 (R 的 integrate單維積分算是堪用)
• 資料視覺化(就是畫圖啦, 我想不出開源軟體中有誰作得比R好的)

C++ SVD with Armabillo: 雖不滿意, 但可接受

用以下的程式碼:

#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

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 到底想說什麼?

一門英烈

鴉片有益論
孔子和平獎
袓孫皆俊傑
磕頭相輝映

Bayesian posterior distribution 20101125-1

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 新聞部表示

如此配圖, 這立場實在是不言可喻啊

直排測試

朱雀橋邊野草花
烏衣巷口夕陽斜
舊時王謝堂前燕
飛入尋常百姓家

《般若波羅密多心經》
觀自在菩薩，行深般若波羅蜜多時，照見五蘊皆空，度一切苦厄。
舍利子，色不異空，空不異色；色即是空，空即是色。
受、想、行、識，亦復如是。
舍利子，是諸法空相，不生不滅，不垢不淨，不增不減，
是故空中無色，無受、想、行、識；
無眼、耳、鼻、舌、身、意；
無色、聲、香、味、觸、法；
無眼界，乃至無意識界；
無無明，亦無無明盡；
乃至無老死，亦無老死盡。
無苦、集、滅、道，無智亦無得，以無所得故。
菩提薩埵，依般若波羅蜜多故，心無罣礙，
無罣礙故，無有恐怖，遠離顛倒夢想，究竟涅槃。
三世諸佛，依般若波羅蜜多故，得阿耨多羅三藐三菩提。
故知般若波羅蜜多，是大神咒，是大明咒，是無上咒，是無等等咒，能除一切苦，真實不虛。
故說般若波羅蜜多咒，即說咒曰：揭諦揭諦，波羅揭諦，波羅僧揭諦，菩提薩婆訶。