Type | Intent | Optional | 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 |
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