# Various .pdf and .cdf's and something more...
#
import numpy as np
import math
from scipy.stats import norm, poisson

def average(d):
    s = 0
    n = 0
    for x in d:
        s += x
        n += 1
    return s/n

def sd(d):
    s = 0
    n = 0
    mean = average(d)
    for x in d:
        s += (x-mean)**2
        n += 1
    s *= 1/(n-1)        
    return math.sqrt(s)



def Phi(x):
    return norm.cdf(x,0,1)
    
def phi(x):
    return norm.pdf(x,0,1)  
    
def cdfNormal(x, mu, sigma):
    return Phi((x-mu)/sigma) 
    
def pdfNormal(x, mu, sigma):
    return    (0.3989422804014327/sigma)*math.exp(-(x-mu)**2/(2*sigma*sigma))


def B(r, mu , sigma): # Backorder function
    return sigma * phi((r - mu) / sigma) + (mu - r) * (1 - Phi((r - mu) / sigma))

    
def cdfWeibul(x,alpha,lmbda):
    return 1 - math.exp(-((x*lmbda)**alpha))

def RWeibul(x,alpha,lmbda):
    return math.exp(-((x*lmbda)**alpha))
    
def pdfWeibul(x,alpha,lmbda):
    return lmbda*alpha*(lmbda*x)**(alpha-1)* math.exp(-((x*lmbda)**alpha))    
    
def cdfInverseGauss(x,alpha,beta):
    if x <= 0:
        return 0
    hlp = Phi(-math.sqrt(x*beta)/alpha-math.sqrt(beta/x))
    if hlp > 0: # Increase precision when multiplying with
        hlp = math.exp(math.log(hlp)+2*beta/alpha)
    return Phi(math.sqrt(x*beta)/alpha-math.sqrt(beta/x)) + hlp
    
def pdfInverseGauss(x,alpha,beta):
    if x <= 0:
        return 0
    return math.sqrt(beta/(2*math.pi*x*x*x))*math.exp(-beta*(x-alpha)*(x-alpha)/(2*alpha*alpha*x))

def cdfFirstHitW(x,mu,sigma,l):
    return cdfInverseGauss(x,l/mu,(l*l/(sigma*sigma)))

def pdfPoisson(x,lmbda):
    return poisson.pmf(x, lmbda)


def cdfPoisson(x,lmbda):
    return poisson.cdf(x, lmbda)