Fits a statistical models to the data, using the robust procedure based on maximum mean discrepancy (MMD) minimization introduced and studied in briol2019statistical,MMD2;textualregMMD.
mmd_est(x, model, par1, par2, kernel, bdwth, control= list())MMD_est returns an object of class
"estMMD".
The functions summary can be used to print the results.
An object of class estMMD is a list containing the following components:
Model estimated
In models where the first parameter is fixed, this is the value par1 fixed by the user. In models where the first parameter is estimated, this is the initialization of the optimization procedure
In models where the second parameter is fixed, this is the value par2 fixed by the user. In models where the second parameter is estimated, this is the initialization of the optimization procedure
Kernel used in the MMD
Bandwidth used. That is, either the value specified by the user, either the bandwidth computedby the median heuristic
Number of steps in the burnin of the GD or SGD algorithm
Number of steps in the GD or SGD algorithm
Stepize parameter in GD or SGD
Parameter used in adagrad to avoid numerical errors in the computation of the step-size
Optimization method used
Error message (if any)
Estimated parameter(s)
The full trajectory of the optimization algorithm (GD or SGD)
Takes the value "est"
Data. Must be a vector for univariate models, a matrix of dimension n by d, where n is the sample size and d the dimension of the model.
Parametric model to be fitted to the data. No default. See details for the list of available models.
First parameter of the model. In models where the first parameter is fixed, it is necessary to provide a value for par1. In models where the first parameter is estimated, par1 can be used to provide an alternative to the default initialization of the optimization algorithms.
Second parameter of the model (if any). In models where the second parameter is fixed, it is necessary to provide a value for par2. In models where the first parameter is estimated,par2 can be used to provide an alternative to the default initialization of the optimization algorithms.
Kernel to be used in the MMD. Available options for kernel are "Gaussian" (Gaussian kernel), "Laplace" (Laplace, or exponential, kernel) and "Cauchy" (Cauchy kernel). By default, kernel="Gaussian"
Bandwidth parameter for the kernel. bdwth must be a strictly positive real number. By default, the value of bdwth is chosen using the median heuristic garreau2017largeregMMD.
A list of control parameters for the numerical optimization of the objective function. See details.
Available options for model are:
beta"Beta distribution with pdf \(~x^{a-1}(1-x)^(b-1)\) on \([0,1]\), par1\(=a\) and par2\(=b\) are both estimated.
binomial"Binomial distribution with pmf \(~p^{x}(1-p)^{N-x}\) on \(\{0,1,...,N\}\), par1\(=N\) and par2\(=p\) are both estimated. Note that in this case, if the user specifies a value for \(N\), it is used as an upper bound rather than an initialization.
binomial.prob"Binomial distribution with pmf \(~p^{x}(1-p)^{N-x}\) on \(\{0,1,...,N\}\), par1\(=N\) is fixed and must be specified by the user while par2\(=p\) is estimated.
binomial.size"Binomial distribution with pmf \(~p^{x}(1-p)^{N-x}\) on \(\{0,1,...,N\}\), par1\(=N\) is estimated while par2\(=p\) fixed and must be specified by the user. Note that in this case, if the user specifies a value for \(N\), it is used as an upper bound rather than an initialization.
Cauchy"Cauchy distribution with pdf \(~1/(1+(x-m)^2)\), par1\(=m\) is estimated.
continuous.uniform.loc"Uniform distribution with pdf \(~1\) on \([m-L/2,m+L/2]\), par1\(=m\) is estimated while par2\(=L\) is fixed and must be specified by the user.
continuous.uniform.upper"Uniform distribution with pdf \(~1\) on \([a,b]\), par1\(=a\) is fixed and must be specified by the user while par2\(=b\) is estimated.
continuous.uniform.lower.upper"Uniform distribution with pdf \(~1\) on \([a,b]\), par1\(=a\) and par2\(=b\) are estimated.
Dirac"Dirac mass at point \(a\) on the reals, par1\(=a\) is estimated.
discrete.uniform"Uniform distribution with pmf \(~1\) on \(\{1,2,..,M\}\), par1\(=M\) is estimated. Note that in this case, if the user specifies a value for \(M\), it is used as an upper bound rather than an initialization.
exponential"Exponential distribution with pdf \(~\exp(-b x)\) on positive reals \(R_+\), par1\(=b\) is estimated.
gamma"Gamma distribution with pdf \(~x^{a-1}\exp(-b x)\) on positive reals \(R_+\), par1\(=a>=0.5\) and par2\(=b\) are estimated.
gamma.shape"Gamma distribution with pdf \(~x^{a-1}\exp(-b x)\) on positive reals \(R_+\), par1\(=a>=0.5\) is estimated while par2\(=b\) is fixed and must be specified by the user.
gamma.rate"Gamma distribution with pdf \(~x^{a-1}\exp(-b x)\) on positive reals \(R_+\), par1\(=a>=0.5\) is fixed and must be specified by the user while par2\(=b\) is estimated.
Gaussian"Gaussian distribution with pdf\(~\exp(-(x-m)^2/2s^2)\) on reals \(R\), par1\(=m\) and par2\(=s\) are estimated.
Gaussian.loc"Gaussian distribution with pdf \(~\exp(-(x-m)^2/2s^2)\) on reals \(R\), par1\(=m\) is estimated while par2\(=s\) is fixed and must be specified by the user.
Gaussian.scale"Gaussian distribution with pdf \(~\exp(-(x-m)^2/2s^2)\) on reals \(R\), par1\(=m\) is fixed and must be specified by the user while par2\(=s\) is estimated.
geometric"Geometric distribution with pmf \(~p(1-p)^x\) on \(\{0,1,2,...\}\), par1\(=p\) is estimated.
multidim.Dirac"Dirac mass at point \(a\) on \(R^d\), par1\(=a\) (\(d\)-dimensional vector) is estimated.
multidim.Gaussian"Gaussian distribution with pdf \(~\exp(-(x-m)'U'U(x-m)\) on \(R^d\), par1\(=m\) (\(d\)-dimensional vector) and par2\(=U\) (\(d\)-\(d\) matrix) are estimated.
multidim.Gaussian.loc"Gaussian distribution with pdf \(~\exp(-\|x-m\|^2/2s^2)\) on \(R^d\), par1\(=m\) (\(d\)-dimensional vector) is estimated while par2\(=s\) is fixed.
multidim.Gaussian.scale"Gaussian distribution with pdf \(~\exp(-(x-m)'U'U(x-m)\) on \(R^d\), par1\(=m\) (\(d\)-dimensional vector) is fixed and must be specified by the user while and par2\(=U\) (\(d\)-\(d\) matrix) is estimated.
Pareto"Pareto distribution with pmf \(~1/x^{a+1}\) on the reals \(>1\), par1\(=a\) is estimated.
Poisson"Poisson distribution with pmf \(~b^x/x!\) on \(\{0,1,2,...\}\), par1\(=b\) is estimated.
The control argument is a list that can supply any of the following components:
Length of the burn-in period in GD or SGD. burnin must be a non-negative integer and defaut burnin==\(500\).
Number of iterations performed after the burn-in period in GD or SGD. nsetps must be an integer strictly larger than 2 and by default nsteps=\(1000\)
Stepsize parameter. An adaptive gradient step is used (adagrad), but it is possible to pre-multiply it by stepsize. It must be strictly positive number and by default stepsize=\(1\)
Parameter used in adagrad to avoid numerical errors in the computation of the step-size. epsilon must be a strictly positive real number and by default epsilon=\(10^{-4}\).
Optimization method to be used: "EXACT" for exact, "GD" for gradient descent and "SGD" for stochastic gradient descent. Not all methods are available for all models. By default, exact is preferred to GD which is prefered to SGD.
Additional details and examples in the package publication MMD-packageregMMD.
#simulate data
x = rnorm(50,0,1.5)
# estimate the mean and variance (assuming the data is Gaussian)
Est = mmd_est(x, model="Gaussian")
# print a summary
summary(Est)
# estimate the mean (assuming the data is Gaussian with known standard deviation =1.5)
Est2 = mmd_est(x, model="Gaussian.loc", par2=1.5)
# print a summary
summary(Est2)
# estimate the standard deviation (assuming the data is Gaussian with known mean = 0)
Est3 = mmd_est(x, model="Gaussian.scale", par1=0)
# print a summary
summary(Est3)
# test of the robustness
x[42] = 100
mean(x)
# estimate the mean and variance (assuming the data is Gaussian)
Est4 = mmd_est(x, model="Gaussian")
summary(Est4)
Run the code above in your browser using DataLab