About

Implementation of the IA2RMS algorithm for univariate densities defined for real values.

In this implementation, the sequence of proposal densities is composed of two exponential tails and uniform or linear non-overlapping piecewise densities in between.

For more information, see ourICASSP conference paper and IEEE Transactions and Signal Processing Paper. This code is maintained by the authors Luca Martino, Jesse Read, and David Luengo.

Important note: This is a second preliminary MATLAB version, not completely optimized, and thus running times should not be still compared directly to inbuilt MATLAB functions (that have been partially written in Java or C and are thus faster).

Download

Download the latest release of IA2RMS here

Or checkout the most recent code directly from the repository

User Parameters

Main Fns:

  • IA2RMS.m
Input:
  • f
    target distribution (normalizing constant not necessary) -- user-supplied parameter
  • S
    vector of initial support points
  • M
    number of desired iterations on the chain
  • type
    type of proposal construction (0/1) = (order zero, constant pieces / order one, linear pieces). In both cases there is exponential decay in the tails.
Output:
  • x
    generated samples

Example Usage

Example 1 (generic pdf often used in Matlab tests):

>> f = @(x) exp( -x.^2/2).*(1+(sin(3*x)).^2).* (1+(cos(5*x).^2));
>> S=[-1 0 1];
>> M=1000; % number of samples
>> type=0; % type of proposal construction (order 0)
>> x=A2RMS(f,S,M,type);

Example 2 (Laplacian pdf [two tails]):

>> f = @(x) exp(-abs(x));
>> S=[-1 0 1];
>> x=A2RMS(f,S,1000,0);

Example 3 (Exponential pdf [one tail]):

>> f = @(x) exp(-x) .* (x > 0);
>> S=[-1 0 1];
>> x=A2RMS(f,S,1000,0);

Example 4 (mixture of Gaussians):

function f = my_target(x) 

    w=[0.3 0.3 0.4];     % weights 
    mu=[-5 1 7];         % means 
    v=[1 1 1];           % variances
    for i=1:length(x) 
        f(i)=sum(w.*(1./sqrt(2*pi*v)).*exp(-(x(i)-mu).^2./(2.*v))); 
    end
    

>> f = @(x) my_target(x);
>> S=[-1 0 1];
>> x=A2RMS(f,S,1000,0);

Example 5 (Lévy pdf [heavy tail]):

>> f = @(x) 1./(x.^(3/2)) .* exp(-1./x) .* (x > 0);
>> S=[-1 0.1 1];
>> x=A2RMS(f,S,1000,0);
UPM UPM UPM UPM UC3M