gibbs: Gibbs sampling with Metropolis steps and multivariate...
In gibbs.met: Naive Gibbs Sampling with Metropolis Steps
Description
Usage
Arguments
Value
Examples
Description
The function gibbs_met performs Gibbs sampling with each 1-dimensional distribution sampled with Metropolis update using Gaussian proposal distribution centered at the previous state. The function met_gaussian updates the whole state with Metropolis method using independent Gaussian proposal distribution centered at the previous state. The sampling is carried out without considering any special tricks for improving efficiency. The functions are written for routine applications in moderate-dimensional problems.
Usage
1
2
3
4
gibbs_met(log_f,no_var,ini_value,
iters,iters_per.iter=1,iters_met,stepsizes_met, ...)
met_gaussian(log_f, no_var, ini_value,
iters,iters_per.iter=1,stepsizes_met, ...)
Arguments
log_f
the log of the density function from which one wants to sample.
no_var
the number of variables to be sampled.
ini_value
the initial value.
iters,iters_per.iter
Run iters super-transition, each consisting of iters_per.iter single Markov chain update (using Gibbs sampling or Metropolis sampling). Only the state at the end of each super-transition is saved and returned.
iters_met
iterations of updating each 1-dim conditional distribution with Metropolis method in Gibbs sampling using the function gibbs_met.
stepsizes_met
a vector of length no_var, with stepsizes_met[i] being the standard deviation of Gaussian proposal for updating 'i'th variable, which is used either in Gibbs sampling for sampling from the 'i'th conditional distribution, or used in multivariate Metropolis sampling (using the function met_gaussian) for 'i'th variable.
...
extra arguments needed to compute log_f.
Value
a matrix with dim (iters +1) * no_var is returned, with each row for an iteration
Examples
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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
################################################################################
## demonstration by sampling from bivariate normal distributions
################################################################################
## the function computing the log density function of multivariate normal
## x --- a vector, the p.d.f at x will be computed
## mu --- the mean vector of multivariate normal distribution
## A --- the inverse covariance matrix of multivariate normal distribution
log_pdf_mnormal <- function(x, mu, A)
{
0.5 * (-length(mu)*log(2*pi)+sum(log(svd(A)$d))-sum(t(A*(x-mu))*(x-mu)) )
}
## sampling from a bivariate normal distribution with correlation 0.1,
## both marginal standard deviations 1, mean vector (0,5)
A <- solve(matrix(c(1,0.1,0.1,1),2,2))
mc_mvn <- gibbs_met(log_f=log_pdf_mnormal,no_var=2,
ini_value=c(0,0),iters=400,iters_met=2,
stepsizes_met=c(0.5,0.5), mu=c(0,5), A = A)
postscript("mc_mvn_lowcor.eps",width=7,height=8,horiz=FALSE)
par(mfrow=c(2,2), oma=c(0,0,1,0))
## looking at the trace of Markov chain in the first 100 iterations
plot(mc_mvn[1:100,1],mc_mvn[1:100,2],type="b",pch=20,
main="Markov chain trace of both variables")
## looking at the trace of Markov chain for a variable
plot(mc_mvn[,1],type="b",pch=20, main="Markov chain trace of the 1st variable")
## looking at the QQ plot of the samples for a variable
qqnorm(mc_mvn[-(1:50),1],main="Normal QQ plot of the 1st variable")
## looking at the ACF of the samples for a variable
acf(mc_mvn[,1],main="ACF plot of the 1st variable")
title(main="Gibbs sampling for a bivariate normal with correlation 0.1",
outer=TRUE)
dev.off()
## checking the correlation of samples
cat("The sample correlation is",cor(mc_mvn[-(1:50),1],mc_mvn[-(1:50),2]),"\n")
################################################################################
## demonstration by sampling from a mixture bivariate normal distribution
################################################################################
## the function computing the log density function of mixture multivariate normal
## x --- a vector, the p.d.f at x will be computed
## mu1,mu2 --- the mean vectors of multivariate normal distributions
## A1,A2 --- the inverse covariance matrice of multivariate normal distributions
## mixture proportion is 0.5
log_pdf_twonormal <- function(x, mu1, A1, mu2, A2)
{ log_sum_exp(c(log_pdf_mnormal(x,mu1,A1),log_pdf_mnormal(x,mu2,A2))
)
}
log_sum_exp <- function(lx)
{ ml <- max(lx)
ml + log(sum(exp(lx-ml)))
}
## set parameters defining a mixture bivariate distribution
A1 <- solve(matrix(c(1,0.1,0.1,1),2,2))
A2 <- solve(matrix(c(1,0.1,0.1,1),2,2))
mu1 <- c(0,0)
mu2 <-c(4,4)
## performing Gibbs sampling
mc_mvn <- gibbs_met(log_f=log_pdf_twonormal,no_var=2,ini_value=c(0,0),
iters=400,iters_met=2,stepsizes_met=c(0.5,0.5),
mu1=mu1,mu2=mu2,A1=A1,A2=A2)
postscript("mc_mvn_closemix.eps",width=7,height=8,horiz=FALSE)
par(mfrow=c(2,2), oma=c(0,0,2,0))
## looking at the trace of Markov chain in the first 100 iterations
plot(mc_mvn[,1],mc_mvn[,2],type="b",pch=20,
main="Markov chain trace of both variables")
## looking at the trace of Markov chain for a variable
plot(mc_mvn[,1],type="b",pch=20, main="Markov chain trace of the 1st variable")
## looking at the trace of Markov chain for a variable
plot(mc_mvn[,2],type="b",pch=20, main="Markov chain trace of the 2nd variable")
## looking at the ACF of the samples for a variable
acf(mc_mvn[,1],main="ACF plot of the 1st variable")
title(main="Gibbs sampling for a mixture of two bivariate normal distributions
with locations (0,0) and (4,4)", outer=TRUE)
dev.off()
## checking the correlation of samples
cat("The sample correlation is",cor(mc_mvn[-(1:50),1],mc_mvn[-(1:50),2]),"\n")
###############################################################################
## Sampling from a mixture bivariate normal distribution with Metropolis method
###############################################################################
## set parameters defining a mixture bivariate distribution
A1 <- solve(matrix(c(1,0.1,0.1,1),2,2))
A2 <- solve(matrix(c(1,0.1,0.1,1),2,2))
mu1 <- c(0,0)
mu2 <-c(6,6)
## performing Gibbs sampling
mc_mvn <- met_gaussian(log_f=log_pdf_twonormal,no_var=2,ini_value=c(0,0),
iters=400,iters_per.iter=2,stepsizes=c(1,1),
mu1=mu1,mu2=mu2,A1=A1,A2=A2)
postscript("mc_mvn_farmix_met.eps",width=7,height=8,horiz=FALSE)
par(mfrow=c(2,2), oma=c(0,0,2,0))
## looking at the trace of Markov chain in the first 100 iterations
plot(mc_mvn[,1],mc_mvn[,2],type="b",pch=20,
main="Markov chain trace of both variables")
## looking at the trace of Markov chain for a variable
plot(mc_mvn[,1],type="b",pch=20, main="Markov chain trace of the 1st variable")
## looking at the trace of Markov chain for a variable
plot(mc_mvn[,2],type="b",pch=20, main="Markov chain trace of the 2nd variable")
## looking at the ACF of the samples for a variable
acf(mc_mvn[,1],main="ACF plot of the 1st variable")
title(main="Sampling with Metropolis method for a mixture of two bivariate normal
distributions with locations (0,0) and (6,6)", outer=TRUE)
dev.off()
## checking the correlation of samples
cat("The sample correlation is",cor(mc_mvn[-(1:50),1],mc_mvn[-(1:50),2]),"\n")
gibbs.met documentation built on May 2, 2019, 1:46 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Description Usage Arguments Value Examples
Description
The function gibbs_met performs Gibbs sampling with each 1-dimensional distribution sampled with Metropolis update using Gaussian proposal distribution centered at the previous state. The function met_gaussian updates the whole state with Metropolis method using independent Gaussian proposal distribution centered at the previous state. The sampling is carried out without considering any special tricks for improving efficiency. The functions are written for routine applications in moderate-dimensional problems.
Usage
1 2 3 4 | gibbs_met(log_f,no_var,ini_value,
iters,iters_per.iter=1,iters_met,stepsizes_met, ...)
met_gaussian(log_f, no_var, ini_value,
iters,iters_per.iter=1,stepsizes_met, ...)
|
Arguments
log_f |
the log of the density function from which one wants to sample. |
no_var |
the number of variables to be sampled. |
ini_value |
the initial value. |
iters,iters_per.iter |
Run |
iters_met |
iterations of updating each 1-dim conditional distribution with Metropolis method in Gibbs sampling using the function |
stepsizes_met |
a vector of length |
... |
extra arguments needed to compute |
Value
a matrix with dim (iters +1) * no_var is returned, with each row for an iteration
Examples
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 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 | ################################################################################
## demonstration by sampling from bivariate normal distributions
################################################################################
## the function computing the log density function of multivariate normal
## x --- a vector, the p.d.f at x will be computed
## mu --- the mean vector of multivariate normal distribution
## A --- the inverse covariance matrix of multivariate normal distribution
log_pdf_mnormal <- function(x, mu, A)
{
0.5 * (-length(mu)*log(2*pi)+sum(log(svd(A)$d))-sum(t(A*(x-mu))*(x-mu)) )
}
## sampling from a bivariate normal distribution with correlation 0.1,
## both marginal standard deviations 1, mean vector (0,5)
A <- solve(matrix(c(1,0.1,0.1,1),2,2))
mc_mvn <- gibbs_met(log_f=log_pdf_mnormal,no_var=2,
ini_value=c(0,0),iters=400,iters_met=2,
stepsizes_met=c(0.5,0.5), mu=c(0,5), A = A)
postscript("mc_mvn_lowcor.eps",width=7,height=8,horiz=FALSE)
par(mfrow=c(2,2), oma=c(0,0,1,0))
## looking at the trace of Markov chain in the first 100 iterations
plot(mc_mvn[1:100,1],mc_mvn[1:100,2],type="b",pch=20,
main="Markov chain trace of both variables")
## looking at the trace of Markov chain for a variable
plot(mc_mvn[,1],type="b",pch=20, main="Markov chain trace of the 1st variable")
## looking at the QQ plot of the samples for a variable
qqnorm(mc_mvn[-(1:50),1],main="Normal QQ plot of the 1st variable")
## looking at the ACF of the samples for a variable
acf(mc_mvn[,1],main="ACF plot of the 1st variable")
title(main="Gibbs sampling for a bivariate normal with correlation 0.1",
outer=TRUE)
dev.off()
## checking the correlation of samples
cat("The sample correlation is",cor(mc_mvn[-(1:50),1],mc_mvn[-(1:50),2]),"\n")
################################################################################
## demonstration by sampling from a mixture bivariate normal distribution
################################################################################
## the function computing the log density function of mixture multivariate normal
## x --- a vector, the p.d.f at x will be computed
## mu1,mu2 --- the mean vectors of multivariate normal distributions
## A1,A2 --- the inverse covariance matrice of multivariate normal distributions
## mixture proportion is 0.5
log_pdf_twonormal <- function(x, mu1, A1, mu2, A2)
{ log_sum_exp(c(log_pdf_mnormal(x,mu1,A1),log_pdf_mnormal(x,mu2,A2))
)
}
log_sum_exp <- function(lx)
{ ml <- max(lx)
ml + log(sum(exp(lx-ml)))
}
## set parameters defining a mixture bivariate distribution
A1 <- solve(matrix(c(1,0.1,0.1,1),2,2))
A2 <- solve(matrix(c(1,0.1,0.1,1),2,2))
mu1 <- c(0,0)
mu2 <-c(4,4)
## performing Gibbs sampling
mc_mvn <- gibbs_met(log_f=log_pdf_twonormal,no_var=2,ini_value=c(0,0),
iters=400,iters_met=2,stepsizes_met=c(0.5,0.5),
mu1=mu1,mu2=mu2,A1=A1,A2=A2)
postscript("mc_mvn_closemix.eps",width=7,height=8,horiz=FALSE)
par(mfrow=c(2,2), oma=c(0,0,2,0))
## looking at the trace of Markov chain in the first 100 iterations
plot(mc_mvn[,1],mc_mvn[,2],type="b",pch=20,
main="Markov chain trace of both variables")
## looking at the trace of Markov chain for a variable
plot(mc_mvn[,1],type="b",pch=20, main="Markov chain trace of the 1st variable")
## looking at the trace of Markov chain for a variable
plot(mc_mvn[,2],type="b",pch=20, main="Markov chain trace of the 2nd variable")
## looking at the ACF of the samples for a variable
acf(mc_mvn[,1],main="ACF plot of the 1st variable")
title(main="Gibbs sampling for a mixture of two bivariate normal distributions
with locations (0,0) and (4,4)", outer=TRUE)
dev.off()
## checking the correlation of samples
cat("The sample correlation is",cor(mc_mvn[-(1:50),1],mc_mvn[-(1:50),2]),"\n")
###############################################################################
## Sampling from a mixture bivariate normal distribution with Metropolis method
###############################################################################
## set parameters defining a mixture bivariate distribution
A1 <- solve(matrix(c(1,0.1,0.1,1),2,2))
A2 <- solve(matrix(c(1,0.1,0.1,1),2,2))
mu1 <- c(0,0)
mu2 <-c(6,6)
## performing Gibbs sampling
mc_mvn <- met_gaussian(log_f=log_pdf_twonormal,no_var=2,ini_value=c(0,0),
iters=400,iters_per.iter=2,stepsizes=c(1,1),
mu1=mu1,mu2=mu2,A1=A1,A2=A2)
postscript("mc_mvn_farmix_met.eps",width=7,height=8,horiz=FALSE)
par(mfrow=c(2,2), oma=c(0,0,2,0))
## looking at the trace of Markov chain in the first 100 iterations
plot(mc_mvn[,1],mc_mvn[,2],type="b",pch=20,
main="Markov chain trace of both variables")
## looking at the trace of Markov chain for a variable
plot(mc_mvn[,1],type="b",pch=20, main="Markov chain trace of the 1st variable")
## looking at the trace of Markov chain for a variable
plot(mc_mvn[,2],type="b",pch=20, main="Markov chain trace of the 2nd variable")
## looking at the ACF of the samples for a variable
acf(mc_mvn[,1],main="ACF plot of the 1st variable")
title(main="Sampling with Metropolis method for a mixture of two bivariate normal
distributions with locations (0,0) and (6,6)", outer=TRUE)
dev.off()
## checking the correlation of samples
cat("The sample correlation is",cor(mc_mvn[-(1:50),1],mc_mvn[-(1:50),2]),"\n")
|
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.