Description Usage Arguments Details Author(s) Examples
A class to implement a Sequential Monte Carlo sampler for a dynamic Probit regression model.
1 2 |
formula |
A formula with the response and covariates in the model. See Details. |
mu0 |
Mean vector of the prior distribution over the Probit regression parameters. |
Sigma0 |
Covariance matrix of the prior distribution over the Probit regression parameters. |
Q |
Covariance matrix of the update distribution over the Probit regression parameters. |
N |
Total number of particles. |
data |
An optional |
resampleC |
value between 0 and 1, indicating when to resample |
output |
specifying what is returned. See details. |
The function DynProbit initializes a particle filter for a dynamic Probit regression model, and runs it over the data specified in the formula. The response can be either a factor, in which the first level is taken as “failure” and the remaining ones as “success”, or a binary vector.
The prior distriution over the parameters at $t=1$ is a (multivariate) Normal distribution with mean mu0
and covariance Sigma0
.
If output = means
, a matrix with the posterior means at each time point is returned.
If output = filter
, a list is returned with the particle values and their corresponding weights for each time-point t.
If output = smoother
, a list is returned with the particle paths of the N particles and the weights at t=T. Note that this
is a naive implementation of a particle smoother, and the results are not expected to be good for initial time points. This
option is to be used mainly for didacting reasons.
Maarten Speekenbrink
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 | set.seed(120912)
nt <- 100
x <- cbind(rnorm(nt,sd=1),rnorm(nt,sd=1))
theta <- cbind(
rnorm(1,mean=0,sd=sqrt(.2))+c(0,cumsum(rnorm(nt-1,sd=sqrt(.01)))),
rnorm(1,mean=0,sd=sqrt(.2))+c(0,cumsum(rnorm(nt-1,sd=sqrt(.01)))),
rnorm(1,mean=0,sd=sqrt(.2))+c(0,cumsum(rnorm(nt-1,sd=sqrt(.01))))
)
y <- rbinom(nt,prob=1-pnorm(rowSums(cbind(1,x)*theta)),size=1)
means <- DynProbit(y~x[,1] + x[,2],mu0=c(0,0,0),Sigma0=.2*diag(3),Q=.01*diag(3),N=2000)
## Not run:
layout(matrix(1:3,ncol=3))
plot(c(1,100),c(min(c(means[,1],theta[,1])),max(c(means[,1],theta[,1]))),type="n")
lines(1:100,theta[,1],lty=2)
lines(1:100,means[,1],type="l",lty=1)
plot(c(1,100),c(min(c(means[,2],theta[,2])),max(c(means[,2],theta[,2]))),type="n")
lines(1:100,theta[,2],lty=2)
lines(1:100,means[,2],type="l",lty=1)
plot(c(1,100),c(min(c(means[,3],theta[,3])),max(c(means[,3],theta[,3]))),type="n")
lines(1:100,theta[,3],lty=2)
lines(1:100,means[,3],type="l",lty=1)
## End(Not run)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.