lundi 22 juin 2020

Issues with Python scipy rv_continuous implementation

I'm trying to create a subclass of rv_continuous with a custom distribution for which I can calculate the pdf through a number of functions.

Here's what I've done so far

import numpy as np
from scipy.stats import rv_continuous

ancillary functions

def func1(xx, a_, b_, rho, m, sigma):
    return a_ + b_*(rho*(xx-m) + np.sqrt((xx-m)*(xx-m) + sigma*sigma))

def func2(xx, a_, b_, rho, m, sigma):
    sig2 = sigma*sigma
    return b_*(rho*np.sqrt((xx-m)*(xx-m)+sig2)+xx-m)/(np.sqrt((xx-m)*(xx-m)+sig2))

def func3(xx, a_, b_, rho, m, sigma):
    sig2 = sigma*sigma
    return b_*sig2/(np.sqrt((xx-m)*(xx-m)+sig2)*((xx-m)*(xx-m)+sig2))

def func4(xx, a_, b_, rho, m, sigma):
    w = func1(xx, a_, b_, rho, m, sigma)
    w1 = func2(xx, a_, b_, rho, m, sigma)
    w2 = func3(xx, a_, b_, rho, m, sigma)
    return (1.-0.5*xx*w1/w)*(1.0-0.5*xx*w1/w) - 0.25*w1*w1*(0.25 + 1./w) + 0.5*w2

def func5(xx, a_, b_, rho, m, sigma):
    vsqrt = np.sqrt(func1(xx, a_, b_, rho, m, sigma))
    return -xx/vsqrt - 0.5*vsqrt

density function eventually

def density(xx, a_, b_, rho, m, sigma):
    dm = func5(xx, a_, b_, rho, m, sigma)
    return func4(xx, a_, b_, rho, m, sigma)*np.exp(-0.5*dm*dm)/np.sqrt(2.*np.pi*func1(xx, a_, b_, rho, m, sigma))

a set of parameters

Params = 1.0073, 0.3401026, -0.8, 0.000830, 0.5109564

check the pdf from function

xmin, xmax, nbPoints = -10., 10., 2000
x_real = np.linspace(xmin, xmax, nbPoints)

den_from_func = density(x_real, *Params)

now construct my distribution class

class density_gen(rv_continuous):
    def _pdf(self, x, a_hat, b_hat, rho, m, sigma):
        return density(x, a_hat, b_hat, rho, m, sigma)

instantiation

my_density = density_gen(name='density_gen')

my_density.a, my_density.b, my_density.numargs

As I've specified _pdf I should have a working distribution instance

this works

pdf = my_density._pdf(x_real, *Params)

cdf works too albeit it's extremely slow

cdf = my_density._cdf(x_real, *Params)
my_density._cdf(0.1, *Params)

but for all the other methods I get nans, for instance

my_density.mean(*Params)    
my_density.ppf(0.01, *Params)

What I'm doing wrong here?




Aucun commentaire:

Enregistrer un commentaire