カテゴリー別アーカイブ: 統計学・プログラミング

ロジスティック写像

本格的にJAVAを学ぼうと思い、何か楽しい課題はないかなぁと考えていて思いついたのがカオス理論でよく取り上げられるロジスティック写像

前にRでプログラム書いた筈だけど見当たらないので、再度作成。といっても今回はこのサイトのものを参照に書いた。更にシンプルにして入力するパラメターの数を半分に。

具体的にはxの初期値を自動的に計算し、また、計算された値の内プロットに用いる値の数を指定ではなく自動的に算出する(最後の三分の一)とした。

Rでウィキペディアにある様なプロットに出来るだけ近づけることが一番大変だった。プロットに関してはコレコレを参照。

コード自体はもっと単純に書ける人はいるだろうけど、JAVAを学ぶ動機のひとつとしてなのでこれでいいかな?

こんな感じに出来上がりました。ただの自己満足です。

rmin <- 2.4 #minimum value of r
rmax <- 4.0 #maximum value of r
r <- seq(rmin, rmax, by = 0.001) # range of r. increment of the sequence = 0.001
N <- 1000 # no. of iterations

logistic.map <- function(r, N){
x<-0
x[1] <- runif(1) #initial value of x
for (i in 1:N){
x[i+1] <- r*x[i]*(1-x[i])
}
x[c((N-N/3):N)] #use only the final 1/3 of iterations
}

x <- sapply(r, logistic.map, N=N)
x <- as.vector(x)
r <-sort(rep(r, N/3+1))

plot(x ~ r, pch=”.”, col=rgb(0,0,0,0.05), xaxs=”i”, yaxs=”i”, xlim=c(rmin, rmax), ylim=c(0,1), axes=FALSE, ylab=””)
axis(1, at=seq(2.4, 4.0, by=0.2)) #axis bottom
axis(2, las=1, at=seq(0.0, 1.0, by=0.2)) #axis left
mtext(“x”, side=2, line=2, las=2) #horizontal “x” label on y axis

BUGS

BUGSのモデルを載せるの忘れてました。。

model
{
for(i in 1:N)
{
lwage[i] ~ dnorm (mu[i], tau)
mu[i] <- b[1] + b[2]*exper[i] + b[3]*expersq[i] + b[4]*educ[i] + b[5]*age[i] + b[6]*kidslt6[i] + b[7]*kidsge6[i] + b[8]*city[i]
}

for(k in 1:8)
{
b[k] ~ dnorm(0, 0.0001)
}

sigma ~ dunif(0, 1000)
tau <-1/(sigma*sigma)
}

回帰分析

久し振りにlinear regressionをやってみた。備忘録として。
Stata、RとWinBUGSで。ベイズ統計についてはいずれ書きたいと思います。まぁ、分析の方が主だから哲学的な所にはあまりこだわってはいないけど。ただ、ベイズ統計の方がpost-estimationがより楽しい点はあるかな。
データはStataのmroz。モデルはここを参照。
*****Stata
set more off
set memory 500m
cd "E:\"
use mroz, clear
reg lwage exper expersq educ age kidslt6 kidsge6 city
*****R
setwd("E:/")
library(foreign)
mroz <-read.dta("mroz.dta")
N <- nrow(mroz)
lwage <- mroz$lwage
exper <- mroz$experexpersq <- mroz$expersq
educ <- mroz$educ
age <- mroz$age
kidslt6 <- mroz$kidslt6
kidsge6 <- mroz$kidsge6
city <- mroz$city
model <- lm(lwage ~ exper+expersq+educ+age+kidslt6+kidsge6+city)
model
*****BUGSlibrary(R2WinBUGS)
data <-list("N", "lwage", "exper", "expersq", "educ", "age", "kidslt6", "kidsge6", "city")
params <- c("b")
inits <- function () {list(b=runif(8))}
BUGSModel <- bugs (data, inits, params, "Model.txt", n.chains=5, n.iter=10000, debug=T)