gamma

518 days ago by jason3

class BetaLP(): def __init__(self,mu,sigma,alpha1,beta1,lambda1,c1,alpha2,beta2,lambda2,c2): for v in ('mu','sigma','alpha1','beta1','lambda1','c1','alpha2','beta2','lambda2','c2'): exec('self.'+v+'='+v) if (not min(self.alpha1,self.alpha2,self.beta1,self.beta2)>0 or not min(self.c1,self.c2)>=0 or not min(self.lambda1,self.lambda2)>0 or not max(self.lambda1,self.lambda2)<3): raise ValueError(HTMLPrint('in BetaLP.init(): the parameters do not '+ 'satisfy the constraints '+ 'alpha_i>0, beta_i>0, c_i>=0, lambda_i \in (0,3)','error')) if self.lambda1==1 or self.lambda2==1 or self.lambda1==2 or self.lambda2==2: raise ValueError(HTMLPrint('in BetaLP.init(): paramater case '+ 'lambda_i \in {1,2} not implemented','error')) self.descriptiveName="Kuznetsov beta family of LP's" y=var('y') self.B=(gamma(x)*gamma(y)/gamma(x+y)).function(x,y) self.psi=(diff(gamma(x),x)/gamma(x)).function(x) self.ggamma=self.c1*self.B(self.alpha1,1-self.lambda1)/self.beta1+self.c2*self.B(self.alpha2,1-self.lambda2)/self.beta2 self.rho=(self.c1*self.B(self.alpha1,1-self.lambda1)*(self.psi(1+self.alpha1-self.lambda1)-self.psi(self.alpha1))/self.beta1^2 -self.c2*self.B(self.alpha2,1-self.lambda2)*(self.psi(1+self.alpha2- self.lambda2)-self.psi(self.alpha2))/self.beta2^2 -self.mu) def getCharacteristicExponent(self): i=sqrt(-1) ###just in case it was redefined somewhere else expr=self.sigma^2*x^2/2+i*self.rho*x-self.c1*self.B(self.alpha1-i*x/self.beta1,1-self.lambda1)/self.beta1 expr+=-self.c2*self.B(self.alpha2+i*x/self.beta2,1-self.lambda2)/self.beta2+self.ggamma return expr.function(x) def getkthMomentAt(self,k,t): charFunc=self.getCharacteristicExponent() i=sqrt(-1) f=(exp(-t*charFunc(-i*x))).function(x) for _ in range(k): f=f.derivative() #ff=symbolic_expression(str(f(x))).function(x) ###temp. workaround return f 
       
BetaLP1=BetaLP(mu=0,sigma=3/10,alpha1=1,beta1=1,lambda1=5/2,c1=1,alpha2=1,beta2=1,lambda2=5/2,c2=1) f=BetaLP1.getkthMomentAt(1,1) 
       
mygammaexpr=f.operands()[0].operands()[0].operands()[1].operands()[0] 
       
mygammaexpr(x=0) 
       
Traceback (click to the left of this block for traceback)
...
ValueError: argument must be >= -1
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "_sage_input_6.py", line 10, in <module>
    exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("bXlnYW1tYWV4cHIoeD0wKQ=="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))
  File "", line 1, in <module>
    
  File "/tmp/tmpiz4apw/___code___.py", line 3, in <module>
    exec compile(u'mygammaexpr(x=_sage_const_0 )
  File "", line 1, in <module>
    
  File "expression.pyx", line 3473, in sage.symbolic.expression.Expression.__call__ (sage/symbolic/expression.cpp:15647)
  File "ring.pyx", line 640, in sage.symbolic.ring.SymbolicRing._call_element_ (sage/symbolic/ring.cpp:6537)
  File "expression.pyx", line 3324, in sage.symbolic.expression.Expression.substitute (sage/symbolic/expression.cpp:15026)
  File "pynac.pyx", line 1047, in sage.symbolic.pynac.py_doublefactorial (sage/symbolic/pynac.cpp:9456)
ValueError: argument must be >= -1
gamma(-x-1/2)-mygammaexpr 
       
0
0
mygammaexpr(x=0) 
       
-2*sqrt(pi)
-2*sqrt(pi)