In [3]:

    

from matplotlib.pyplot import plot


    

In [4]:

    

from numpy import pi, exp, tan, arange, vectorize, cos, sqrt, sin, abs


    

In [5]:

    

from numpy import power as pow


    

In [6]:

    

from scipy import interpolate


    

In [7]:

    

def hgauss_simple(x, event, var):
    """
    I keep here the original idea with the simple implementation, so that if i want to /simplify/ the formula using trig
    equivalences, optimize it or experiment with it, i'll do it in another place, retaining the original idea here as a reference.
    
    This is an experimentation with event modeling. In a deadline [0, dl] one can have multiple events, so that Events(t) has
    sum(Ni*dirac(t-ti)), whrere ti is the time when some Ni events happened.
    
    One can then interpret a single event dirac as a probability density with variance->inf, and loosen up that requirement on the
    Events(t) to get the total prob density, keeping in mind to normalize the area by sum(Ni), total number of single events foreach
    event prob distribution
    
    One problem with this is that the tails of the single distributions are not confined in the deadline, so i decided to project
    the distribution on an hyperplane. Doing so the area of the distr changes, but one thing i learned is that i dont have to
    normalize by the integral of the new distr, i just scale back the transformation using its jacobian preserving the original
    area.
    """
    J = pi/(dl*pow(sin(pi*x/dl), 2))
    xh = tan((x-dl/2.0)*pi/dl) + dl/2.0
    # mul gauss with the jacobian of the transformation, so that it's area stays the same
    return var*J*exp(-pow((xh-event)*var, 2))


    

In [8]:

    

dl, sigma, dt = 5.0, 0.1, 0.01


    

In [9]:

    

x = arange(0, dl+dt, dt)


    

In [10]:

    

def hgauss_fst(x, event, cvar):
    # first attempt to normalize variation coeff
    var = (1+2*pow(sin(pi*event/dl), 2))*cvar
    J = var*pi/(dl*pow(sin(pi*x/dl), 2))
    xh = -1.0/tan(x*pi/dl) + dl/2.0
    return J*exp(-pow((xh-event)*var, 2))


    

In [11]:

    

def hgauss_snd(x, event, var):
    # second attempt, useless but interesting
    inverse_J = var*dl*pow(sin(pi*x/dl), 2)/pi
    xh = tan((x-dl/2)*pi/dl) + dl/2
    return var*exp(-pow((xh-event)*inverse_J, 2))


    

In [12]:

    

def hgauss_norm(x, event, cvar):
    # last attempt to normalize variation coeff
    if x == 0 or x == dl or event == cvar: return 0.0
    var = dl/abs(1.0/tan(pi*(event-cvar)/dl)-1.0/tan(pi*(event+cvar)/dl))
    J = var*pi/(dl*pow(sin(pi*x/dl), 2))
    xh = -1.0/tan(x*pi/dl) + dl/2.0
    return J*exp(-pow((xh-event)*var, 2))


    

In [13]:

    

f = vectorize(hgauss_norm)


    

In [14]:

    

for event in arange(0, dl+0.5, 0.5):
    plot(x, f(x, event, 0.3))