廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

r2stlパッケージの紹介:R から3Dプリンターに出力

r2stl

ぼくもよく理解してないんだが、ともあれ、stl形式にすれば MakerBot Desktop から3Dプリンターに出力できるみたいだ。

二変量正規分布の密度関数つくって「相関係数いくつくらいだと思う?」っていうクイズ出したり、尤度関数の模型つくって「最尤推定っていうのは要はこの点を探してるんですよ〜」とかやったらおもしろいかもしれない(が、それをおもしろいと感じるやつが日本に何人いるんだろう)

ともあれ、まずは r2stlパッケージをインストール
CRAN - Package r2stl

例1 volcano

library(r2stl)
#dim(volcano) #87行 61列
x <- 1:nrow(volcano)
y <- 1:ncol(volcano)
r2stl(x, y, volcano, filename="volcano.stl", show.persp=TRUE)

例2 二変量正規分布

二変量正規分布の密度関数は、平均=0、分散=1の場合でも
f(x,y)=\frac{1}{2 \pi \sqrt{1- r^2} } \exp\left(-\frac{1}{2(1- r^2)}(x^2 -2  r   y + y^2) \right)
という長い式なので、一回くらい手で書いてみたらなんかの思いでになるでしょう。
r は相関係数)

#library(r2stl)
dnorm2 <- function(x, y, r=0.7){
  det <- 1- r^2
  return(
    1/(2 * pi * sqrt(det)) * exp(-(x^2 -2 * r * y + y^2)/(2*det) )
  )
}
x <- seq(-3,3,length=100)
#seqは-3から3まで長さ100のベクトル(数列)を作ってるだけなので
#図形をもっと滑らかにしたいなら
#length=10000とか適当な数字を与えてやればいい
y <- x
z <- outer(x, y, dnorm2)
#
r2stl(x, y, z, filename="dnorm2.stl", show.persp=TRUE)

例3 尤度関数

正規分布のパラメーターの対数尤度関数は、途中の計算式は省略して結果だけ書くと、

\log L(\mu, \sigma^2 ; \bf{x} )=\frac{1}{2} \left(-n \log(2 \pi)- n \log(\sigma^2)-\frac{\sum{(x_i-\mu)^2}}{\sigma^2}\right)

なので、以下のようになる。

library(r2stl)
ll = function(x){
  n =length(x)
 return(function(mu, sigma2){
  -(-n*log(2*pi)- n*log(sigma2)-(sum((x-mu)^2))/sigma2)/2
  #対数尤度関数の値は負になるのでマイナスをかけてる
  #なのでこの関数を「最小」にする点が最尤推定値になる
  })
}
x1=rnorm(100,mean=3,sd=sqrt(9))
opt = ll(x1)
xs <- seq(1, 5, length= 100)
ys <- xs
z <- outer(xs,ys, opt)
#
r2stl(xs, ys, z, filename="llnorm.stl", show.persp=TRUE)

デモ


参考文献

入門・演習 数理統計

入門・演習 数理統計

R/S‐PLUSによる統計解析入門

R/S‐PLUSによる統計解析入門