おっぱい関数

「おっぱい関数」を材料として3次元表示についてまとめる。

「おっぱい関数」は元々誰が作ったか不詳だが，検索すると
Julia + Excel　http://www.hirax.net/diaryweb/2012/03/01.html
Maxima　https://plus.google.com/u/0/117731699005945352808/posts/Yswj5tswtNd
などが引っかかる。変形したものもある。
Python　http://kemeconoajito.blog88.fc2.com/blog-entry-175.html
R によるものは，http://siguniang.wordpress.com/2011/11/09/3d-surface-plot-with-latticergl/
にある。

library(lattice)
 
g <- function(x, y) {
 1/8 * (6 * exp(-((2/3 * abs(x) - 1)^2 + (2/3 * y)^2) - 1/3 * (2/3 * y + 1/2)^3) +
        2/3 * exp(-2.818^11 * ((abs(2/3 * x) - 1)^2 + (2/3 * y)^2)^2) +
        2/3 * y - (2/3 * x)^4)
}
m <- 50
x <- seq(-3, 3, length.out=m)
y <- seq(-3, 3, length.out=m)
grid   <- expand.grid(x=x, y=y)
grid$z <- g(grid$x, grid$y)
wireframe(z ~ x * y, grid, shade=TRUE)
library(rgl)
grid$col <- ifelse(sqrt((abs(grid$x) - 3/2) ^ 2 + grid$y ^2 ) < 0.45, 'hotpink', 'lightpink')# nipple
open3d()
rgl.surface(x, y, grid$z, color=grid$col)

この関数部分を小生が書き直した

n <- 100
x0 <- y0 <- seq(-3, 3, length=n+1)
xy <- expand.grid(x=x0, y=y0)
x <-  xy[,1]
y <-  xy[,2]
t0 <- (2 * abs(x) / 3 - 1) ^ 2
t1 <- 2 * y / 3
t2 <- t1 ^ 2
t4 <- -t0 - t2 - (t1 + 0.5) ^ 3 / 3
z <- 0.125 * (6 * exp(t4) +
     2 * exp(-exp(11) * (t0 + t2) ^ 2) / 3 +
     2.5 * t1 -16 * x ^ 4 / 81)
library(rgl)
plot3d(x, y, z, col = "burlywood1", aspect = c(1, 1, 0.6), axes=FALSE, xlab="", ylab="", zlab="")

# 以下の，ふたつの方法もやってみる
open3d()
rgl.surface(x0, y0, z, color="burlywood1")
library(lattice)
wireframe(z ~ x * y, xy, shade=TRUE, aspect=c(1, 0.5), xlab="", ylab="", zlab="")

# 以下で出力される csv ファイルを Excel で 3d 表示
write(matrix(z, n+1), file="oppai2.csv", ncolumns=n+1, sep=",")

変形したものも書いておく。
http://kemeconoajito.blog88.fc2.com/blog-entry-175.html

n <- 100
xy <- expand.grid(seq(-3, 3, length=n+1), seq(-3, 3, length=n+1))
x <-  xy[,1]
y <-  xy[,2]
t0 <- (2 * abs(x) / 3 - 1) ^ 2
t1 <- 2 * y / 3
t2 <- t1 ^ 2
t3 <- (t0 + t2) ^ 2
t4 <- -t0 - t2 -  (t1 + 0.5) ^ 3 / 3
z <- 0.125 * (7 * exp(t4) + 1.5 * exp(-exp(10) * t3) +
     exp(-exp(5) * t3) + 2.5 * t1 - 16 * x ^ 4 / 81)
library(rgl)
plot3d(x, y, z, col = "burlywood1", aspect = c(1, 1, 0.6), axes=FALSE, xlab="", ylab="", zlab="")

ついでに見つかった，2次元の「おっぱい関数」も R で書いたものを載せておく。

