calculs_skku_exp Subroutine

private subroutine calculs_skku_exp(bounds, lg, ssk, sku)

Note

Function to calculate the skewness and kurtosis of an exponential series.
The principle is the same as calculs_skku_tan, however it fits better some particular series quite binary (roughly two heights).

Arguments

Type IntentOptional Attributes Name
real(kind=R8), intent(in), dimension(1:2) :: bounds

interval limits [-(1/bounds(1)-1), +(1/bounds(2)-1)]
integer(kind=I4), intent(in) :: lg

vec size
real(kind=R8), intent(out) :: ssk

theoretical Ssk
real(kind=R8), intent(out) :: sku

theoretical Sku

Calls

proc~~calculs_skku_exp~~CallsGraph proc~calculs_skku_exp calculs_skku_exp proc~add_expo add_expo proc~calculs_skku_exp->proc~add_expo proc~expo expo proc~add_expo->proc~expo

Called by

proc~~calculs_skku_exp~~CalledByGraph proc~calculs_skku_exp calculs_skku_exp proc~calculs_skku_generique calculs_skku_generique proc~calculs_skku_generique->proc~calculs_skku_exp proc~fitness_skku_anal fitness_skku_anal proc~fitness_skku_anal->proc~calculs_skku_exp proc~cost_func_skku cost_func_skku proc~cost_func_skku->proc~fitness_skku_anal

