In [5]:

    

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, lognorm
from scipy.interpolate import UnivariateSpline


    

In [78]:

    

x = np.linspace(0, 100000, 1000)
dist = [
    (0.75, lognorm(s=1, scale=10000)),
    (0.1, norm(loc=50000, scale=5000)),
    (0.05, norm(loc=10000, scale=1500)),
    (0.05, norm(loc=30000, scale=5000)),
    (0.01, norm(loc=80000, scale=2000)),
    (0.005, norm(loc=90000, scale=500))
]
Fx = sum([p * d.cdf(x) for p, d in dist]) / sum(p for p, d in dist)
fx = sum([p * d.pdf(x) for p, d in dist]) / sum(p for p, d in dist)



Fx = bra.cdf(x)
fx = bra.pdf(x)
dfx = np.pad(fx[1:]-fx[:-1],[1,0],'constant',constant_values=float("NaN"))
ddfx = np.pad(dfx[1:]-dfx[:-1],[0,1],'constant',constant_values=float("NaN"))
abs_ddfx = np.abs(ddfx)


    

In [121]:

    

##########################################
plt.rcParams["figure.figsize"] = (12,15)
fig, ax = plt.subplots(3, 1)
##########################################
K = 100

smoothing = 1
direction = None
adjustment = 2
for iterations in range(100): # maximum iterations
    spline_fx = UnivariateSpline(x, fx, k=1, s=smoothing)
    knots = spline_fx.get_knots()
    if (len(knots)==K):
        break
    elif (len(knots) > K):
        if (direction == -1):
            adjustment /= 2
        smoothing *= adjustment
        direction = 1
    elif (len(knots) < K):
        if (direction == 1):
            adjustment /= 2
        smoothing /= adjustment
        direction = -1

ax[0].plot(x,Fx)
ax[0].plot(knots,spline_fx.antiderivative()(knots),color="r")
ax[0].scatter(knots,spline_fx.antiderivative()(knots),color="r")

ax[1].plot(x,fx)
ax[1].plot(x,spline_fx(x),color="r")
ax[1].scatter(knots,spline_fx(knots),color="r")

mean = sum(((x[1:] - x[:-1])*(x[1:] + x[:-1])/2*spline_fx((x[1:] + x[:-1])/2)))
ax[2].plot(spline_fx.antiderivative()(x[1:]), np.cumsum(((x[1:] - x[:-1])*(x[1:] + x[:-1])/2*spline_fx((x[1:] + x[:-1])/2)))/mean)
ax[2].plot(spline_fx.antiderivative()(knots[1:]), np.cumsum(((knots[1:] - knots[:-1])*(knots[1:] + knots[:-1])/2*spline_fx((knots[1:] + knots[:-1])/2)))/mean, "r")

#ax[2].plot(spline_fx.antiderivative()(knots), np.cumsum(knots)/mean,"r")


    

    
Out[121]:





[<matplotlib.lines.Line2D at 0x11c2ae400>]






    

In [114]:

    

knots


    

    
Out[114]:





array([      0.        ,     100.1001001 ,     200.2002002 ,
           300.3003003 ,     400.4004004 ,     500.5005005 ,
           600.6006006 ,     800.8008008 ,    1001.001001  ,
          1201.2012012 ,    1401.4014014 ,    1601.6016016 ,
          1801.8018018 ,    2002.002002  ,    2202.2022022 ,
          2402.4024024 ,    2602.6026026 ,    2802.8028028 ,
          3003.003003  ,    3203.2032032 ,    3603.6036036 ,
          4004.004004  ,    4404.4044044 ,    4804.8048048 ,
          6306.30630631,    6706.70670671,    7107.10710711,
          7507.50750751,    7907.90790791,    8308.30830831,
          8708.70870871,    9109.10910911,    9409.40940941,
          9809.80980981,   10210.21021021,   10610.61061061,
         11011.01101101,   11811.81181181,   12212.21221221,
         12512.51251251,   12912.91291291,   13313.31331331,
         13713.71371371,   14114.11411411,   14914.91491491,
         15715.71571572,   16516.51651652,   17317.31731732,
         18818.81881882,   20420.42042042,   21921.92192192,
         23523.52352352,   25025.02502503,   28228.22822823,
         29829.82982983,   31331.33133133,   32932.93293293,
         34434.43443443,   36036.03603604,   37537.53753754,
         39139.13913914,   40740.74074074,   42342.34234234,
         43843.84384384,   46946.94694695,   47747.74774775,
         48548.54854855,   49349.34934935,   50050.05005005,
         50850.85085085,   51651.65165165,   53253.25325325,
         56356.35635636,   57957.95795796,   59459.45945946,
         61061.06106106,   62562.56256256,   64164.16416416,
         65765.76576577,   68868.86886887,   71971.97197197,
         75075.07507508,   76676.67667668,   78278.27827828,
         79079.07907908,   79879.87987988,   80680.68068068,
         81381.38138138,   82182.18218218,   82982.98298298,
         84484.48448448,   86086.08608609,   87587.58758759,
         88388.38838839,   88788.78878879,   88988.98898899,
         89189.18918919,   89589.58958959,   89789.78978979,
         89989.98998999,   90190.19019019,   90390.39039039,
         90690.69069069,   90890.89089089,   91091.09109109,
         91491.49149149,   92292.29229229,   93793.79379379,  100000.        ])

