Why I can't build TM from non standard normal distribution to Gumbel distribution ?

Question

Hi,

Thank you for your software and all the support provided.

I am trying to use your TransportMaps library to solve some of my examples but I encounter some difficulties. I went back to a simple example which is not working.

I am trying to build a map between a gaussian distribution and a gumbel distribution (mu=-63, beta=3) (which is close to a posterior distribution I want to sample). When I use a standard normal gaussian distribution (rho) as reference distribution it’s working well. But when I use another gaussian distribution (mean=-60, stdd=8) (rho2) I can’t reach a good approximation building the map whatever the order of the map I use (see file attached). For me it will be easier to build the map between closer distributions. I can’t understand what I am doing wrong. 

I would very appreciate if you can find some time to answer me.

Best regards,

Paul.

Here the python script I use:

import logging
import warnings
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats
import SpectralToolbox.Spectral1D as S1D
import TransportMaps as TM
import TransportMaps.Functionals as FUNC
import TransportMaps.Maps as MAPS
import scipy.io
import scipy.linalg as scila
import SpectralToolbox.SpectralND as SND

warnings.simplefilter("ignore")
TM.setLogLevel(logging.INFO)

import TransportMaps.Distributions as DIST

    
class GumbelDistribution(DIST.Distribution):
    def __init__(self, mu, beta):
        super(GumbelDistribution,self).__init__(1)
        self.mu = mu
        self.beta = beta
        self.dist = stats.gumbel_r(loc=mu, scale=beta)
    def pdf(self, x, params=None):
        return self.dist.pdf(x).flatten()
    def log_pdf(self, x, params=None):
        return self.dist.logpdf(x).flatten()
    def grad_x_log_pdf(self, x, params=None):
        m = self.mu
        b = self.beta
        z = (x-m)/b
        return (np.exp(-z)-1.)/b
    def hess_x_log_pdf(self, x, params=None):
        m = self.mu
        b = self.beta
        z = (x-m)/b
        return (-np.exp(-z)/b**2.)[:,:,np.newaxis]

mu = -63.
beta = 3.
pi = GumbelDistribution(mu,beta)

M=np.array([[-50],[-70],[-60]])

xnew=np.linspace(-70,-50,201)[:,np.newaxis]

mu = np.array([-60])
sigma = np.array([[8]])

rho=DIST.StandardNormalDistribution(1)
rho2=DIST.GaussianDistribution(mu,sigma)

x = np.linspace(-70., -50., 100).reshape((100,1))

order=7
T = TM.Default_IsotropicIntegratedSquaredTriangularTransportMap(
    1, order, 'full')
T2 = TM.Default_IsotropicIntegratedSquaredTriangularTransportMap(
    1, order, 'full')

push_rho=DIST.PushForwardTransportMapDistribution(T,rho)
push_rho2=DIST.PushForwardTransportMapDistribution(T2,rho2)


qtype = 3      # Gauss quadrature
qparams =[20] # Quadrature order
reg = None     # No regularization
tol = 1e-10    # Optimization tolerance
ders = 2       # Use gradient and Hessian

log = push_rho.minimize_kl_divergence(
    pi, qtype=qtype, qparams=qparams, regularization=reg,
    tol=tol, ders=ders)

log = push_rho2.minimize_kl_divergence(
    pi, qtype=qtype, qparams=qparams, regularization=reg,
    tol=tol, ders=ders)

plt.figure()
plt.plot(x, push_rho.pdf(x),'b',label=r'$push rho$');
plt.plot(x, push_rho2.pdf(x),'r',label=r'$push rho2$');
plt.plot(x, pi.pdf(x),'g',label=r'$\pi$');
plt.legend();

 

Comments

Hi, 

I am trying to use your TransportMaps library to track an airplane using ground based radar. However, I got an error : NameError: free variable 'gfk' referenced before assignment in enclosing scope.

I would very much appreciate if you can find some time to answer me.

Best regards,

My Contact Page

Answer 1

hi Paul,

The issue you're encountering is due to the basis functions that the code uses to represent maps. We use Hermite polynomials, Hermite functions, and/or radial basis functions (depending on the options you pick). When you shift the reference to have mean -60, these basis functions are not scaled and shifted properly to be particularly expressive or controlled in the region of high reference probability. In other words, our map basis is well tuned to a reference with mean roughly zero and variance roughly one.

To get around this, we recommend first using a linear map to recenter/rescale everything. A linear map is robust to translation and scaling, because that is essentially its action. Then everything should work.

best,

youssef

Answer 2

Hi Paul, Youssef's comment is spot on. Here is the code implementing his suggestion:

import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats
import TransportMaps as TM
import TransportMaps.Functionals as FUNC
import TransportMaps.Maps as MAPS

import TransportMaps.Distributions as DIST
    
class GumbelDistribution(DIST.Distribution):
    def __init__(self, mu, beta):
        super(GumbelDistribution,self).__init__(1)
        self.mu = mu
        self.beta = beta
        self.dist = stats.gumbel_r(loc=mu, scale=beta)
    def pdf(self, x, params=None):
        return self.dist.pdf(x).flatten()
    def log_pdf(self, x, params=None):
        return self.dist.logpdf(x).flatten()
    def grad_x_log_pdf(self, x, params=None):
        m = self.mu
        b = self.beta
        z = (x-m)/b
        return (np.exp(-z)-1.)/b
    def hess_x_log_pdf(self, x, params=None):
        m = self.mu
        b = self.beta
        z = (x-m)/b
        return (-np.exp(-z)/b**2.)[:,:,np.newaxis]

mu = -63.
beta = 3.
pi = GumbelDistribution(mu,beta)

M=np.array([[-50],[-70],[-60]])

xnew=np.linspace(-70,-50,201)[:,np.newaxis]

mu = np.array([-60])
sigma = np.array([[8]])

rho=DIST.StandardNormalDistribution(1)

# Linear map
L = MAPS.LinearTransportMap(mu, sigma)
# New target (pullback of pi through L)
pull_L_pi = DIST.PullBackTransportMapDistribution(L, pi)

x = np.linspace(-70., -50., 100).reshape((100,1))

order=7
T = TM.Default_IsotropicIntegratedSquaredTriangularTransportMap(
    1, order, 'full')

push_rho=DIST.PushForwardTransportMapDistribution(T,rho)


qtype = 3      # Gauss quadrature
qparams =[20] # Quadrature order
reg = None     # No regularization
tol = 1e-6    # Optimization tolerance
ders = 2       # Use gradient and Hessian

# Solve D_KL(T_\sharp \rho || L^\sharp \pi)
log = push_rho.minimize_kl_divergence(
    pull_L_pi, qtype=qtype, qparams=qparams, regularization=reg,
    tol=tol, ders=ders)

# Obtain \pi \approx L_\sharp T_\sharp \rho
push_LT_rho = DIST.PushForwardTransportMapDistribution(L, push_rho)

plt.figure()
plt.plot(x, push_LT_rho.pdf(x),'b',label=r'pushrho');
plt.plot(x, pi.pdf(x),'g',label=r'π');
plt.legend();

To know more on composition of maps, you can see this page of the tutorial.

Best,

 Daniele

Comments

Thank you very much for your both answers. Now I understand that, I can try to use it on more complex applications.

Best,

Paul