## PROGRAMS THE GIBBS SAMPLING - pag 64-67
# Define the number of observations & simulations
n <- 30
# Load the data
#y <- rnorm(n,4,10) # This would generate random values
y <- c(1.2697,7.7637,2.2532,3.4557,4.1776,6.4320,-3.6623,7.7567,5.9032,7.2671,
-2.3447,8.0160,3.5013,2.8495,0.6467,3.2371,5.8573,-3.3749,4.1507,4.3092,
11.7327,2.6174,9.4942,-2.7639,-1.5859,3.6986,2.4544,-0.3294,0.2329,5.2846)
# Defines the hyper-parameters to build the full conditionals
ybar <- mean(y)
mu_0 <- 0
sigma2_0 <- 10000
alpha_0 <- 0.01
beta_0 <- 0.01
# Initialises the parameters
mu <- tau <- numeric()
sigma2 <- 1/tau
mu[1] <- rnorm(1,0,3)
tau[1] <- runif(1,0,3)
sigma2[1] <- 1/tau[1]
# Gibbs sampling (samples from the full conditionals)
nsim <- 1000
for (i in 2:nsim) {
sigma_n <- sqrt(1/(1/sigma2_0 + n/sigma2[i-1]))
mu_n <- (mu_0/sigma2_0 + n*ybar/sigma2[i-1])*sigma_n^2
mu[i] <- rnorm(1,mu_n,sigma_n)
alpha_n <- alpha_0+n/2
beta_n <- beta_0 + sum((y-mu[i])^2)/2
tau[i] <- rgamma(1,alpha_n,beta_n)
sigma2[i] <- 1/tau[i]
}
## Creates a bivariate contour, using the package mvtnorm
require(mvtnorm)
theta <- c(mean(mu),mean(sqrt(sigma2)))
sigma <- c(var(mu),var(sqrt(sigma2)))
rho <- cor(mu,sqrt(sigma2))
ins <- c(-10,10)
x1 <- seq(theta[1]-abs(ins[1]),theta[1]+abs(ins[1]),length.out=1000)
x2 <- seq(theta[2]-abs(ins[2]),theta[2]+abs(ins[2]),length.out=1000)
all <- expand.grid(x1, x2)
Sigma<-matrix(c(sigma[1], sigma[1]*sigma[2]*rho, sigma[1]*sigma[2]*rho, sigma[2]), nrow=2, ncol=2)
f.x <- dmvnorm(all, mean = theta, sigma = Sigma)
f.x2 <- matrix(f.x, nrow=length(x1), ncol=length(x2))
##Plots of the results at different simulation lengths
lw <- c(1,1.6)
# First 100 iterations
plot(mu,sqrt(sigma2),col="white",pch=20,cex=.3,xlim=range(mu),ylim=range(sqrt(sigma2)),xlab=expression(mu),ylab=expression(sigma),
main="After 100 iterations")
for (i in 1:99) {
points(c(mu[i],mu[i+1]),c(sqrt(sigma2)[i],sqrt(sigma2)[i]),t="l",col="grey60")
points(c(mu[i+1],mu[i+1]),c(sqrt(sigma2)[i],sqrt(sigma2)[i+1]),t="l",col="grey60")
}
text(jitter(mu[1:100]),jitter(sqrt(sigma2[1:100])),1:100,col="grey20")
contour(x = x1, y = x2, z = f.x2, nlevels=5,add=T,col="black",drawlabels=FALSE,lwd=lw[1])
=============================================================
똑같은 걸 이번에는 jags로
=============================
library('rjags')
y <- c(1.2697,7.7637,2.2532,3.4557,4.1776,6.4320,-3.6623,7.7567,5.9032,7.2671,
-2.3447,8.0160,3.5013,2.8495,0.6467,3.2371,5.8573,-3.3749,4.1507,4.3092,
11.7327,2.6174,9.4942,-2.7639,-1.5859,3.6986,2.4544,-0.3294,0.2329,5.2846)
model_string <- "model {
for (i in 1:length(y)) {
y[i] ~ dnorm(mu, tau)
}
mu ~ dnorm(0,0.0001)
sigma ~ dlnorm(0,0.0625)
tau <- 1/ pow(sigma, 2)
}"
model <- jags.model(textConnection(model_string),
data = list(y = y),
n.chains = 3,
n.adapt = 100)
update(jags, 100)
output<-coda.samples(model,
variable.names=c('mu', 'sigma'),
n.iter=100)
plot(output)
2013년 11월 21일 목요일
2013년 11월 18일 월요일
2013년 11월 15일 금요일
installation of perl modules
# cpan cpan shell -- CPAN exploration and modules installation (v1.9205) ReadLine support available (maybe install Bundle::CPAN or Bundle::CPANxxl?) cpan[1]> install "Email::Reply";
2013년 11월 7일 목요일
2013년 11월 6일 수요일
2013년 10월 31일 목요일
link
section .data
ldlibry dd 0
pro dd 0
dll db "user32.dll",0
myFtion db "MessageBoxA",0
MSG db "ASM GetProcAddress",0
extern _LoadLibraryA@4
extern _FreeLibrary@4
extern _GetProcAddress@8
extern _ExitProcess@4
global _start
section .text
_start:
push dll ;push user32.dll
call _LoadLibraryA@4 ;Call the API.
mov [ldlibry],eax ;eax hold return address. So eax=LoadLibrary("user32.dll") and now ldlibry=LoadLibrary("user32.dll")
;now we need to call GetProcAddress
push myFtion ;The API name we are going to call
push eax ;LoadLibrary("user32.dll")
call _GetProcAddress@8 ;GetProcAddress(LoadLibrary("user32.dll"),"MessageBoxA"). Again eax holding the return address
push 0x0 ;MB_OK
push MSG ;TITLE="ASM GetProcAddress"
push MSG ;Messgage="ASM GetProcAddress"
push 0 ;Reserved=0
call eax ;Call MessageBoxA through GetProcAddress.
push dword [ldlibry] ; ldlibry holding the LoadLibrary("user32.dll"). Again load to Free up.
call _FreeLibrary@4 ;Call the Windows api FreeLibrary()
;We should exit the process otherwise it may cause "access violation"
push 0 ;load 0 to stack
call _ExitProcess@4 ;Call ExitProcess
----------------
C:\work>nasm -f win32 asmcall.asm
C:\work>link asmcall.obj /entry:start /subsystem:console kernel32.lib
Microsoft (R) Incremental Linker Version 6.00.8168
Copyright (C) Microsoft Corp 1992-1998. All rights reserved.
피드 구독하기:
글 (Atom)