In [119]:

    

print(spline_fx(knots).tolist())


    

    




[-7.00170692036934e-08, 1.0947698984041765e-07, 7.875501460248058e-07, 2.233225641141353e-06, 4.367748843994068e-06, 6.982827630431993e-06, 9.883968818968741e-06, 1.5993319418932447e-05, 2.1942528496166424e-05, 2.7364251973301945e-05, 3.2128566113229755e-05, 3.6216509175664014e-05, 3.966022097706759e-05, 4.251482104287394e-05, 4.484330961707371e-05, 4.6708800806492374e-05, 4.8170734661740725e-05, 4.9282533112208486e-05, 5.009365098256787e-05, 5.067629256107429e-05, 5.116587136598674e-05, 5.0977455222565017e-05, 5.033101037603566e-05, 4.940432804104205e-05, 4.4841847526280174e-05, 4.390070838596992e-05, 4.3267392246011555e-05, 4.306643318209113e-05, 4.332385382364866e-05, 4.3965082540462756e-05, 4.479142209702402e-05, 4.546798330093391e-05, 4.568981883260911e-05, 4.5345383510587105e-05, 4.4078525351016345e-05, 4.191109313518956e-05, 3.904687471932377e-05, 3.252140510902292e-05, 2.9513644896766766e-05, 2.753462220785759e-05, 2.530556930014244e-05, 2.3542659288529663e-05, 2.215642320691538e-05, 2.1011777367467665e-05, 1.9268586562917076e-05, 1.7879275582799446e-05, 1.6654970376875214e-05, 1.5539696647190365e-05, 1.3791048963842694e-05, 1.239297204631472e-05, 1.1484348391330815e-05, 1.0901119465855798e-05, 1.0654664030267313e-05, 1.0324907975428293e-05, 9.87725690344552e-06, 9.176337489298979e-06, 8.144915170642634e-06, 7.041514461643501e-06, 5.926420376977343e-06, 5.108488237536243e-06, 4.6378114381702e-06, 4.691422613398233e-06, 5.29975665888446e-06, 6.322827783867369e-06, 8.903601448247196e-06, 9.397491962230248e-06, 9.781534764117548e-06, 9.971335828620402e-06, 9.979494468344458e-06, 9.785421885967264e-06, 9.4282676909164e-06, 8.169227217014803e-06, 4.902913993394706e-06, 3.4422624814364983e-06, 2.4166576826473992e-06, 1.680272923194885e-06, 1.2598415127817301e-06, 9.999554742229113e-07, 8.470307578838619e-07, 7.029123589058478e-07, 6.060482084282032e-07, 5.992898877586231e-07, 9.86553056650152e-07, 1.8955671012470675e-06, 2.3440040388912675e-06, 2.5383376168642756e-06, 2.4066299507543764e-06, 2.0581729628358483e-06, 1.5536291746176355e-06, 1.0447784961565498e-06, 5.08037767359809e-07, 3.666041477246561e-07, 3.3687368072879814e-07, 3.320362668250251e-07, 5.162154736703169e-07, 8.405150275396369e-07, 1.3875312798730582e-06, 3.249784697355953e-06, 4.1231092405514735e-06, 4.47750237125129e-06, 4.180521324243001e-06, 3.359390077813789e-06, 1.8795783087788202e-06, 1.1286108548541142e-06, 6.425614522424119e-07, 3.097494475974095e-07, 2.829659551176327e-07, 2.6915233247840224e-07, 2.1763938469149517e-07]

In [ ]: