Description Usage Arguments Details Value Author(s) References Examples

View source: R/bayesian.multiband.R

`bayesian.multiband`

produces posterior samples of all the model parameters including the time delay between doubly-lensed multi-band light curves.

1 2 3 4 5 6 7 8 9 10 | ```
bayesian.multiband(data.band1, data.band2, n.bands = 2,
theta.ini = c(0.01, 0.01, 100, 100, 0.5),
delta.ini, delta.uniform.range = c(-500, 500), delta.proposal.scale = 1,
tau.proposal.scale = 1, tau.prior.shape = 1, tau.prior.scale = 1,
sigma2.proposal.scale = 0.5, sigma2.prior.shape = 1, sigma2.prior.scale = 1e-7,
rho.proposal.scale = 0.1, beta.prior.diag = 10 * c(0.1, 0.01, 1e-3, 1e-5)^2,
micro = 3, timesc = 100, adaptive.frequency = 100,
adaptive.delta.factor = 0.1, adaptive.tau.factor = 0.1,
adaptive.sigma2.factor = 0.1, adaptive.rho.factor = 0.1,
sample.size, warmingup.size)
``` |

`data.band1` |
A ( |

`data.band2` |
A ( |

`n.bands` |
The number of bands used to obtain the data. For now, only two bands are allowed. The authors plan to modify the code to allow more bands in the near future. |

`theta.ini` |
A vector for initial values of the OU processes for two bands, i.e., ( |

`delta.ini` |
Initial values of the time delay. |

`delta.uniform.range` |
The range of the Uniform prior distribution for the time delay. Default range is set to |

`delta.proposal.scale` |
The proposal scale of the Metropolis step for the time delay. Default is 1. |

`tau.proposal.scale` |
The proposal scale of the Metropolis-Hastings step for |

`tau.prior.shape` |
The shape parameter of the Inverse-Gamma hyper-prior distribution for |

`tau.prior.scale` |
The scale parameter of the Inverse-Gamma hyper-prior distribution for |

`sigma2.proposal.scale` |
The proposal scale of the Metropolis-Hastings step for |

`sigma2.prior.shape` |
The shape parameter of the Inverse-Gamma hyper-prior distribution for |

`sigma2.prior.scale` |
The scale parameter of the Inverse-Gamma hyper-prior distribution for |

`rho.proposal.scale` |
The proposal scale of the Metropolis-Hastings step for |

`beta.prior.diag` |
The diagonal elements of the covariance matrix in the multivariate Gaussian prior for |

`micro` |
A non-negative integer less than or equal to 3. It determines the order of a polynomial regression model that accounts for the long-term trend of microlensing effect. Default is 3. |

`timesc` |
It scales the observation time for fitting polynomial of microlensing, i.e., time / timesc, so that the coefficients are not too small. Default is 100. |

`adaptive.frequency` |
The adaptive MCMC is applied for every specified frequency. If it is specified as 500, the adaptive MCMC is applied to every 500th iterstion. Default is 100. |

`adaptive.delta.factor` |
The factor, exp( |

`adaptive.tau.factor` |
The factor, exp( |

`adaptive.sigma2.factor` |
The factor, exp( |

`adaptive.rho.factor` |
The factor, exp( |

`sample.size` |
The number of the posterior samples of each model parameter. |

`warmingup.size` |
The number of burn-in samples for MCMC. |

The function `bayesian.multiband`

produces posterior samples of the model parameters, where the time delay (delta) is of primary interest. For now, this function only supports doubly-lensed data observed in two bands. The authors plan to generalize this code to account for more than two bands and more than two lens.

Please note that when astronomical time series data are loaded on R by `read.table`

, `read.csv`

, etc., some decimal places of the the observation times are automatically rounded because R's default is to load seven effective digits. For example, R will load the observation time 51075.412789 as 51075.41. This default will produce many ties in observation times even though there is actually no tie in observation times. To prevent this, please type "`options(digits = 11)`

" before loading the data if the observation times are in seven effective digits.

The outcomes of `bayesian.multiband`

are composed of:

`delta` |
A vector for |

`beta` |
An |

`rho` |
A vector for |

`sigma` |
An |

`tau` |
An |

`tau.accept.rate` |
The acceptance rate of the MCMC for |

`sigma.accept.rate` |
The acceptance rate of the MCMC for |

`delta.accept.rate` |
The acceptance rate of the MCMC for the time delay |

`rho.accept.rate` |
The acceptance rate of the MCMC for |

Zhirui Hu and Hyungsuk Tak

Hyungsuk Tak, Kaisey Mandel, David A. van Dyk, Vinay L. Kashyap, Xiao-Li Meng, and Aneta Siemiginowska (2017). "Bayesian Estimates of Astronomical Time Delays between Gravitationally Lensed Stochastic Light Curves," The Annals of Applied Statistics, 11 (3), 1309-1348. Hyungsuk Tak, Xiao-Li Meng, and David A. van Dyk (2018), "A Repelling-Attracting Metropolis Algorithm for Multimodality", Journal of Computational and Graphical Statistics, 27 (3), 479-490. Zhirui Hu and Hyungsuk Tak (2020+), "Modeling Stochastic Variability in Multi-Band Time Series Data," arXiv:2005.08049.

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 | ```
# Loading datasets
data(simple.band1)
data(simple.band2)
# Doubly-lensed quasar data set observed in two bands
# Each data set contains doubly-lensed light curves observed in one band.
head(simple.band1)
head(simple.band2)
# The length of each data set (i.e., number of observation times)
# do not need to be the same.
dim(simple.band1)
dim(simple.band2)
###############################################
# Time delay estimation via Bayesian approach #
###############################################
# Cubic microlensing model (m = 3)
output <- bayesian.multiband(data.band1 = simple.band1,
data.band2 = simple.band2, n.bands = 2,
theta.ini = c(0.01, 0.01, 100, 100, 0.5),
delta.ini = 100, delta.uniform.range = c(-500, 500),
tau.proposal.scale = 1, tau.prior.shape = 1, tau.prior.scale = 1,
sigma2.proposal.scale = 0.5, sigma2.prior.shape = 1, sigma2.prior.scale = 1e-7,
rho.proposal.scale = 0.1, beta.prior.diag = 10 * c(0.1, 0.01, 1e-3, 1e-5)^2,
micro = 3, timesc = 100, adaptive.frequency = 100,
adaptive.delta.factor = 0.1, adaptive.tau.factor = 0.1,
adaptive.sigma2.factor = 0.1, adaptive.rho.factor = 0.1,
sample.size = 100, warmingup.size = 100)
names(output)
# hist(output$delta, 20)
# plot(output$delta, type = "l")
# acf(output$delta)
# output$delta.accept.rate
``` |

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.