Source Code

   subroutine calculs_skku_exp(bounds, lg, ssk, sku)
   !================================================================================================
   !<@note Function to calculate the skewness and kurtosis of an **exponential** series.<br/>
   !< The principle is the same as [[calculs_skku_tan]], however it fits better some particular
   !< series quite binary (roughly two heights).
   !<
   !<@endnote
   !------------------------------------------------------------------------------------------------
   implicit none
   real   (kind=R8), intent(in), dimension(1:2) :: bounds  !! *interval limits* [-(1/bounds(1)-1), +(1/bounds(2)-1)]
   integer(kind=I4), intent(in)                 :: lg      !! *vec size*
   real   (kind=R8), intent(out)                :: ssk     !! *theoretical Ssk*
   real   (kind=R8), intent(out)                :: sku     !! *theoretical Sku*

      real(kind=R8) :: xa, xb, xk, mu, si, sk, ku, a, b, pente
      real(kind=R8) :: h, hh, b1, b2, alp, bet
      real(kind=R8) :: exp1b, exp1a, exp2b, exp2a, exp3b, exp3a, exp4b, exp4a, tmp1a, tmp1b, tmp2a, tmp2b, tmp3a, tmp3b, tmp4a, tmp4b
      integer(kind=I4) :: deb, fin

      pente = UN

      !---------------------------------------------
      ! WXMaxima file
      !---------------------------------------------
      !kill(all);
      !load (f90)$
      !
      !f11(x):=-1+exp(+xk*x)$
      !g11(x):=+1-exp(-xk*x)$
      !
      !f21(x):=f11(x)-mu$
      !g21(x):=g11(x)-mu$
      !
      !f31(x):=f21(x)/si$
      !g31(x):=g21(x)/si$
      !
      !assume(v>0.)$
      !assume(u>0.)$
      !/* [wxMaxima: input   end   ] */
      !
      !/* [wxMaxima: input   start ] */
      !I11:integrate(f11(x),x,-u,0)$
      !J11:integrate(g11(x),x,+0,v)$
      !I11:subst(1./xa-1.,u,I11)$
      !J11:subst(1./xb-1.,v,J11)$
      !I11:I11+J11$
      !I11:expand(trigsimp(I11))$
      !f90(I11);
      !/* [wxMaxima: input   end   ] */
      !
      !/* [wxMaxima: input   start ] */
      !I21:integrate(f21(x)^2,x,-u,0)$
      !J21:integrate(g21(x)^2,x,+0,v)$
      !I21:subst(1./xa-1.,u,I21)$
      !J21:subst(1./xb-1.,v,J21)$
      !I21:I21+J21$
      !I21:expand(trigsimp(I21))$
      !f90(I21);
      !/* [wxMaxima: input   end   ] */
      !
      !/* [wxMaxima: input   start ] */
      !I31:integrate(f31(x)^3,x,-u,0)$
      !J31:integrate(g31(x)^3,x,+0,v)$
      !I31:subst(1./xa-1.,u,I31)$
      !J31:subst(1./xb-1.,v,J31)$
      !I31:I31+J31$
      !I31:expand(trigsimp(I31))$
      !f90(I31);
      !/* [wxMaxima: input   end   ] */
      !
      !/* [wxMaxima: input   start ] */
      !I41:integrate(f31(x)^4,x,-u,0)$
      !J41:integrate(g31(x)^4,x,+0,v)$
      !I41:subst(1./xa-1.,u,I41)$
      !J41:subst(1./xb-1.,v,J41)$
      !I41:I41+J41$
      !I41:expand(trigsimp(I41))$
      !f90(I41);
      !---------------------------------------------

      a = max(bounds(1), 1.e-6_R8)
      b = max(bounds(2), 1.e-6_R8)

      hh = (-2._R8+UN/a+UN/b)/(lg-1)
      h  = hh

      xa = a
      xb = b
      xk = pente

      b1 = -(UN-a)/a
      b2 = +(UN-b)/b

      alp = -(b2-lg*b1)/(b2-b1)
      bet =      (lg-1)/(b2-b1)

      deb = 1
      fin = lg

      tmp1a = 1*(xk-  xk/xa) ; tmp1a = max(-0.9*HIG_E8, tmp1a) ; tmp1a = min(+0.9*HIG_E8, tmp1a)
      tmp1b = 1*(xk-  xk/xb) ; tmp1b = max(-0.9*HIG_E8, tmp1b) ; tmp1b = min(+0.9*HIG_E8, tmp1b)

      tmp2a = 2*(xk-  xk/xa) ; tmp2a = max(-0.9*HIG_E8, tmp2a) ; tmp2a = min(+0.9*HIG_E8, tmp2a)
      tmp2b = 2*(xk-  xk/xb) ; tmp2b = max(-0.9*HIG_E8, tmp2b) ; tmp2b = min(+0.9*HIG_E8, tmp2b)

      tmp3a = 3*(xk-  xk/xa) ; tmp3a = max(-0.9*HIG_E8, tmp3a) ; tmp3a = min(+0.9*HIG_E8, tmp3a)
      tmp3b = 3*(xk-  xk/xb) ; tmp3b = max(-0.9*HIG_E8, tmp3b) ; tmp3b = min(+0.9*HIG_E8, tmp3b)

      tmp4a = 4*(xk-  xk/xa) ; tmp4a = max(-0.9*HIG_E8, tmp4a) ; tmp4a = min(+0.9*HIG_E8, tmp4a)
      tmp4b = 4*(xk-  xk/xb) ; tmp4b = max(-0.9*HIG_E8, tmp4b) ; tmp4b = min(+0.9*HIG_E8, tmp4b)

      exp1a = exp(tmp1a)
      exp1b = exp(tmp1b)
      exp2a = exp(tmp2a)
      exp2b = exp(tmp2b)
      exp3a = exp(tmp3a)
      exp3b = exp(tmp3b)
      exp4a = exp(tmp4a)
      exp4b = exp(tmp4b)

      mu = exp1b/xk-exp1a/xk+1/xb-1/xa
      mu = (UN/h)*mu +add_expo(1, deb, fin, alp, bet, mu=0._R8, si=1._R8)
      mu = mu/lg
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      si = -2*mu*exp1b/xk+2*exp1b/xk-exp2b/xk/2.0+2*mu*exp1a/xk+2*exp1a/xk-exp2a/xk/2.0      &  !
           -3/xk+mu**2/xb-2*mu/xb+1/xb+mu**2/xa+2*mu/xa+1/xa-2*mu**2-2                          !
      si = (UN/h)*si +add_expo(2, deb, fin, alp, bet, mu, si=1._R8)
      si = si/lg
      si = sqrt(si)
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      sk = 3*mu**2*exp1b/(si**3*xk)-6*mu*exp1b/(si**3*xk)+3*exp1b/(si**3*xk)+3.0*mu*exp2b/(2.0*si**3*xk)          &  !
           +(-3.0)*exp2b/(2.0*si**3*xk)+exp3b/(si**3*xk)/3.0-3*mu**2*exp1a/(si**3*xk)-6*mu*exp1a/(si**3*xk)       &  !
           -3*exp1a/(si**3*xk)+3.0*mu*exp2a/(2.0*si**3*xk)+3.0*exp2a/(2.0*si**3*xk)-exp3a/(si**3*xk)/3.0          &  !
           +9*mu/(si**3*xk)-mu**3/(si**3*xb)+3*mu**2/(si**3*xb)-3*mu/(si**3*xb)+1/(si**3*xb)-mu**3/(si**3*xa)     &  !
           -3*mu**2/(si**3*xa)-3*mu/(si**3*xa)-1/(si**3*xa)+2*mu**3/si**3+6*mu/si**3                                 !
      sk = (UN/h)*sk +add_expo(3, deb, fin, alp, bet, mu, si)
      sk = sk/lg
      !. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
      ku = -4*mu**3*exp1b/(si**4*xk)+12*mu**2*exp1b/(si**4*xk)-12*mu*exp1b/(si**4*xk)+4*exp1b/(si**4*xk)          &  !
           -3*mu**2*exp2b/(si**4*xk)+6*mu*exp2b/(si**4*xk)-3*exp2b/(si**4*xk)+(-4.0)*mu*exp3b/(3.0*si**4*xk)      &  !
           +4.0*exp3b/(3.0*si**4*xk)-exp4b/(si**4*xk)/4.0+4*mu**3*exp1a/(si**4*xk)+12*mu**2*exp1a/(si**4*xk)      &  !
           +12*mu*exp1a/(si**4*xk)+4*exp1a/(si**4*xk)-3*mu**2*exp2a/(si**4*xk)-6*mu*exp2a/(si**4*xk)              &  !
           -3*exp2a/(si**4*xk)+4.0*mu*exp3a/(3.0*si**4*xk)+4.0*exp3a/(3.0*si**4*xk)-exp4a/(si**4*xk)/4.0          &  !
           -18*mu**2/(si**4*xk)+(-25.0)/(6.0*si**4*xk)+mu**4/(si**4*xb)-4*mu**3/(si**4*xb)+6*mu**2/(si**4*xb)     &  !
           -4*mu/(si**4*xb)+1/(si**4*xb)+mu**4/(si**4*xa)+4*mu**3/(si**4*xa)+6*mu**2/(si**4*xa)                   &  !
           +4*mu/(si**4*xa)+1/(si**4*xa)-2*mu**4/si**4-12*mu**2/si**4-2/si**4                                        !
      ku = (UN/h)*ku +add_expo(4, deb, fin, alp, bet, mu, si)
      ku = ku/lg

      ssk = sk
      sku = ku

   return
   endsubroutine calculs_skku_exp