get_moments Subroutine

public subroutine get_moments(data, ave, sdev, var, skew, curt)

Arguments

Type IntentOptional Attributes Name
real(kind=8), intent(in), dimension(:) :: data
real(kind=8), intent(out), optional :: ave
real(kind=8), intent(out), optional :: sdev
real(kind=8), intent(out), optional :: var
real(kind=8), intent(out), optional :: skew
real(kind=8), intent(out), optional :: curt

Called by

proc~~get_moments~~CalledByGraph proc~get_moments get_moments proc~get_curt get_curt proc~get_curt->proc~get_moments proc~get_mean get_mean proc~get_mean->proc~get_moments proc~get_sd get_sd proc~get_sd->proc~get_moments proc~get_skew get_skew proc~get_skew->proc~get_moments proc~get_var get_var proc~get_var->proc~get_moments

Source Code

  subroutine get_moments(data,ave,sdev,var,skew,curt)
    real(8), intent(out),optional     :: ave,sdev,var,skew,curt
    real(8)                           :: ave_,sdev_,var_,skew_,curt_
    real(8), dimension(:), intent(in) :: data
    integer                           :: n
    real(8)                           :: ep,adev
    real(8), dimension(size(data))    :: p,s
    n=size(data)
    if (n <= 1) then
       print*,'data.size must be at least 2'
       stop
    endif
    ave_=sum(data(:))/n
    s(:)=data(:)-ave_
    ep=sum(s(:))
    adev=sum(abs(s(:)))/n
    p(:)=s(:)*s(:)
    var_=sum(p(:))
    p(:)=p(:)*s(:)
    skew_=sum(p(:))
    p(:)=p(:)*s(:)
    curt_=sum(p(:))
    var_=(var_-ep**2/n)/(n-1)
    sdev_=sqrt(var_)
    if (var_ /= 0.0) then
       skew_=skew_/(n*sdev_**3)
       curt_=curt_/(n*var_**2)-3.0d0
    else
       skew_=0.d0
       curt_=0.d0
    endif
    if(present(ave))ave=ave_
    if(present(sdev))sdev=sdev_
    if(present(var))var=var_
    if(present(skew))skew=skew_
    if(present(curt))curt=curt_
  end subroutine get_moments