Supported by Dr. Osamu Ogasawara and  providing  . |
|
Last data update: 2014.03.03 |
Adaptive Mixture of Student-t DistributionsDescriptionFunction which performs the fitting of an adaptive mixture of Student-t distributions to approximate a target density through its kernel function UsageAdMit(KERNEL, mu0, Sigma0 = NULL, control = list(), ...) Arguments
DetailsThe argument As a first example, consider the kernel function proposed by Gelman-Meng (1991): k(x1,x2) = exp( 0.5*[A*x1^2*x2^2 + x1^2 + x2^2 - 2*B*x1*x2 - 2*C1*x1 - 2*C2*x2] ) where commonly used values are A=1, B=0, C1=3 and C2=3. A vectorized implementation of this function might be:
GelmanMeng <- function(x, A = 1, B = 0, C1 = 3, C2 = 3, log = TRUE)
{
if (is.vector(x))
x <- matrix(x, nrow = 1)
r <- -.5 * (A * x[,1]^2 * x[,2]^2 + x[,1]^2 + x[,2]^2
- 2 * B * x[,1] * x[,2] - 2 * C1 * x[,1] - 2 * C2 * x[,2])
if (!log)
r <- exp(r)
as.vector(r)
}
This way, we may supply a point (x1,x2)
for As a second example, consider the following (simple) econometric model: y_t ~ i.i.d. N(mu,sigma^2) t=1,...,T where mu is the mean value and sigma is the standard deviation. Our purpose is to estimate theta=(mu,sigma) within a Bayesian framework, based on a vector y of T observations; the kernel thus consists of the product of the prior and the likelihood function. As previously mentioned, the kernel function should be vectorized, i.e., treat a (Nx2) matrix of points theta for which the kernel will be evaluated. Using the common (Jeffreys) prior p(theta)=1/sigma for sigma>0, a vectorized implementation of the kernel function might be:
KERNEL <- function(theta, y, log = TRUE)
{
if (is.vector(theta))
theta <- matrix(theta, nrow = 1)
## sub function which returns the log-kernel for a given
## thetai value (i.e., a given row of theta)
KERNEL_sub <- function(thetai)
{
if (thetai[2] > 0) ## check if sigma>0
{ ## if yes, compute the log-kernel at thetai
r <- - log(thetai[2])
+ sum(dnorm(y, thetai[1], thetai[2], TRUE))
}
else
{ ## if no, returns -Infinity
r <- -Inf
}
as.numeric(r)
}
## 'apply' on the rows of theta (faster than a for loop)
r <- apply(theta, 1, KERNEL_sub)
if (!log)
r <- exp(r)
as.numeric(r)
}
Since this kernel function also depends on the vector y, it
must be passed to To gain even more speed, implementation of The argument
ValueA list with the following components:
where H (>=1) is the number of components in the adaptive
mixture of Student-t distributions and d (>=1) is
the dimension of the first argument in
NoteFurther details and examples of the R package Further details on the core algorithm are given in Hoogerheide (2006), Hoogerheide, Kaashoek, van Dijk (2007) and Hoogerheide, van Dijk (2008). The adaptive mixture Please cite the package in publications. Use Author(s)David Ardia for the R port,
Lennart F. Hoogerheide and Herman K. van Dijk for the ReferencesArdia, D., Hoogerheide, L.F., van Dijk, H.K. (2009a). AdMit: Adaptive Mixture of Student-t Distributions. The R Journal 1(1), pp.25–30. http://journal.r-project.org/2009-1/ Ardia, D., Hoogerheide, L.F., van Dijk, H.K. (2009b). Adaptive Mixture of Student-t Distributions as a Flexible Candidate Distribution for Efficient Simulation: The R Package AdMit. Journal of Statistical Software 29(3), pp.1–32. http://www.jstatsoft.org/v29/i03/ Gelman, A., Meng, X.-L. (1991). A Note on Bivariate Distributions That Are Conditionally Normal. The American Statistician 45(2), pp.125–126. Hoogerheide, L.F. (2006). Essays on Neural Network Sampling Methods and Instrumental Variables. PhD thesis, Tinbergen Institute, Erasmus University Rotterdam (NL). ISBN: 9051708261. (Book nr. 379 of the Tinbergen Institute Research Series.) Hoogerheide, L.F., Kaashoek, J.F., van Dijk, H.K. (2007). On the Shape of Posterior Densities and Credible Sets in Instrumental Variable Regression Models with Reduced Rank: An Application of Flexible Sampling Methods using Neural Networks. Journal of Econometrics 139(1), pp.154–180. doi: 10.1016/j.jeconom.2006.06.009. Hoogerheide, L.F., van Dijk, H.K. (2008). Possibly Ill-Behaved Posteriors in Econometric Models: On the Connection between Model Structures, Non-elliptical Credible Sets and Neural Network Simulation Techniques. Tinbergen Institute discussion paper 2008-036/4. http://www.tinbergen.nl/discussionpapers/08036.pdf See Also
Examples
## NB : Low number of draws for speedup. Consider using more draws!
## Gelman and Meng (1991) kernel function
GelmanMeng <- function(x, A = 1, B = 0, C1 = 3, C2 = 3, log = TRUE)
{
if (is.vector(x))
x <- matrix(x, nrow = 1)
r <- -.5 * (A * x[,1]^2 * x[,2]^2 + x[,1]^2 + x[,2]^2
- 2 * B * x[,1] * x[,2] - 2 * C1 * x[,1] - 2 * C2 * x[,2])
if (!log)
r <- exp(r)
as.vector(r)
}
## Run AdMit (with default values)
set.seed(1234)
outAdMit <- AdMit(GelmanMeng, mu0 = c(0.0, 0.1), control = list(Ns = 1e4))
print(outAdMit)
## Run AdMit (using importance sampling to estimate
## the modes and the scale matrices)
set.seed(1234)
outAdMit <- AdMit(KERNEL = GelmanMeng,
mu0 = c(0.0, 0.1),
control = list(IS = TRUE, Ns = 1e4))
print(outAdMit)
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(AdMit)
Loading required package: mvtnorm
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/AdMit/AdMit.Rd_%03d_medium.png", width=480, height=480)
> ### Name: AdMit
> ### Title: Adaptive Mixture of Student-t Distributions
> ### Aliases: AdMit
> ### Keywords: htest
>
> ### ** Examples
>
> ## NB : Low number of draws for speedup. Consider using more draws!
> ## Gelman and Meng (1991) kernel function
> GelmanMeng <- function(x, A = 1, B = 0, C1 = 3, C2 = 3, log = TRUE)
+ {
+ if (is.vector(x))
+ x <- matrix(x, nrow = 1)
+ r <- -.5 * (A * x[,1]^2 * x[,2]^2 + x[,1]^2 + x[,2]^2
+ - 2 * B * x[,1] * x[,2] - 2 * C1 * x[,1] - 2 * C2 * x[,2])
+ if (!log)
+ r <- exp(r)
+ as.vector(r)
+ }
>
> ## Run AdMit (with default values)
> set.seed(1234)
> outAdMit <- AdMit(GelmanMeng, mu0 = c(0.0, 0.1), control = list(Ns = 1e4))
> print(outAdMit)
$CV
[1] 4.4425871 1.3448258 0.8856950 0.8380885
$mit
$mit$p
cmp1 cmp2 cmp3 cmp4
0.4525493 0.1338627 0.2672455 0.1463425
$mit$mu
k1 k2
cmp1 0.3819661 2.61803340
cmp2 3.8276502 0.20336840
cmp3 1.8030488 1.05601824
cmp4 2.5878805 0.05964092
$mit$Sigma
k1k1 k1k2 k2k1 k2k2
cmp1 0.2291798 -0.40000024 -0.40000024 1.57082073
cmp2 0.8477147 -0.08618971 -0.08618971 0.07277385
cmp3 0.2885380 -0.10007347 -0.10007347 0.21978570
cmp4 0.7388323 -0.17056222 -0.17056222 0.21741696
$mit$df
[1] 1
$summary
H METHOD.mu TIME.mu METHOD.p TIME.p CV
1 1 BFGS 0.003 NONE 0.000 4.4425871
2 2 BFGS 0.008 NLMINB 0.006 1.3448258
3 3 BFGS 0.024 NLMINB 0.016 0.8856950
4 4 BFGS 0.027 NLMINB 0.033 0.8380885
>
> ## Run AdMit (using importance sampling to estimate
> ## the modes and the scale matrices)
> set.seed(1234)
> outAdMit <- AdMit(KERNEL = GelmanMeng,
+ mu0 = c(0.0, 0.1),
+ control = list(IS = TRUE, Ns = 1e4))
> print(outAdMit)
$CV
[1] 4.4425871 1.0830595 0.9504815 0.9039732
$mit
$mit$p
cmp1 cmp2 cmp3 cmp4
0.4772211 0.2456065 0.1478713 0.1293012
$mit$mu
k1 k2
cmp1 0.3819661 2.6180334
cmp2 3.3555320 0.3555571
cmp3 2.1254006 0.1176695
cmp4 2.7367404 0.1234779
$mit$Sigma
k1k1 k1k2 k2k1 k2k2
cmp1 0.2291798 -0.4000002 -0.4000002 1.5708207
cmp2 1.7567372 -0.5057741 -0.5057741 0.3152036
cmp3 0.3314020 -0.4409544 -0.4409544 0.7345853
cmp4 0.5203117 -0.4649378 -0.4649378 0.4992338
$mit$df
[1] 1
$summary
H METHOD.mu TIME.mu METHOD.p TIME.p CV
1 1 BFGS 0.001 NONE 0.000 4.4425871
2 2 IS 0.05-1 0.003 NLMINB 0.004 1.0830595
3 2 IS 0.05-0.25 0.003 NLMINB 0.004 1.1666006
4 2 IS 0.05-4 0.003 NLMINB 0.005 1.3681444
5 2 IS 0.15-1 0.003 NLMINB 0.005 1.1898718
6 2 IS 0.15-0.25 0.003 NLMINB 0.004 1.1663085
7 2 IS 0.15-4 0.003 NLMINB 0.005 1.6186371
8 2 IS 0.3-1 0.003 NLMINB 0.005 1.3583240
9 2 IS 0.3-0.25 0.003 BFGS 0.023 1.2629412
10 2 IS 0.3-4 0.003 NLMINB 0.006 1.9138975
11 3 IS 0.05-1 0.002 NLMINB 0.011 0.9504815
12 3 IS 0.05-0.25 0.002 NLMINB 0.010 0.9885004
13 3 IS 0.05-4 0.002 NLMINB 0.011 1.0051189
14 3 IS 0.15-1 0.002 NLMINB 0.013 0.9844874
15 3 IS 0.15-0.25 0.002 NLMINB 0.014 0.9526902
16 3 IS 0.15-4 0.002 NLMINB 0.043 1.0851996
17 3 IS 0.3-1 0.002 NLMINB 0.013 1.0187545
18 3 IS 0.3-0.25 0.002 NLMINB 0.012 0.9540494
19 3 IS 0.3-4 0.002 NLMINB 0.040 1.0864449
20 4 IS 0.05-1 0.002 NLMINB 0.022 0.9080267
21 4 IS 0.05-0.25 0.002 NLMINB 0.055 0.9063889
22 4 IS 0.05-4 0.002 NLMINB 0.062 0.9522559
23 4 IS 0.15-1 0.002 NLMINB 0.023 0.9456680
24 4 IS 0.15-0.25 0.002 NLMINB 0.025 0.9039732
25 4 IS 0.15-4 0.002 NLMINB 0.050 0.9538416
26 4 IS 0.3-1 0.002 NLMINB 0.049 0.9547578
27 4 IS 0.3-0.25 0.002 NLMINB 0.038 0.9058338
28 4 IS 0.3-4 0.002 NLMINB 0.049 0.9528613
>
>
>
>
>
> dev.off()
null device
1
>
|