-
Notifications
You must be signed in to change notification settings - Fork 0
/
ptII_quan_Bayes_RandomWalk_helpfuncs.r
93 lines (83 loc) · 2.78 KB
/
ptII_quan_Bayes_RandomWalk_helpfuncs.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
### (C) 2005-2023 by Leo Guertler
### R-code supplement
### to the book
###
### "Subjektive Ansichten und objektive Betrachtungen"
###
### written by Gürtler & Huber (2023)
###
### All R-code is published under the GPL v3 license:
###
### https://www.gnu.org/licenses/gpl-3.0.en.html
###
### except for 'borrowed' code - see links and references.
### For this R-code the original license of the respective
### authors is valid.
###
### R-code published on
###
### https://osdn.net/projects/mixedmethod-rcode
### https://github.com/abcnorio/mixedmethod-rcode
# file:
# ptII_quan_Bayes_RandomWalk_helpfuncs.r
# location:
# chap. 6 [6.13]
# Marko Chain Monte Carlo Simulationen — MCMC
# HELPER FUNCTIONS
###### function to calculate simple random walk
randomwalk <- function(nsim=1000, seed=NA, p=c(0.5,0.5,0.5), disti=c(10,10,10), initial=c(0,0,0), graph=TRUE, D=2,
phi=30, theta=45, type3d="l", ...)
{
require(scatterplot3d)
if(!is.na(seed)) set.seed(seed)
walk <- matrix(data=NA, nrow=nsim, ncol=4)
colnames(walk) <- c("no","x","y","z")
walk[,"no"] <- 1:nsim
# initial value
walk[1,c("x","y","z")] <- initial
for(i in 2:nsim)
{
ux <- runif(1)
uy <- runif(1)
uz <- runif(1)
udx <- round(runif(1,0,disti[1]))
udy <- round(runif(1,0,disti[2]))
udz <- round(runif(1,0,disti[3]))
ifelse(ux <= p[1], walkx <- +1, walkx <- -1)
ifelse(uy <= p[2], walky <- +1, walky <- -1)
ifelse(uz <= p[3], walkz <- +1, walkz <- -1)
walk[i,"x"] <- walk[i-1,"x"] + walkx*udx
walk[i,"y"] <- walk[i-1,"y"] + walky*udy
walk[i,"z"] <- walk[i-1,"z"] + walkz*udz
}
if(graph)
{
if(D == 3)
{
#
nrz <- nrow(walk)
ncz <- ncol(walk)
jet.colors <- colorRampPalette( c("blue", "green") )
# Generate the desired number of colors from this palette
nbcol <- dim(walk)[1]
color <- jet.colors(nbcol)
# Compute the z-value at the facet centres
z <- outer(walk[,"x"],walk[,"y"], function(a,b) a*b^2)
zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
# Recode facet z-values into color indices
facetcol <- cut(zfacet, nbcol)
# persp(walk[,c("x","y","z")], phi=phi, theta=theta, shade=.1, col=color[facetcol], border="yellow")
scatterplot3d(walk[,c("x","y","z")], col.axis="blue", col.grid="lightblue", main="Random walk (3D)", pch=20,
angle=30, highlight.3d=FALSE, grid=TRUE, color=color, box=FALSE, type=type3d, bg="violetred3")
} else #D = 2
{
plot(walk[,c("x","y")], bty="n", type="l", col="violetred3", pre.plot=grid(), main="Random walk (2D)")
points(initial[1:2], pch=8, cex=1.5, col="blue")
}
}
return(walk)
}
# call:
# walk <- randomwalk()
# walk <- randomwalk(D=3)
########################## END OF FUNCTION