# https://twitter.com/math_neko/status/227672102487068672/photo/1
y <- seq(.Machine$double.eps, 1.2, length=300)
x <- 3*y*log(y)-1/36*exp(-(36*y-36/exp(1))^4)
plot(1, 1, type="n", xlim=c(-1.2, 1), ylim=c(-0.1, 1))
lines(x, y)

# http://goo.gl/rtZS7
n <- 5000
curve(sqrt(1-(x-1)^2), n=n, xlim=c(-2.5, 2.5), ylim=c(-0.1,  1.25), asp=1)
curve(sqrt(1-(x+1)^2), n=n, add=TRUE)
curve(sqrt(0.01-(x-1)^2)+1, n=n, add=TRUE)
curve(sqrt(0.01-(x+1)^2)+1, n=n, add=TRUE)

# http://i.imgur.com/IHAhR.jpg
t <- seq(-1.4, 1.5, length=200)
y <- 2*t^3-0.5*t-1
x <- 0.9*(exp(-t^4+t^2+t)+0.2*(exp(-100*((y-0.2)^2)))-0.5*exp(-(y-3.5)^2)-0.15*exp(-(y-5)^2))
plot(x, y, type="l", asp=1)

ダメ出し：ちゃんと関数があるんだから，モ～

久しぶりの大物

「多次元度数分布を取る」より

プログラム自体は，やっと <- の前後に空白をおくなどで少しは読みやすくなったけど，まだ，){ や if( や for( があったり，カンマの後に空白を置かないなど相変わらず，読みにくい。

3次元クロス集計をする自前の関数を書いたというところだけど，リストや for 文を惜しげもなく使っているので，処理時間がべらぼうにかかる。

> set.seed(123456789)
> X <- matrix(rnorm(100000*3),ncol=3) + matrix(rnorm(100000*3,mean=3,sd=2),ncol=3)
>
> ns <- c(5,4,3)
> system.time(mult.hist.out <- multi.hist(X,ns = ns))
   ユーザ   システム       経過  
    21.665      0.346     21.797

これでは困る。R にはちゃんとした関数 xtabs があるのだから，それ使おうヨ。

ただ，xtabs は，元のデータがカテゴリー化されているというのが前提なので，カテゴリー化は cut 関数を使って自分で行う。

もとのプログラムも cut 関数で使う breaks と同じ働きをする変数の設定はしているが，ちょっと小細工をしている。しかし，「*0.0001」はいいかげんなものだ。意味不明のオマジナイノジュモンは使わない。

    if(is.null(bins)){
        range.X <- apply(X,2,range)
        range.X[2,] <- range.X[2,] + min(range.X[2,]-range.X[1,])*0.0001
        bins <- list()
        for(i in 1:length(X[1,])){
            bins[[i]] <- seq(from=range.X[1,i],to=range.X[2,i],length=ns[i]+1)
        }
    }

正統的なプログラムにしよう。
div 関数が元のデータのカテゴリー化を行う。cut 関数で使う breaks は，一番最初の区間の左端は含まれないので（マニュアル嫁） seq(range.x[1], range.x[2], length=ns+1) だけだと最小値の行き場所がなくなる（NA になる）ので，.Machine$double.eps を引いてやる。

mx <- function(x) {
    x <- data.frame(x)
    div <- function(x, ns) {
        range.x <- range(x)
        breaks <- seq(range.x[1], range.x[2], length=ns+1)
        breaks[1] <- breaks[1]-.Machine$double.eps
        return(as.integer(cut(x, breaks=breaks)))
    }    
    return(xtabs(~., mapply(div, x, ns)))
}

このプログラムの実行時間は

> system.time(mx(X))
   ユーザ   システム       経過  
     0.810      0.034      0.842

となり，26倍ほど速い。

自分でプログラムするにしても，もう少し速いプログラムを書ける（ようになってネ）。

ちなみに，元記事のプログラム名は multi.hist となっているが，hist はヒストグラム（度数分布図）。
度数分布表は frequency table，クロス集計表は contingency table とかいうから，名前を変えたほうがよい。