gpfso | R Documentation |
This function implements the G-PFSO (Global Particle Filter Stochastic Optimization) algorithm of \insertCitegerber2020online2;textualPFoptim for minimzing either the function θ\mapsto E[\mathrm{fn}(θ,Y)] from i.i.d. realizations y_1,...,y_n of Y or the function θ\mapsto∑_{i=1}^n \mathrm{fn}(θ,y_i), where θ is a vector of dimension d.
gpfso(obs, N, fn, init, numit, resampling=c("SSP", "STRAT", "MULTI"), ..., control= list())
obs |
Either a vector of observations or a matrix of observations (the number of rows being the sample size). |
N |
Number of particles. The parameter |
fn |
function for a single observation. If theta is an |
init |
Either a vector of size d or a function used to sample the initial particles such that |
numit |
Number of iterations of the algorithm. If |
resampling |
Resampling algorithm to be used. Resamping should be either "SSP" (SSP resampling), "STRAT" (stratified resampling) or "MULTI" (multinomial resampling). |
... |
Further arguments to be passed to |
control |
A |
Note that arguments after ...
must be matched exactly.
G-PFSO computes two estimators of the minimizer of the objective function, namely the estimators \bar{θ}^N_{\mathrm{numit}} and \tilde{θ}^N_{\mathrm{numit}}. The former is defined by \bar{θ}^N_{\mathrm{numit}}=\frac{1}{\mathrm{numit}}∑_{t=1}^{\mathrm{numit}}\tilde{θ}_t^N and converges to a particular element of the search space at a faster rate than the latter, but the latter estimator can find more quickly a small neighborhood of the minimizer of the objective function.
By default the sequence (t_p)_{p≥ 0} is taken as
t_p=t_{p-1}+\lceil \max\big( A t_{p-1}^{\varrho}\log(t_{p-1}),B\big) \rceil
with A=B=1, \varrho=0.1 and t_0=5. The value of A,B,\varrho and t_0 can be changed using the control argument (see below).
The control
argument is a list that can supply any of the following components:
Variance parameter of the distribution used by default to sample the initial particles.
Parameter α of the learning rate t^{-α}, which must be a strictly positive real number. By default, alpha=0.5
.
Scale matrix used to sample the particles. Sigma
must be either a d by d covariance matrix or a strictly positive real number. In this latter case the scale matrix used to sample the particles is diag(Sigma , d )
. By default, Sigma=1
.
If trace=TRUE then the value of \tilde{θ}_{t} and of the effective sample size ESS_t for all t=1,…,\mathrm{numit} are returned. By default, trace=FALSE.
If indep=TRUE and Sigma
is a diagonal matrix or a scalar then the Student's t-distributions have independent components. By default, indep=FALSE and if Sigma
is a not a diagonal matrix this parameter is ignored.
Parameter A of the sequence (t_p)_{p≥q 0} used by default (see above). This parameter must be strictly positive.
Parameter B of the sequence (t_p)_{p≥q 0} used by default (see above). This parameter must non-negative.
Parameter varrho of the sequence (t_p)_{p≥q 0} used by default (see above). This parameter must be in the interval (0,1).
Parameter t_0 of the the sequence (t_p)_{p≥q 0} used by default (see above). This parameter must be a non-negative integer.
Number of degrees of freedom of the Student's t-distributions used at time t\in(t_p)_{≥ 0} to generate the new particles. By default nu=10
A resamling step is performed when ESS_t<=N c_\mathrm{ess}. This parameter must be in the interval (0,1] and by default c_ess
=0.7.
A list with the following components:
B_par |
Value of \bar{θ}^N_{\mathrm{numit}} |
T_par |
Value of \tilde{θ}^N_{\mathrm{numit}} |
T_hist |
Value of \tilde{θ}^N_{t} for t=1,...,\mathrm{numit} (only if trace=TRUE) |
ESS |
Value of the effective sample for t=1,...,\mathrm{numit} (only if trace=TRUE) |
#Definition of fn fn_toy<-function(theta, obs){ test<-rep(0,nrow(theta)) test[theta[,2]>0]<-1 ll<-rep(-Inf,nrow(theta)) ll[test==1]<-dnorm(obs,mean=theta[test==1,1], sd=theta[test==1,2],log=TRUE) return(-ll) } #Generate data y_1,...,y_n n<-10000 #sample size theta_star<-c(0,1) #true parameter value y<-rnorm(n,mean=theta_star[1], sd=theta_star[2]) d<-length(theta_star) #Define init funciton to be used pi0<-function(N){ return(cbind(rnorm(N,0,5), rexp(N))) } ##Example 1: Maximum likelihood estimation in the Gaussian model ##true value of the MLE mle<-c(mean(y),sd(y)) ## use gpfso to compute the MLE Est<-gpfso(y, N=100, fn=fn_toy, init=pi0, numit=20000, control=list(trace=TRUE)) ## print \bar{\theta}^N_{numit} and \tilde{\theta}^N_{numit} print(Est$B_par) print(Est$T_par) ##assess convergence par(mfrow=c(1,2)) for(k in 1:2){ plot(Est$T_hist[,k],type='l', xlab="iteration", ylab="estimated value") lines(cumsum(Est$T_hist[,k])/1:length(Est$T_hist[,k]),type='l', col='red') abline(h=mle[k]) } ##Example 2: Expected log-likelihood estimation in the Gaussian model ## Estimation of theta_star using gpfso Est<-gpfso(y, N=100, fn=fn_toy, init=pi0, control=list(trace=TRUE)) ## print \bar{\theta}^N_{numit} and \tilde{\theta}^N_{numit} print(Est$B_par) print(Est$T_par) ##assess convergence par(mfrow=c(1,2)) for(k in 1:2){ plot(Est$T_hist[,k],type='l', xlab="iteration", ylab="estimated value") lines(cumsum(Est$T_hist[,k])/1:length(Est$T_hist[,k]),type='l', col='red') abline(h=theta_star[k]) }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.