kosugitti

par(mfrow=c(2,2))
x <- seq(-4,4,0.1)
plot(exp(x),main="eの累乗",xaxt="n",xlab="")
axis(side=1, at=1:length(x), labels=x) # ラベルの描画
plot(1/exp(x),main="ひっくり返す",xaxt="n",xlab="")
axis(side=1, at=1:length(x), labels=x) # ラベルの描画
plot(1/(1+exp(-x)),main="マイナスを掛ける",xaxt="n",xlab="")
axis(side=1, at=1:length(x), labels=x) # ラベルの描画
plot(1/(1+2.7^(-x)),main="近似値は2.7",xaxt="n",xlab="")
axis(side=1, at=1:length(x), labels=x) # ラベルの描画

kosugitti

# 正規分布のできるまで
plot(exp(-x),main="eの累乗",xaxt="n",xlab="")
axis(side=1, at=1:length(x), labels=x) # ラベルの描画
plot(exp(-(x^2)/2),main="二乗して二で割る",xaxt="n",xlab="")
axis(side=1, at=1:length(x), labels=x) # ラベルの描画
plot(1/sqrt(2*pi)*exp(-(x^2)/2),main="1/sqrt(2*pi)を掛ける",xaxt="n",xlab="")
axis(side=1, at=1:length(x), labels=x) # ラベルの描画
plot(1/sqrt(2*pi*3)*exp(-((x+1)^2)/(2*3)),main="平均-1，SD3にしてみた",xaxt="n",xlab="")
axis(side=1, at=1:length(x), labels=x) # ラベルの描画

kosugitti

# 参考
# http://d.hatena.ne.jp/isseing333/20100508/1273325218

library(ggplot2)

#---数式で記述
gauss.dencity <- function(x) 1/sqrt(2*pi)*exp(-x^2/2)
plot(gauss.dencity, -5, 5, xaxs="i",yaxs="i", bty="l", ylim=c(0, 0.5), las=1)

#---ダミーデータで描く

kosugitti

library(mvtnorm)
library(scatterplot3d)

############################################ 無相関
cor.mat <- matrix(c(1,0,0,1),nrow=2) 
x100 <- rmvnorm(n=100,mean=c(0,0),sigma=cor.mat)
x1000 <- rmvnorm(n=1000,mean=c(0,0),sigma=cor.mat)
x10000 <- rmvnorm(n=10000,mean=c(0,0),sigma=cor.mat)

kosugitti

library(MASS)
library(psych)
N <- 1000
cor.mat <- matrix(c(1,0.5,0.4,0.7,0.8,
                    0.5,1,0.3,0.4,0.4,
                    0.4,0.3,1,0.6,0.3,
                    0.7,0.4,0.6,1,0.9,
                    0.8,0.4,0.3,0.9,1),ncol=5)
# 乱数発生
dat <- as.data.frame(mvrnorm(N,mu=c(0,0,0,0,0),Sigma=cor.mat))

kosugitti

# ダミーデータによる理解
# 標準正規乱数を発生させて数え上げる
count.n.rand <- function(N,rg=seq(-4,4,0.1),mu=0,sd=1){
  X <- rnorm(N,mu,sd)
  lt <- list()
  for(i in 1:length(rg)){
    lt[i] <- length(X[X<rg[i]]) #データがある値を超えた数を数える
  }
  return(unlist(lt))
}

kosugitti

theta <- seq(-4,4,0.1)

# 1PL logistic model
# 関数
OnePL.func <- function(x,b) 1/(1+exp(-1.7*(x-b)))

# データ生成
dat <- data.frame(theta,b0=OnePL.func(theta,0),b1=OnePL.func(theta,1),bn1=OnePL.func(theta,-1))

# データ整形

kosugitti

# 2PL logistic model
# 関数
TwoPL.func <- function(x,a=1,b=0) 1/(1+exp(-1.7*a*(x-b)))

# データ生成...bは0で統一
dat <- data.frame(theta,
                  a1=TwoPL.func(theta,a=1,b=0),
                  a05=TwoPL.func(theta,a=0.5,b=0),
                  a2=TwoPL.func(theta,a=2,b=0))

kosugitti

####### bはmuに，aはsdに対応している
plot.dat <- data.frame(cbind(theta,
                             TwoPL.func(theta,a=1,b=0),   # 標準
                             TwoPL.func(theta,a=2,b=1),   # a=2,b=1
                             count.n.rand(N=50,sd=1,mu=0,rg=theta)/50,
                             count.n.rand(N=100,sd=1/2,mu=1,rg=theta)/100))
names(plot.dat) <- c("theta","default","a2b1","N50","N100")
# データの整形
library(tidyr)
plot.dat <- tidyr::gather(plot.dat,key=models,value=val,-theta)

kosugitti

# 3PL logistic model
# 関数
ThreePL.func <- function(x,a=1,b=0,c=0) c+(1-c)/(1+exp(-1.7*a*(x-b)))

# データ生成...bは0で統一
dat <- data.frame(theta,
                  ThreePL.func(theta,a=1,b=0,c=0.0),
                  ThreePL.func(theta,a=1,b=1,c=0.2),
                  ThreePL.func(theta,a=2,b=-1,c=0.4))