From 979902d8008eabbb643d86823d2e32e60b1f46d0 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sun, 20 Dec 2020 21:12:24 -0500 Subject: [PATCH 01/49] initial commit --- CMakeLists.txt | 108 ++++---- Makefile.manual | 108 ++++++-- stdlib_stats_distribution_normal.fypp | 366 ++++++++++++++++++++++++++ 3 files changed, 520 insertions(+), 62 deletions(-) create mode 100644 stdlib_stats_distribution_normal.fypp diff --git a/CMakeLists.txt b/CMakeLists.txt index 6143b5ab5..33144ad07 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,52 +1,70 @@ -cmake_minimum_required(VERSION 3.14.0) -project(stdlib Fortran) -enable_testing() - -include(${CMAKE_SOURCE_DIR}/cmake/stdlib.cmake) - -# --- compiler options -if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU) - add_compile_options(-fimplicit-none) - add_compile_options(-ffree-line-length-132) - add_compile_options(-Wall) - add_compile_options(-Wextra) - add_compile_options(-Wimplicit-procedure) - add_compile_options(-Wconversion-extra) - # -pedantic-errors triggers a false positive for optional arguments of elemental functions, - # see test_optval and https://gcc.gnu.org/bugzilla/show_bug.cgi?id=95446 - add_compile_options(-pedantic-errors) - if(CMAKE_Fortran_COMPILER_VERSION VERSION_GREATER_EQUAL 8.0) - add_compile_options(-std=f2018) - else() - add_compile_options(-std=f2008ts) - endif() -elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Intel) - add_compile_options(-warn declarations,general,usage,interfaces,unused) - add_compile_options(-standard-semantics) - if(CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 18.0) - add_compile_options(-stand f15) - else() - add_compile_options(-stand f18) - endif() -elseif(CMAKE_Fortran_COMPILER_ID STREQUAL PGI) - add_compile_options(-Mdclchk) -endif() +### Pre-process: .fpp -> .f90 via Fypp + +# Create a list of the files to be preprocessed +set(fppFiles + stdlib_bitsets.fypp + stdlib_bitsets_64.fypp + stdlib_bitsets_large.fypp + stdlib_io.fypp + stdlib_linalg.fypp + stdlib_linalg_diag.fypp + stdlib_optval.fypp + stdlib_stats.fypp + stdlib_stats_corr.fypp + stdlib_stats_cov.fypp + stdlib_stats_mean.fypp + stdlib_stats_moment.fypp + stdlib_stats_var.fypp + stdlib_quadrature.fypp + stdlib_quadrature_trapz.fypp + stdlib_quadrature_simps.fypp + stdlib_stats_distribution_PRNG.fypp + stdlib_stats_distribution_uniform.fypp + stdlib_stats_distribution_normal.fypp +) -# --- compiler feature checks -include(CheckFortranSourceCompiles) -include(CheckFortranSourceRuns) -check_fortran_source_compiles("error stop i; end" f18errorstop SRC_EXT f90) -check_fortran_source_compiles("real, allocatable :: array(:, :, :, :, :, :, :, :, :, :); end" f03rank SRC_EXT f90) -check_fortran_source_runs("use, intrinsic :: iso_fortran_env, only : real128; real(real128) :: x; x = x+1; end" f03real128) +# Custom preprocessor flags if(DEFINED CMAKE_MAXIMUM_RANK) - set(CMAKE_MAXIMUM_RANK ${CMAKE_MAXIMUM_RANK}) + set(fyppFlags "-DMAXRANK=${CMAKE_MAXIMUM_RANK}") +elseif(f03rank) + set(fyppFlags) +else() + set(fyppFlags "-DVERSION90") endif() -# --- find preprocessor -find_program(FYPP fypp) -if(NOT FYPP) - message(FATAL_ERROR "Preprocessor fypp not found!") +fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) + +set(SRC + stdlib_ascii.f90 + stdlib_error.f90 + stdlib_kinds.f90 + stdlib_logger.f90 + stdlib_system.F90 + ${outFiles} +) + +add_library(fortran_stdlib ${SRC}) + +set(LIB_MOD_DIR ${CMAKE_CURRENT_BINARY_DIR}/mod_files/) +set_target_properties(fortran_stdlib PROPERTIES + Fortran_MODULE_DIRECTORY ${LIB_MOD_DIR}) +target_include_directories(fortran_stdlib PUBLIC + $ + $ +) + +if(f18errorstop) + target_sources(fortran_stdlib PRIVATE f18estop.f90) +else() + target_sources(fortran_stdlib PRIVATE f08estop.f90) endif() -add_subdirectory(src) +add_subdirectory(tests) + +install(TARGETS fortran_stdlib + RUNTIME DESTINATION bin + ARCHIVE DESTINATION lib + LIBRARY DESTINATION lib + ) +install(DIRECTORY ${LIB_MOD_DIR} DESTINATION include) diff --git a/Makefile.manual b/Makefile.manual index b8280b102..725b4aec9 100644 --- a/Makefile.manual +++ b/Makefile.manual @@ -1,24 +1,98 @@ -# Fortran stdlib Makefile +SRC = f18estop.f90 \ + stdlib_ascii.f90 \ + stdlib_bitsets.f90 \ + stdlib_bitsets_64.f90 \ + stdlib_bitsets_large.f90 \ + stdlib_error.f90 \ + stdlib_io.f90 \ + stdlib_kinds.f90 \ + stdlib_linalg.f90 \ + stdlib_linalg_diag.f90 \ + stdlib_logger.f90 \ + stdlib_optval.f90 \ + stdlib_quadrature.f90 \ + stdlib_quadrature_trapz.f90 \ + stdlib_stats.f90 \ + stdlib_stats_mean.f90 \ + stdlib_stats_moment.f90 \ + stdlib_stats_var.f90 \ + stdlib_stats_distribution_PRNG.f90\ + stdlib_stats_distribution_uniform.f90 \ + stdlib_stats_distribution_normal.f90 + +LIB = libstdlib.a -FC = gfortran -FFLAGS = -Wall -Wextra -Wimplicit-interface -fPIC -g -fcheck=all -FYPPFLAGS= -export FC -export FFLAGS -export FYPPFLAGS -.PHONY: all clean test +OBJS = $(SRC:.f90=.o) +MODS = $(OBJS:.o=.mod) +SMODS = $(OBJS:.o=*.smod) -all: - $(MAKE) -f Makefile.manual --directory=src - $(MAKE) -f Makefile.manual --directory=src/tests +.PHONY: all clean -test: - $(MAKE) -f Makefile.manual --directory=src/tests test - @echo - @echo "All tests passed." +all: $(LIB) + +$(LIB): $(OBJS) + ar rcs $@ $(OBJS) clean: - $(MAKE) -f Makefile.manual clean --directory=src - $(MAKE) -f Makefile.manual clean --directory=src/tests + $(RM) $(LIB) $(OBJS) $(MODS) $(SMODS) + +%.o: %.f90 + $(FC) $(FFLAGS) -c $< + +%.f90: %.fypp + fypp $(FYPPFLAGS) $< $@ + +# Fortran module dependencies +f18estop.o: stdlib_error.o +stdlib_bitsets.o: stdlib_kinds.o +stdlib_bitsets_64.o: stdlib_bitsets.o +stdlib_bitsets_large.o: stdlib_bitsets.o +stdlib_error.o: stdlib_optval.o +stdlib_io.o: \ + stdlib_error.o \ + stdlib_optval.o \ + stdlib_kinds.o +stdlib_linalg_diag.o: stdlib_kinds.o +stdlib_logger.o: stdlib_ascii.o stdlib_optval.o +stdlib_optval.o: stdlib_kinds.o +stdlib_quadrature.o: stdlib_kinds.o +stdlib_stats_mean.o: \ + stdlib_optval.o \ + stdlib_kinds.o \ + stdlib_stats.o +stdlib_stats_moment.o: \ + stdlib_optval.o \ + stdlib_kinds.o \ + stdlib_stats.o +stdlib_stats_var.o: \ + stdlib_optval.o \ + stdlib_kinds.o \ + stdlib_stats.o +stdlib_stats_distribution_PRNG.o: stdlib_kinds.o +stdlib_stats_distribution_uniform.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_PRNG.o +stdlib_stats_distribution_normal.o \ + stdlib_kinds.o + stdlib_error.o \ + stdlib_stats_distribution.PRNG.o \ + stdlib_stats_distribution.uniform.o + +# Fortran sources that are built from fypp templates +stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp +stdlib_bitsets_large.f90: stdlib_bitsets_large.fypp +stdlib_bitsets.f90: stdlib_bitsets.fypp +stdlib_io.f90: stdlib_io.fypp +stdlib_linalg.f90: stdlib_linalg.fypp +stdlib_linalg_diag.f90: stdlib_linalg_diag.fypp +stdlib_quadrature.f90: stdlib_quadrature.fypp +stdlib_stats.f90: stdlib_stats.fypp +stdlib_stats_mean.f90: stdlib_stats_mean.fypp +stdlib_stats_moment.f90: stdlib_stats_moment.fypp +stdlib_stats_var.f90: stdlib_stats_var.fypp +stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp +stdlib_stats_distribution_uniform.f90: stdlib_stats_distribution_uniform.fypp +stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp \ No newline at end of file diff --git a/stdlib_stats_distribution_normal.fypp b/stdlib_stats_distribution_normal.fypp new file mode 100644 index 000000000..de55412ab --- /dev/null +++ b/stdlib_stats_distribution_normal.fypp @@ -0,0 +1,366 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +Module stdlib_stats_distribution_normal + use stdlib_kinds + use stdlib_error, only : error_stop + use stdlib_stats_distribution_PRNG, only : dist_rand + use stdlib_stats_distribution_uniform, only : uni=>uniform_distribution_rvs + + implicit none + private + + real(dp), parameter :: HALF = 0.5_dp, ONE = 1.0_dp, TWO = 2.0_dp + integer, save :: kn(0:127) + real(dp), save :: wn(0:127), fn(0:127) + logical, save :: zig_norm_initialized = .false. + + public :: normal_distribution_rvs + public :: normal_distribution_pdf + public :: normal_distribution_cdf + + interface normal_distribution_rvs + !! Version experimental + !! + !! Normal Distribution Random Variates + !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! description)) + !! + module procedure norm_dist_rvs_0_rsp !0 dummy varaible + + #:for k1, t1 in RC_KINDS_TYPES + module procedure norm_dist_rvs_${t1[0]}$${k1}$ !2 dummy variables + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + module procedure norm_dist_rvs_array_${t1[0]}$${k1}$ !3 dummy variables + #:endfor + end interface normal_distribution_rvs + + interface normal_distribution_pdf + !! Version experimental + !! + !! Normal Distribution Probability Density Function + !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! description)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure norm_dist_pdf_${t1[0]}$${k1}$ + #:endfor + end interface normal_distribution_pdf + + interface normal_distribution_cdf + !! Version experimental + !! + !! Normal Distribution Cumulative Distribution Function + !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! description)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure norm_dist_cdf_${t1[0]}$${k1}$ + #:endfor + end interface normal_distribution_cdf + + + contains + + subroutine zigset + ! Marsaglia & Tsang generator for random normals & random exponentials. + ! Translated from C by Alan Miller (amiller@bigpond.net.au) + ! + ! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating + ! random variables', J. Statist. Software, v5(8). + ! + ! This is an electronic journal which can be downloaded from: + ! http://www.jstatsoft.org/v05/i08 + ! + ! N.B. It is assumed that all integers are 32-bit. + ! + ! Latest version - 1 January 2001 + ! + real(dp), parameter :: M1 = 2147483648.0_dp + real(dp) :: dn = 3.442619855899_dp, tn, & + vn = 0.00991256303526217_dp, q + integer :: i + + tn = dn + ! tables for random normals + q = vn * exp(HALF * dn * dn) + kn(0) = int((dn / q) * M1, kind = int32) + kn(1) = 0 + wn(0) = q / M1 + wn(127) = dn / M1 + fn(0) = ONE + fn(127) = exp( -HALF * dn * dn ) + do i = 126, 1, -1 + dn = sqrt( -TWO * log( vn / dn + exp( -HALF * dn * dn ) ) ) + kn(i+1) = int((dn / tn) * M1, kind = int32) + tn = dn + fn(i) = exp(-HALF * dn * dn) + wn(i) = dn / M1 + end do + zig_norm_initialized = .true. + return + end subroutine zigset + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function norm_dist_rvs_0_${t1[0]}$${k1}$( ) result(res) + ! Standard normal random vairate (0,1) + ! + ${t1}$ :: res + ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r + ${t1}$ :: x, y + integer :: hz, iz + + if( .not. zig_norm_initialized ) call zigset + iz = 0 + ! original algorithm use 32bit + hz = dist_rand(iz) + + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + res = hz * wn(iz) + else + L1: do + L2: if( iz == 0 ) then + do + x = -log( uni( ) ) * rr + y = -log( uni( ) ) + if( y + y >= x * x ) exit + end do + res = r + x + if( hz <= 0 ) res = -res + exit L1 + end if L2 + x = hz * wn(iz) + if( fn(iz) + uni( ) * (fn(iz-1) - fn(iz)) < & + exp(-HALF * x * x) ) then + res = x + exit L1 + end if + + !original algorithm use 32bit + hz = dist_rand(iz) + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + res = hz * wn(iz) + exit L1 + end if + end do L1 + end if + return + end function norm_dist_rvs_0_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function norm_dist_rvs_${t1[0]}$${k1}$(loc, scale) & + result(res) + ! Normal random variate (loc, scale) + ! + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res + ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r + ${t1}$ :: x, y + integer :: hz, iz + + if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & + //" parameter must be non-zero") + if( .not. zig_norm_initialized ) call zigset + iz = 0 + ! original algorithm use 32bit + hz = dist_rand(iz) + + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + res = hz * wn(iz) + else + L1: do + L2: if( iz == 0 ) then + do + x = -log( uni( ) ) * rr + y = -log( uni( ) ) + if( y + y >= x * x ) exit + end do + res = r + x + if( hz <= 0 ) res = -res + exit L1 + end if L2 + x = hz * wn(iz) + if( fn(iz) + uni( ) * (fn(iz-1) - fn(iz)) < & + exp(-HALF * x * x) ) then + res = x + exit L1 + end if + + !original algorithm use 32bit + hz = dist_rand(iz) + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + res = hz * wn(iz) + exit L1 + end if + end do L1 + end if + res = res * scale + loc + return + end function norm_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function norm_dist_rvs_${t1[0]}$${k1}$(loc, scale) & + result(res) + ! Normal distributed complex. The real part and imaginary part are & + ! independent of each other. + ! + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res + real(${k1}$) :: tr, ti + + tr = norm_dist_rvs_r${k1}$(real(loc), real(scale)) + ti = norm_dist_rvs_r${k1}$(aimag(loc), aimag(scale)) + res = cmplx(tr, ti) + return + end function norm_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + function norm_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & + result(res) + ${t1}$, intent(in) :: loc, scale + integer, intent(in) :: array_size + ${t1}$, allocatable :: res(:) + ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r + ${t1}$ :: x, y, re + integer :: hz, iz, i + + if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & + //" parameter must be non-zero") + if( .not. zig_norm_initialized ) call zigset + allocate(res(array_size)) + do i = 1, array_size + iz = 0 + ! original algorithm use 32bit + hz = dist_rand(iz) + + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + re = hz * wn(iz) + else + L1: do + L2: if( iz == 0 ) then + do + x = -log( uni( ) ) * rr + y = -log( uni( ) ) + if( y + y >= x * x ) exit + end do + re = r + x + if( hz <= 0 ) re = -re + exit L1 + end if L2 + x = hz * wn(iz) + if( fn(iz) + uni( ) * (fn(iz-1) - fn(iz)) < & + exp(-HALF * x * x) ) then + re = x + exit L1 + end if + + !original algorithm use 32bit + hz = dist_rand(iz) + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + re = hz * wn(iz) + exit L1 + end if + end do L1 + end if + res(i) = re * scale + loc + end do + return + end function norm_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + function norm_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & + result(res) + ${t1}$, intent(in) :: loc, scale + integer, intent(in) :: array_size + integer :: i + ${t1}$, allocatable :: res(:) + real(${k1}$) :: tr, ti + + allocate(res(array_size)) + do i = 1, array_size + tr = norm_dist_rvs_r${k1}$(real(loc), real(scale)) + ti = norm_dist_rvs_r${k1}$(aimag(loc), aimag(scale)) + res(i) = cmplx(tr, ti) + end do + return + end function norm_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function norm_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) & + result(res) + ! Normal distributed probability function + ! + ${t1}$, intent(in) :: x, loc, scale + real :: res + ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$)) + + if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & + //" parameter must be non-zero") + res = exp(- 0.5_${k1}$ * (x - loc) * (x - loc) / (scale * scale)) / & + (sqrt_2_Pi * scale) + return + end function norm_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function norm_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) & + result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + + res = norm_dist_pdf_r${k1}$(real(x), real(loc), real(scale)) + res = res * norm_dist_pdf_r${k1}$(aimag(x), aimag(loc), aimag(scale)) + return + end function norm_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function norm_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) & + result(res) + ! Normal random cumulative distribution function + ! + ${t1}$, intent(in) :: x, loc, scale + real :: res + ${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$) + + if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & + //" parameter must be non-zero") + res = (1.0_${k1}$ + erf((x - loc) / (scale * sqrt_2))) / 2.0_${k1}$ + return + end function norm_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function norm_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) & + result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + + res = norm_dist_cdf_r${k1}$(real(x), real(loc), real(scale)) + res = res * norm_dist_cdf_r${k1}$(aimag(x), aimag(loc), aimag(scale)) + return + end function norm_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + +end module stdlib_stats_distribution_normal \ No newline at end of file From 658aa44449736b4d808c834ba653d5c9b1411b0a Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sun, 20 Dec 2020 21:13:48 -0500 Subject: [PATCH 02/49] Delete stdlib_stats_distribution_normal.fypp --- stdlib_stats_distribution_normal.fypp | 366 -------------------------- 1 file changed, 366 deletions(-) delete mode 100644 stdlib_stats_distribution_normal.fypp diff --git a/stdlib_stats_distribution_normal.fypp b/stdlib_stats_distribution_normal.fypp deleted file mode 100644 index de55412ab..000000000 --- a/stdlib_stats_distribution_normal.fypp +++ /dev/null @@ -1,366 +0,0 @@ -#:include "common.fypp" -#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES -Module stdlib_stats_distribution_normal - use stdlib_kinds - use stdlib_error, only : error_stop - use stdlib_stats_distribution_PRNG, only : dist_rand - use stdlib_stats_distribution_uniform, only : uni=>uniform_distribution_rvs - - implicit none - private - - real(dp), parameter :: HALF = 0.5_dp, ONE = 1.0_dp, TWO = 2.0_dp - integer, save :: kn(0:127) - real(dp), save :: wn(0:127), fn(0:127) - logical, save :: zig_norm_initialized = .false. - - public :: normal_distribution_rvs - public :: normal_distribution_pdf - public :: normal_distribution_cdf - - interface normal_distribution_rvs - !! Version experimental - !! - !! Normal Distribution Random Variates - !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# - !! description)) - !! - module procedure norm_dist_rvs_0_rsp !0 dummy varaible - - #:for k1, t1 in RC_KINDS_TYPES - module procedure norm_dist_rvs_${t1[0]}$${k1}$ !2 dummy variables - #:endfor - - #:for k1, t1 in RC_KINDS_TYPES - module procedure norm_dist_rvs_array_${t1[0]}$${k1}$ !3 dummy variables - #:endfor - end interface normal_distribution_rvs - - interface normal_distribution_pdf - !! Version experimental - !! - !! Normal Distribution Probability Density Function - !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# - !! description)) - !! - #:for k1, t1 in RC_KINDS_TYPES - module procedure norm_dist_pdf_${t1[0]}$${k1}$ - #:endfor - end interface normal_distribution_pdf - - interface normal_distribution_cdf - !! Version experimental - !! - !! Normal Distribution Cumulative Distribution Function - !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# - !! description)) - !! - #:for k1, t1 in RC_KINDS_TYPES - module procedure norm_dist_cdf_${t1[0]}$${k1}$ - #:endfor - end interface normal_distribution_cdf - - - contains - - subroutine zigset - ! Marsaglia & Tsang generator for random normals & random exponentials. - ! Translated from C by Alan Miller (amiller@bigpond.net.au) - ! - ! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating - ! random variables', J. Statist. Software, v5(8). - ! - ! This is an electronic journal which can be downloaded from: - ! http://www.jstatsoft.org/v05/i08 - ! - ! N.B. It is assumed that all integers are 32-bit. - ! - ! Latest version - 1 January 2001 - ! - real(dp), parameter :: M1 = 2147483648.0_dp - real(dp) :: dn = 3.442619855899_dp, tn, & - vn = 0.00991256303526217_dp, q - integer :: i - - tn = dn - ! tables for random normals - q = vn * exp(HALF * dn * dn) - kn(0) = int((dn / q) * M1, kind = int32) - kn(1) = 0 - wn(0) = q / M1 - wn(127) = dn / M1 - fn(0) = ONE - fn(127) = exp( -HALF * dn * dn ) - do i = 126, 1, -1 - dn = sqrt( -TWO * log( vn / dn + exp( -HALF * dn * dn ) ) ) - kn(i+1) = int((dn / tn) * M1, kind = int32) - tn = dn - fn(i) = exp(-HALF * dn * dn) - wn(i) = dn / M1 - end do - zig_norm_initialized = .true. - return - end subroutine zigset - - #:for k1, t1 in REAL_KINDS_TYPES - impure elemental function norm_dist_rvs_0_${t1[0]}$${k1}$( ) result(res) - ! Standard normal random vairate (0,1) - ! - ${t1}$ :: res - ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r - ${t1}$ :: x, y - integer :: hz, iz - - if( .not. zig_norm_initialized ) call zigset - iz = 0 - ! original algorithm use 32bit - hz = dist_rand(iz) - - iz = iand( hz, 127 ) - if( abs( hz ) < kn(iz) ) then - res = hz * wn(iz) - else - L1: do - L2: if( iz == 0 ) then - do - x = -log( uni( ) ) * rr - y = -log( uni( ) ) - if( y + y >= x * x ) exit - end do - res = r + x - if( hz <= 0 ) res = -res - exit L1 - end if L2 - x = hz * wn(iz) - if( fn(iz) + uni( ) * (fn(iz-1) - fn(iz)) < & - exp(-HALF * x * x) ) then - res = x - exit L1 - end if - - !original algorithm use 32bit - hz = dist_rand(iz) - iz = iand( hz, 127 ) - if( abs( hz ) < kn(iz) ) then - res = hz * wn(iz) - exit L1 - end if - end do L1 - end if - return - end function norm_dist_rvs_0_${t1[0]}$${k1}$ - - #:endfor - - #:for k1, t1 in REAL_KINDS_TYPES - impure elemental function norm_dist_rvs_${t1[0]}$${k1}$(loc, scale) & - result(res) - ! Normal random variate (loc, scale) - ! - ${t1}$, intent(in) :: loc, scale - ${t1}$ :: res - ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r - ${t1}$ :: x, y - integer :: hz, iz - - if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & - //" parameter must be non-zero") - if( .not. zig_norm_initialized ) call zigset - iz = 0 - ! original algorithm use 32bit - hz = dist_rand(iz) - - iz = iand( hz, 127 ) - if( abs( hz ) < kn(iz) ) then - res = hz * wn(iz) - else - L1: do - L2: if( iz == 0 ) then - do - x = -log( uni( ) ) * rr - y = -log( uni( ) ) - if( y + y >= x * x ) exit - end do - res = r + x - if( hz <= 0 ) res = -res - exit L1 - end if L2 - x = hz * wn(iz) - if( fn(iz) + uni( ) * (fn(iz-1) - fn(iz)) < & - exp(-HALF * x * x) ) then - res = x - exit L1 - end if - - !original algorithm use 32bit - hz = dist_rand(iz) - iz = iand( hz, 127 ) - if( abs( hz ) < kn(iz) ) then - res = hz * wn(iz) - exit L1 - end if - end do L1 - end if - res = res * scale + loc - return - end function norm_dist_rvs_${t1[0]}$${k1}$ - - #:endfor - - #:for k1, t1 in CMPLX_KINDS_TYPES - impure elemental function norm_dist_rvs_${t1[0]}$${k1}$(loc, scale) & - result(res) - ! Normal distributed complex. The real part and imaginary part are & - ! independent of each other. - ! - ${t1}$, intent(in) :: loc, scale - ${t1}$ :: res - real(${k1}$) :: tr, ti - - tr = norm_dist_rvs_r${k1}$(real(loc), real(scale)) - ti = norm_dist_rvs_r${k1}$(aimag(loc), aimag(scale)) - res = cmplx(tr, ti) - return - end function norm_dist_rvs_${t1[0]}$${k1}$ - - #:endfor - - #:for k1, t1 in REAL_KINDS_TYPES - function norm_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & - result(res) - ${t1}$, intent(in) :: loc, scale - integer, intent(in) :: array_size - ${t1}$, allocatable :: res(:) - ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r - ${t1}$ :: x, y, re - integer :: hz, iz, i - - if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & - //" parameter must be non-zero") - if( .not. zig_norm_initialized ) call zigset - allocate(res(array_size)) - do i = 1, array_size - iz = 0 - ! original algorithm use 32bit - hz = dist_rand(iz) - - iz = iand( hz, 127 ) - if( abs( hz ) < kn(iz) ) then - re = hz * wn(iz) - else - L1: do - L2: if( iz == 0 ) then - do - x = -log( uni( ) ) * rr - y = -log( uni( ) ) - if( y + y >= x * x ) exit - end do - re = r + x - if( hz <= 0 ) re = -re - exit L1 - end if L2 - x = hz * wn(iz) - if( fn(iz) + uni( ) * (fn(iz-1) - fn(iz)) < & - exp(-HALF * x * x) ) then - re = x - exit L1 - end if - - !original algorithm use 32bit - hz = dist_rand(iz) - iz = iand( hz, 127 ) - if( abs( hz ) < kn(iz) ) then - re = hz * wn(iz) - exit L1 - end if - end do L1 - end if - res(i) = re * scale + loc - end do - return - end function norm_dist_rvs_array_${t1[0]}$${k1}$ - - #:endfor - - #:for k1, t1 in CMPLX_KINDS_TYPES - function norm_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & - result(res) - ${t1}$, intent(in) :: loc, scale - integer, intent(in) :: array_size - integer :: i - ${t1}$, allocatable :: res(:) - real(${k1}$) :: tr, ti - - allocate(res(array_size)) - do i = 1, array_size - tr = norm_dist_rvs_r${k1}$(real(loc), real(scale)) - ti = norm_dist_rvs_r${k1}$(aimag(loc), aimag(scale)) - res(i) = cmplx(tr, ti) - end do - return - end function norm_dist_rvs_array_${t1[0]}$${k1}$ - - #:endfor - - #:for k1, t1 in REAL_KINDS_TYPES - impure elemental function norm_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) & - result(res) - ! Normal distributed probability function - ! - ${t1}$, intent(in) :: x, loc, scale - real :: res - ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$)) - - if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & - //" parameter must be non-zero") - res = exp(- 0.5_${k1}$ * (x - loc) * (x - loc) / (scale * scale)) / & - (sqrt_2_Pi * scale) - return - end function norm_dist_pdf_${t1[0]}$${k1}$ - - #:endfor - - #:for k1, t1 in CMPLX_KINDS_TYPES - impure elemental function norm_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) & - result(res) - ${t1}$, intent(in) :: x, loc, scale - real :: res - - res = norm_dist_pdf_r${k1}$(real(x), real(loc), real(scale)) - res = res * norm_dist_pdf_r${k1}$(aimag(x), aimag(loc), aimag(scale)) - return - end function norm_dist_pdf_${t1[0]}$${k1}$ - - #:endfor - - #:for k1, t1 in REAL_KINDS_TYPES - impure elemental function norm_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) & - result(res) - ! Normal random cumulative distribution function - ! - ${t1}$, intent(in) :: x, loc, scale - real :: res - ${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$) - - if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & - //" parameter must be non-zero") - res = (1.0_${k1}$ + erf((x - loc) / (scale * sqrt_2))) / 2.0_${k1}$ - return - end function norm_dist_cdf_${t1[0]}$${k1}$ - - #:endfor - - #:for k1, t1 in CMPLX_KINDS_TYPES - impure elemental function norm_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) & - result(res) - ${t1}$, intent(in) :: x, loc, scale - real :: res - - res = norm_dist_cdf_r${k1}$(real(x), real(loc), real(scale)) - res = res * norm_dist_cdf_r${k1}$(aimag(x), aimag(loc), aimag(scale)) - return - end function norm_dist_cdf_${t1[0]}$${k1}$ - - #:endfor - -end module stdlib_stats_distribution_normal \ No newline at end of file From 51fb947e76a0b98888a2064b98bdff00c879b4fa Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sun, 20 Dec 2020 21:15:36 -0500 Subject: [PATCH 03/49] initial commit --- src/CMakeLists.txt | 6 + src/Makefile.manual | 42 +- src/stdlib_stats_distribution_PRNG.fypp | 211 +++++++++ src/stdlib_stats_distribution_normal.fypp | 366 ++++++++++++++++ src/stdlib_stats_distribution_uniform.fypp | 484 +++++++++++++++++++++ 5 files changed, 1094 insertions(+), 15 deletions(-) create mode 100644 src/stdlib_stats_distribution_PRNG.fypp create mode 100644 src/stdlib_stats_distribution_normal.fypp create mode 100644 src/stdlib_stats_distribution_uniform.fypp diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index ea7403663..33144ad07 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,6 +2,9 @@ # Create a list of the files to be preprocessed set(fppFiles + stdlib_bitsets.fypp + stdlib_bitsets_64.fypp + stdlib_bitsets_large.fypp stdlib_io.fypp stdlib_linalg.fypp stdlib_linalg_diag.fypp @@ -15,6 +18,9 @@ set(fppFiles stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp stdlib_quadrature_simps.fypp + stdlib_stats_distribution_PRNG.fypp + stdlib_stats_distribution_uniform.fypp + stdlib_stats_distribution_normal.fypp ) diff --git a/src/Makefile.manual b/src/Makefile.manual index 129dbe553..725b4aec9 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -1,5 +1,8 @@ SRC = f18estop.f90 \ stdlib_ascii.f90 \ + stdlib_bitsets.f90 \ + stdlib_bitsets_64.f90 \ + stdlib_bitsets_large.f90 \ stdlib_error.f90 \ stdlib_io.f90 \ stdlib_kinds.f90 \ @@ -13,10 +16,10 @@ SRC = f18estop.f90 \ stdlib_stats_mean.f90 \ stdlib_stats_moment.f90 \ stdlib_stats_var.f90 \ - stdlib_stats_distribution.f90 \ - stdlib_stats_distribution_rvs.f90 \ - stdlib_stats_distribution_implementation.f90 \ - + stdlib_stats_distribution_PRNG.f90\ + stdlib_stats_distribution_uniform.f90 \ + stdlib_stats_distribution_normal.f90 + LIB = libstdlib.a @@ -43,6 +46,9 @@ clean: # Fortran module dependencies f18estop.o: stdlib_error.o +stdlib_bitsets.o: stdlib_kinds.o +stdlib_bitsets_64.o: stdlib_bitsets.o +stdlib_bitsets_large.o: stdlib_bitsets.o stdlib_error.o: stdlib_optval.o stdlib_io.o: \ stdlib_error.o \ @@ -64,14 +70,21 @@ stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o -stdlib_stats_distribution_rvs.o: stdlib_kinds.o -stdlib_stats_distribution.o: \ - stdlib_error.o \ - stdlib_kinds.o \ - stdlib_stats_distribution_rvs.o \ -stdlib_stats_distribution_imp.o: \ - stdlib_stats_distribution.o +stdlib_stats_distribution_PRNG.o: stdlib_kinds.o +stdlib_stats_distribution_uniform.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_PRNG.o +stdlib_stats_distribution_normal.o \ + stdlib_kinds.o + stdlib_error.o \ + stdlib_stats_distribution.PRNG.o \ + stdlib_stats_distribution.uniform.o + # Fortran sources that are built from fypp templates +stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp +stdlib_bitsets_large.f90: stdlib_bitsets_large.fypp +stdlib_bitsets.f90: stdlib_bitsets.fypp stdlib_io.f90: stdlib_io.fypp stdlib_linalg.f90: stdlib_linalg.fypp stdlib_linalg_diag.f90: stdlib_linalg_diag.fypp @@ -80,7 +93,6 @@ stdlib_stats.f90: stdlib_stats.fypp stdlib_stats_mean.f90: stdlib_stats_mean.fypp stdlib_stats_moment.f90: stdlib_stats_moment.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp -stdlib_stats_distribution_rvs.f90: stdlib_stats_distribution_rvs.fypp -stdlib_stats_distribution.f90: stdlib_stats_distribution.fypp -stdlib_stats_distribution_implementation.f90: \ - stdlib_stats_distribution_implementation.fypp \ No newline at end of file +stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp +stdlib_stats_distribution_uniform.f90: stdlib_stats_distribution_uniform.fypp +stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp \ No newline at end of file diff --git a/src/stdlib_stats_distribution_PRNG.fypp b/src/stdlib_stats_distribution_PRNG.fypp new file mode 100644 index 000000000..90d9f367e --- /dev/null +++ b/src/stdlib_stats_distribution_PRNG.fypp @@ -0,0 +1,211 @@ +#:include "common.fypp" +module stdlib_stats_distribution_PRNG + use stdlib_kinds + implicit none + private + integer, parameter :: MAX_INT_BIT_SIZE = 64 + integer(int64), save :: st(4), si = 614872703977525537_int64 + logical, save :: seed_initialized = .false. + + public :: random_seed + public :: dist_rand + public :: jump + public :: long_jump + + + interface dist_rand + !! Version experimental + !! + !! Generation of random integers with different kinds + !! ([Specification](../page/specs/stdlib_stats_distribution_PRNG.html# + !! description)) + #:for k1, t1 in INT_KINDS_TYPES + module procedure dist_rand_${t1[0]}$${k1}$ + #:endfor + end interface dist_rand + + interface random_seed + !! Version experimental + !! + !! Set seed value for random number generator + !! ([Specification](../page/specs/stdlib_stats_distribution_PRNG.html# + !! description)) + !! + #:for k1, t1 in INT_KINDS_TYPES + module procedure random_distribution_seed_${t1[0]}$${k1}$ + #:endfor + end interface random_seed + + + contains + + #:for k1, t1 in INT_KINDS_TYPES + function dist_rand_${t1[0]}$${k1}$(n) result(res) + !! Random integer generation for various kinds + !! result = [-2^k, 2^k - 1], k = 7, 15, 31, 63, depending on input kind + !! Result is used as bit model number instead of regular arithmetic number + !! + ${t1}$, intent(in) :: n + ${t1}$ :: res + integer :: k + + k = MAX_INT_BIT_SIZE - bit_size(n) + res = shiftr(xoshiro256ss( ), k) + end function dist_rand_${t1[0]}$${k1}$ + + #:endfor + + function xoshiro256ss( ) result (res) + ! Generate random 64-bit integers + ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) + ! + ! This is xoshiro256** 1.0, one of our all-purpose, rock-solid + ! generators. It has excellent (sub-ns) speed, a state (256 bits) that is + ! large enough for any parallel application, and it passes all tests we + ! are aware of. + ! + ! The state must be seeded so that it is not everywhere zero. If you have + ! a 64-bit seed, we suggest to seed a splitmix64 generator and use its + ! output to fill st. + ! + ! Fortran 90 version translated from C by Jim-215-Fisher + ! + integer(int64) :: res, t + + if(.not. seed_initialized) call random_distribution_seed_iint64(si,t) + res = rol64(st(2) * 5 , 7) * 9 + t = shiftl(st(2), 17) + st(3) = ieor(st(3), st(1)) + st(4) = ieor(st(4), st(2)) + st(2) = ieor(st(2), st(3)) + st(1) = ieor(st(1), st(4)) + st(3) = ieor(st(3), t) + st(4) = rol64(st(4), 45) + return + end function xoshiro256ss + + function rol64(x, k) result(res) + integer(int64), intent(in) :: x + integer, intent(in) :: k + integer(int64) :: t1, t2, res + + t1 = shiftr(x, (64 - k)) + t2 = shiftl(x, k) + res = ior(t1, t2) + return + end function rol64 + + + subroutine jump + ! This is the jump function for the xoshiro256ss generator. It is equivalent + ! to 2^128 calls to xoshiro256ss(); it can be used to generate 2^128 + ! non-overlapping subsequences for parallel computations. + ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) + ! + ! Fortran 90 version translated from C by Jim-215-Fisher + integer(int64) :: jp(4) = [1733541517147835066_int64, & + -3051731464161248980_int64, & + -6244198995065845334_int64, & + 4155657270789760540_int64] + integer(int64) :: s1 = 0, s2 = 0, s3 = 0, s4 = 0, c = 1_int64 + integer :: i, j, k + + do i = 1, 4 + do j = 1, 64 + if(iand(jp(i), shiftl(c, j - 1)) /= 0) then + s1 = ieor(s1, st(1)) + s2 = ieor(s2, st(2)) + s3 = ieor(s3, st(3)) + s4 = ieor(s4, st(4)) + end if + k = xoshiro256ss( ) + end do + end do + st(1) = s1 + st(2) = s2 + st(3) = s3 + st(4) = s4 + end subroutine jump + + subroutine long_jump + ! This is the long-jump function for the xoshiro256ss generator. It is + ! equivalent to 2^192 calls to xoshiro256ss(); it can be used to generate + ! 2^64 starting points, from each of which jump() will generate 2^64 + ! non-overlapping subsequences for parallel distributed computations + ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) + ! + ! Fortran 90 version translated from C by Jim-215-Fisher + integer(int64) :: jp(4) = [8566230491382795199_int64, & + -4251311993797857357_int64, & + 8606660816089834049_int64, & + 4111957640723818037_int64] + integer(int64) :: s1 = 0, s2 = 0, s3 = 0, s4 = 0, c = 1_int64 + integer(int32) :: i, j, k + + do i = 1, 4 + do j = 1, 64 + if(iand(jp(i), shiftl(c, j - 1)) /= 0) then + s1 = ieor(s1, st(1)) + s2 = ieor(s2, st(2)) + s3 = ieor(s3, st(3)) + s4 = ieor(s4, st(4)) + end if + k = xoshiro256ss() + end do + end do + st(1) = s1 + st(2) = s2 + st(3) = s3 + st(4) = s4 + end subroutine long_jump + + function splitmix64(s) result(res) + ! Written in 2015 by Sebastiano Vigna (vigna@acm.org) + ! This is a fixed-increment version of Java 8's SplittableRandom + ! generator. + ! See http://dx.doi.org/10.1145/2714064.2660195 and + ! http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html + ! + ! It is a very fast generator passing BigCrush, and it can be useful if + ! for some reason you absolutely want 64 bits of state. + ! + ! Fortran 90 translated from C by Jim-215-Fisher + ! + integer(int64) :: res, int01, int02, int03 + integer(int64), intent(in), optional :: s + data int01, int02, int03/-7046029254386353131_int64, & + -4658895280553007687_int64, & + -7723592293110705685_int64/ + + if(present(s)) si = s + res = si + si = res + int01 + res = ieor(res, shiftr(res, 30)) * int02 + res = ieor(res, shiftr(res, 27)) * int03 + res = ieor(res, shiftr(res, 31)) + return + end function splitmix64 + + #:for k1, t1 in INT_KINDS_TYPES + subroutine random_distribution_seed_${t1[0]}$${k1}$(put, get) + !! Set seed value for random number generator + !! + ${t1}$, intent(in) :: put + ${t1}$, intent(out) :: get + integer(int64) :: tmp + integer :: i + + tmp = splitmix64(int(put, kind = int64)) + do i = 1, 10 + tmp = splitmix64( ) + end do + do i = 1, 4 + tmp = splitmix64( ) + st(i) = tmp + end do + get = int(tmp, kind = ${k1}$) + seed_initialized = .true. + end subroutine random_distribution_seed_${t1[0]}$${k1}$ + + #:endfor +end module stdlib_stats_distribution_PRNG \ No newline at end of file diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp new file mode 100644 index 000000000..de55412ab --- /dev/null +++ b/src/stdlib_stats_distribution_normal.fypp @@ -0,0 +1,366 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +Module stdlib_stats_distribution_normal + use stdlib_kinds + use stdlib_error, only : error_stop + use stdlib_stats_distribution_PRNG, only : dist_rand + use stdlib_stats_distribution_uniform, only : uni=>uniform_distribution_rvs + + implicit none + private + + real(dp), parameter :: HALF = 0.5_dp, ONE = 1.0_dp, TWO = 2.0_dp + integer, save :: kn(0:127) + real(dp), save :: wn(0:127), fn(0:127) + logical, save :: zig_norm_initialized = .false. + + public :: normal_distribution_rvs + public :: normal_distribution_pdf + public :: normal_distribution_cdf + + interface normal_distribution_rvs + !! Version experimental + !! + !! Normal Distribution Random Variates + !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! description)) + !! + module procedure norm_dist_rvs_0_rsp !0 dummy varaible + + #:for k1, t1 in RC_KINDS_TYPES + module procedure norm_dist_rvs_${t1[0]}$${k1}$ !2 dummy variables + #:endfor + + #:for k1, t1 in RC_KINDS_TYPES + module procedure norm_dist_rvs_array_${t1[0]}$${k1}$ !3 dummy variables + #:endfor + end interface normal_distribution_rvs + + interface normal_distribution_pdf + !! Version experimental + !! + !! Normal Distribution Probability Density Function + !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! description)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure norm_dist_pdf_${t1[0]}$${k1}$ + #:endfor + end interface normal_distribution_pdf + + interface normal_distribution_cdf + !! Version experimental + !! + !! Normal Distribution Cumulative Distribution Function + !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! description)) + !! + #:for k1, t1 in RC_KINDS_TYPES + module procedure norm_dist_cdf_${t1[0]}$${k1}$ + #:endfor + end interface normal_distribution_cdf + + + contains + + subroutine zigset + ! Marsaglia & Tsang generator for random normals & random exponentials. + ! Translated from C by Alan Miller (amiller@bigpond.net.au) + ! + ! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating + ! random variables', J. Statist. Software, v5(8). + ! + ! This is an electronic journal which can be downloaded from: + ! http://www.jstatsoft.org/v05/i08 + ! + ! N.B. It is assumed that all integers are 32-bit. + ! + ! Latest version - 1 January 2001 + ! + real(dp), parameter :: M1 = 2147483648.0_dp + real(dp) :: dn = 3.442619855899_dp, tn, & + vn = 0.00991256303526217_dp, q + integer :: i + + tn = dn + ! tables for random normals + q = vn * exp(HALF * dn * dn) + kn(0) = int((dn / q) * M1, kind = int32) + kn(1) = 0 + wn(0) = q / M1 + wn(127) = dn / M1 + fn(0) = ONE + fn(127) = exp( -HALF * dn * dn ) + do i = 126, 1, -1 + dn = sqrt( -TWO * log( vn / dn + exp( -HALF * dn * dn ) ) ) + kn(i+1) = int((dn / tn) * M1, kind = int32) + tn = dn + fn(i) = exp(-HALF * dn * dn) + wn(i) = dn / M1 + end do + zig_norm_initialized = .true. + return + end subroutine zigset + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function norm_dist_rvs_0_${t1[0]}$${k1}$( ) result(res) + ! Standard normal random vairate (0,1) + ! + ${t1}$ :: res + ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r + ${t1}$ :: x, y + integer :: hz, iz + + if( .not. zig_norm_initialized ) call zigset + iz = 0 + ! original algorithm use 32bit + hz = dist_rand(iz) + + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + res = hz * wn(iz) + else + L1: do + L2: if( iz == 0 ) then + do + x = -log( uni( ) ) * rr + y = -log( uni( ) ) + if( y + y >= x * x ) exit + end do + res = r + x + if( hz <= 0 ) res = -res + exit L1 + end if L2 + x = hz * wn(iz) + if( fn(iz) + uni( ) * (fn(iz-1) - fn(iz)) < & + exp(-HALF * x * x) ) then + res = x + exit L1 + end if + + !original algorithm use 32bit + hz = dist_rand(iz) + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + res = hz * wn(iz) + exit L1 + end if + end do L1 + end if + return + end function norm_dist_rvs_0_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function norm_dist_rvs_${t1[0]}$${k1}$(loc, scale) & + result(res) + ! Normal random variate (loc, scale) + ! + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res + ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r + ${t1}$ :: x, y + integer :: hz, iz + + if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & + //" parameter must be non-zero") + if( .not. zig_norm_initialized ) call zigset + iz = 0 + ! original algorithm use 32bit + hz = dist_rand(iz) + + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + res = hz * wn(iz) + else + L1: do + L2: if( iz == 0 ) then + do + x = -log( uni( ) ) * rr + y = -log( uni( ) ) + if( y + y >= x * x ) exit + end do + res = r + x + if( hz <= 0 ) res = -res + exit L1 + end if L2 + x = hz * wn(iz) + if( fn(iz) + uni( ) * (fn(iz-1) - fn(iz)) < & + exp(-HALF * x * x) ) then + res = x + exit L1 + end if + + !original algorithm use 32bit + hz = dist_rand(iz) + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + res = hz * wn(iz) + exit L1 + end if + end do L1 + end if + res = res * scale + loc + return + end function norm_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function norm_dist_rvs_${t1[0]}$${k1}$(loc, scale) & + result(res) + ! Normal distributed complex. The real part and imaginary part are & + ! independent of each other. + ! + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res + real(${k1}$) :: tr, ti + + tr = norm_dist_rvs_r${k1}$(real(loc), real(scale)) + ti = norm_dist_rvs_r${k1}$(aimag(loc), aimag(scale)) + res = cmplx(tr, ti) + return + end function norm_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + function norm_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & + result(res) + ${t1}$, intent(in) :: loc, scale + integer, intent(in) :: array_size + ${t1}$, allocatable :: res(:) + ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r + ${t1}$ :: x, y, re + integer :: hz, iz, i + + if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & + //" parameter must be non-zero") + if( .not. zig_norm_initialized ) call zigset + allocate(res(array_size)) + do i = 1, array_size + iz = 0 + ! original algorithm use 32bit + hz = dist_rand(iz) + + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + re = hz * wn(iz) + else + L1: do + L2: if( iz == 0 ) then + do + x = -log( uni( ) ) * rr + y = -log( uni( ) ) + if( y + y >= x * x ) exit + end do + re = r + x + if( hz <= 0 ) re = -re + exit L1 + end if L2 + x = hz * wn(iz) + if( fn(iz) + uni( ) * (fn(iz-1) - fn(iz)) < & + exp(-HALF * x * x) ) then + re = x + exit L1 + end if + + !original algorithm use 32bit + hz = dist_rand(iz) + iz = iand( hz, 127 ) + if( abs( hz ) < kn(iz) ) then + re = hz * wn(iz) + exit L1 + end if + end do L1 + end if + res(i) = re * scale + loc + end do + return + end function norm_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + function norm_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & + result(res) + ${t1}$, intent(in) :: loc, scale + integer, intent(in) :: array_size + integer :: i + ${t1}$, allocatable :: res(:) + real(${k1}$) :: tr, ti + + allocate(res(array_size)) + do i = 1, array_size + tr = norm_dist_rvs_r${k1}$(real(loc), real(scale)) + ti = norm_dist_rvs_r${k1}$(aimag(loc), aimag(scale)) + res(i) = cmplx(tr, ti) + end do + return + end function norm_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function norm_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) & + result(res) + ! Normal distributed probability function + ! + ${t1}$, intent(in) :: x, loc, scale + real :: res + ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$)) + + if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & + //" parameter must be non-zero") + res = exp(- 0.5_${k1}$ * (x - loc) * (x - loc) / (scale * scale)) / & + (sqrt_2_Pi * scale) + return + end function norm_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function norm_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) & + result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + + res = norm_dist_pdf_r${k1}$(real(x), real(loc), real(scale)) + res = res * norm_dist_pdf_r${k1}$(aimag(x), aimag(loc), aimag(scale)) + return + end function norm_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function norm_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) & + result(res) + ! Normal random cumulative distribution function + ! + ${t1}$, intent(in) :: x, loc, scale + real :: res + ${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$) + + if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & + //" parameter must be non-zero") + res = (1.0_${k1}$ + erf((x - loc) / (scale * sqrt_2))) / 2.0_${k1}$ + return + end function norm_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function norm_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) & + result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + + res = norm_dist_cdf_r${k1}$(real(x), real(loc), real(scale)) + res = res * norm_dist_cdf_r${k1}$(aimag(x), aimag(loc), aimag(scale)) + return + end function norm_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + +end module stdlib_stats_distribution_normal \ No newline at end of file diff --git a/src/stdlib_stats_distribution_uniform.fypp b/src/stdlib_stats_distribution_uniform.fypp new file mode 100644 index 000000000..d092129ce --- /dev/null +++ b/src/stdlib_stats_distribution_uniform.fypp @@ -0,0 +1,484 @@ +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES +#:set ALL_KINDS_TYPES = INT_KINDS_TYPES + RC_KINDS_TYPES +Module stdlib_stats_distribution_uniform + use stdlib_kinds + use stdlib_error, only : error_stop + use stdlib_stats_distribution_PRNG, only : dist_rand + + implicit none + private + + real(dp), parameter :: MESENNE_NUMBER = 1.0_dp / (2.0_dp ** 53 - 1.0_dp) + integer(int64), parameter :: INT_ONE = 1_int64 + + public :: uniform_distribution_rvs + public :: uniform_distribution_pdf + public :: uniform_distribution_cdf + public :: shuffle + + interface uniform_distribution_rvs + !! Version experimental + !! + !! Get uniformly distributed random variate for integer, real and complex + !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# + !! description)) + + module procedure unif_dist_rvs_0_rsp ! 0 dummy variable + + #:for k1, t1 in ALL_KINDS_TYPES + module procedure unif_dist_rvs_1_${t1[0]}$${k1}$ ! 1 dummy variable + #:endfor + + #:for k1, t1 in ALL_KINDS_TYPES + module procedure unif_dist_rvs_${t1[0]}$${k1}$ ! 2 dummy variables + #:endfor + + #:for k1, t1 in ALL_KINDS_TYPES + module procedure unif_dist_rvs_array_${t1[0]}$${k1}$ ! 3 dummy variables + #:endfor + end interface uniform_distribution_rvs + + interface uniform_distribution_pdf + !! Version experiment + !! + !! Get uniform distribution probability density (pdf) for integer, real and complex variables + !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# + !! description)) + + #:for k1, t1 in ALL_KINDS_TYPES + module procedure unif_dist_pdf_${t1[0]}$${k1}$ + #:endfor + end interface uniform_distribution_pdf + + interface uniform_distribution_cdf + !! Version experimental + !! + !! Get uniform distribution cumulative distribution function (cdf) for integer, real and complex variables + !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# + !! description)) + !! + #:for k1, t1 in ALL_KINDS_TYPES + module procedure unif_dist_cdf_${t1[0]}$${k1}$ + #:endfor + end interface uniform_distribution_cdf + + interface shuffle + !! Version experimental + !! + !! Fisher-Yates shuffle algorithm for a rank one array of integer, real and complex variables + !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# + !! description)) + !! + #:for k1, t1 in ALL_KINDS_TYPES + module procedure shuffle_${t1[0]}$${k1}$ + #:endfor + end interface shuffle + + + contains + + #:for k1, t1 in INT_KINDS_TYPES + impure elemental function unif_dist_rvs_1_${t1[0]}$${k1}$(scale) result(res) + ! Uniformly distributed integer in [0, scale] + ! Bitmask with rejection + ! https://www.pcg-random.org/posts/bounded-rands.html + ! + ! Fortran 90 translated from c by Jim-215-fisher + ${t1}$, intent(in) :: scale + ${t1}$ :: res, u, mask, n + integer :: zeros, bits_left, bits + + n = scale + if(n <= 0_${k1}$) call error_stop("Error: Uniform distribution scale" & + //" parameter must be positive") + zeros = leadz(n) + bits = bit_size(n) - zeros + mask = shiftr(not(0_${k1}$), zeros) + L1 : do + u = dist_rand(n) + res = iand(u, mask) + if(res <= n) exit L1 + bits_left = zeros + L2 : do + if(bits_left < bits) exit L2 + u = shiftr(u, bits) + res = iand(u, mask) + if(res <= n) exit L1 + bits_left = bits_left - bits + end do L2 + end do L1 + return + end function unif_dist_rvs_1_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + impure elemental function unif_dist_rvs_${t1[0]}$${k1}$(loc, scale) & + result( res ) + ! Uniformly distributed integer in [loc, loc + scale] + ! + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res + + if(scale == 0_${k1}$) call error_stop("Error: Uniform distribution" & + //" scale parameter must be non-zero") + res = loc + unif_dist_rvs_1_${t1[0]}$${k1}$(scale) + return + end function unif_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function unif_dist_rvs_0_${t1[0]}$${k1}$( ) result(res) + ! Uniformly distributed float in [0,1] + ! Based on the paper by Frederic Goualard, "Generating Random Floating- + ! Point Numbers By Dividing Integers: a Case Study", Proceedings of + ! ICCS 2020, June 20202, Amsterdam, Netherlands + ! + ${t1}$ :: res + integer(int64) :: tmp + + tmp = shiftr(dist_rand(INT_ONE), 11) ! Get random from [0,2^53-1] + res = real(tmp * MESENNE_NUMBER, kind =${k1}$) ! convert to [0,1] + return + end function unif_dist_rvs_0_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function unif_dist_rvs_1_${t1[0]}$${k1}$(scale) result(res) + ! Uniformly distributed float in [0, scale] + ! + ${t1}$, intent(in) :: scale + ${t1}$ :: res + + if(scale == 0._${k1}$) call error_stop("Error: Uniform distribution" & + //" scale parameter must be non-zero") + res = scale * unif_dist_rvs_0_${t1[0]}$${k1}$( ) + return + end function unif_dist_rvs_1_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + impure elemental function unif_dist_rvs_${t1[0]}$${k1}$(loc, scale) & + result(res) + ! Uniformly distributed float in [loc, loc + scale] + ! + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res + + if(scale == 0._${k1}$) call error_stop("Error: Uniform distribution" & + //" scale parameter must be non-zero") + res = loc + scale * unif_dist_rvs_0_${t1[0]}$${k1}$( ) + return + end function unif_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function unif_dist_rvs_1_${t1[0]}$${k1}$(scale) result(res) + ! Uniformly distributed complex in [(0,0i), (scale, i(scale)] + ! The real part and imaginary part are independent of each other, so that + ! the joint distribution is on an unit square [(0,0i), scale,i(scale)] + ! + ${t1}$, intent(in) :: scale + ${t1}$ :: res + integer(int64) :: tmp + real(${k1}$) :: r1, r2, tr, ti + + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & + //" distribution scale parameter must be non-zero") + r1 = unif_dist_rvs_0_r${k1}$( ) + if(real(scale) == 0.0_${k1}$) then + ti = aimag(scale) * r1 + tr = 0.0_${k1}$ + elseif(aimag(scale) == 0.0_${k1}$) then + tr = real(scale) * r1 + ti = 0.0_${k1}$ + else + r2 = unif_dist_rvs_0_r${k1}$( ) + tr = real(scale) * r1 + ti = aimag(scale) * r2 + endif + res = cmplx(tr, ti) + return + end function unif_dist_rvs_1_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + impure elemental function unif_dist_rvs_${t1[0]}$${k1}$(loc, scale) & + result(res) + ! Uniformly distributed complex in [(loc,iloc), (loc + scale, i(loc + scale)] + ! The real part and imaginary part are independent of each other, so that + ! the joint distribution is on an unit square [(loc,iloc), (loc + scale, + ! i(loc + scale))] + ! + ${t1}$, intent(in) :: loc, scale + ${t1}$ :: res + integer(int64) :: tmp + real(${k1}$) :: r1, r2, tr, ti + + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & + //" distribution scale parameter must be non-zero") + r1 = unif_dist_rvs_0_r${k1}$( ) + if(real(scale) == 0.0_${k1}$) then + tr = real(loc) + ti = aimag(loc) + aimag(scale) * r1 + elseif(aimag(scale) == 0.0_${k1}$) then + tr = real(loc) + real(scale) * r1 + ti = aimag(loc) + else + r2 = unif_dist_rvs_0_r${k1}$( ) + tr = real(loc) + real(scale) * r1 + ti = aimag(loc) + aimag(scale) * r2 + endif + res = cmplx(tr, ti) + return + end function unif_dist_rvs_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + function unif_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & + result(res) + ${t1}$, intent(in) :: loc, scale + ${t1}$, allocatable :: res(:) + ${t1}$ :: u, mask, n, nn + integer, intent(in) :: array_size + integer :: i, zeros, bits_left, bits + + n = scale + if(n == 0_${k1}$) call error_stop("Error: Uniform distribution" & + //" scale parameter must be non-zero") + allocate(res(array_size)) + zeros = leadz(n) + bits = bit_size(n) - zeros + mask = shiftr(not(0_${k1}$), zeros) + do i = 1, array_size + L1 : do + u = dist_rand(n) + nn = iand(u, mask) + if(nn <= n) exit L1 + bits_left = zeros + L2 : do + if(bits_left < bits) exit L2 + u = shiftr(u, bits) + nn = iand(u, mask) + if(nn <= n) exit L1 + bits_left = bits_left - bits + end do L2 + end do L1 + res(i) = loc + nn + end do + return + end function unif_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + function unif_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & + result(res) + ${t1}$, intent(in) :: loc, scale + ${t1}$, allocatable :: res(:) + ${t1}$ :: t + integer, intent(in) :: array_size + integer(int64) :: tmp + integer :: i + + + if(scale == 0._${k1}$) call error_stop("Error: Uniform distribution" & + //" scale parameter must be non-zero") + allocate(res(array_size)) + do i = 1, array_size + tmp = shiftr(dist_rand(INT_ONE), 11) + t = real(tmp * MESENNE_NUMBER, kind = ${k1}$) + res(i) = loc + scale * t + enddo + return + end function unif_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + function unif_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & + result(res) + ${t1}$, intent(in) :: loc, scale + ${t1}$, allocatable :: res(:) + real(${k1}$) :: r1, r2, tr, ti + integer, intent(in) :: array_size + integer(int64) :: tmp + integer :: i + + + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & + //" distribution scale parameter must be non-zero") + allocate(res(array_size)) + do i = 1, array_size + tmp = shiftr(dist_rand(INT_ONE), 11) + r1 = real(tmp * MESENNE_NUMBER, kind = ${k1}$) + if(real(scale) == 0.0_${k1}$) then + tr = real(loc) + ti = aimag(loc) + aimag(scale) * r1 + elseif(aimag(scale) == 0.0_${k1}$) then + tr = real(loc) + real(scale) * r1 + ti = aimag(loc) + else + tmp = shiftr(dist_rand(INT_ONE), 11) + r2 = real(tmp * MESENNE_NUMBER, kind = ${k1}$) + tr = real(loc) + real(scale) * r1 + ti = aimag(loc) + aimag(scale) * r2 + endif + res(i) = cmplx(tr, ti) + enddo + return + end function unif_dist_rvs_array_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + elemental function unif_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + + if(scale == 0) then + res = 0.0 + elseif(x < loc .or. x >loc + scale) then + res = 0.0 + else + res = 1. / (scale + 1) + end if + return + end function unif_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + elemental function unif_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + + if(scale == 0.0_${k1}$) then + res = 0.0 + elseif(x <= loc .or. x >= (loc + scale)) then + res = 0.0 + else + res = 1.0 / scale + end if + return + end function unif_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + elemental function unif_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + real(${k1}$) :: tr, ti + + tr = real(loc) + real(scale); ti = aimag(loc) + aimag(scale) + if(scale == (0.0_${k1}$,0.0_${k1}$)) then + res = 0.0 + elseif((real(x) >= real(loc) .and. real(x) <= tr) .and. & + (aimag(x) >= aimag(loc) .and. aimag(x) <= ti)) then + res = 1.0 / (real(scale) * aimag(scale)) + else + res = 0.0 + end if + return + end function unif_dist_pdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in INT_KINDS_TYPES + elemental function unif_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + + if(scale == 0) then + res = 0.0 + elseif(x < loc) then + res = 0.0 + elseif(x >= loc .and. x <= (loc + scale)) then + res = real((x - loc + 1)) / real((scale + 1)) + else + res = 1.0 + end if + return + end function unif_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES + elemental function unif_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + + if(scale == 0.0_${k1}$) then + res = 0.0 + elseif(x < loc) then + res = 0.0 + elseif(x >= loc .and. x <= (loc + scale)) then + res = (x - loc) / scale + else + res = 1.0 + end if + return + end function unif_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES + elemental function unif_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) result(res) + ${t1}$, intent(in) :: x, loc, scale + real :: res + logical :: r1, r2, i1, i2 + + if(scale == (0.0_${k1}$,0.0_${k1}$)) then + res = 0.0 + return + endif + r1 = real(x) < real(loc) + r2 = real(x) > (real(loc) + real(scale)) + i1 = aimag(x) < aimag(loc) + i2 = aimag(x) > (aimag(loc) + aimag(scale)) + if(r1 .or. i1) then + res = 0.0 + elseif((.not. r1) .and. (.not. r2) .and. i2) then + res = (real(x) - real(loc)) / real(scale) + elseif((.not. i1) .and. (.not. i2) .and. r2) then + res = (aimag(x) - aimag(loc)) / aimag(scale) + elseif((.not. r1) .and. (.not. r2) .and. (.not. i1) .and. (.not. i2)) & + then + res = (real(x) - real(loc)) * (aimag(x) - aimag(loc)) / & + (real(scale) * aimag(scale)) + elseif(r2 .and. i2)then + res = 1.0 + end if + return + end function unif_dist_cdf_${t1[0]}$${k1}$ + + #:endfor + + #:for k1, t1 in ALL_KINDS_TYPES + function shuffle_${t1[0]}$${k1}$( list ) result(res) + ${t1}$, intent(in) :: list(:) + ${t1}$, allocatable :: res(:) + ${t1}$ :: tmp + integer :: n, i, j + + n = size(list) + allocate(res(n), source=list) + do i = 1, n - 1 + j = uniform_distribution_rvs(n - i) + i + tmp = res(i) + res(i) = res(j) + res(j) = tmp + end do + return + end function shuffle_${t1[0]}$${k1}$ + + #:endfor +end module stdlib_stats_distribution_uniform \ No newline at end of file From ce8731452a07f6003536861ad6d73e0992cd83a3 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sun, 20 Dec 2020 21:16:19 -0500 Subject: [PATCH 04/49] initial commit --- doc/specs/index.md | 4 + doc/specs/stdlib_stats_distribution_normal.md | 247 ++++++++++++++++++ 2 files changed, 251 insertions(+) create mode 100644 doc/specs/stdlib_stats_distribution_normal.md diff --git a/doc/specs/index.md b/doc/specs/index.md index 91284c2df..a827b392c 100644 --- a/doc/specs/index.md +++ b/doc/specs/index.md @@ -17,7 +17,11 @@ This is and index/directory of the specifications (specs) for each new module/fe - [optval](./stdlib_optval.html) - Fallback value for optional arguments - [quadrature](./stdlib_quadrature.html) - Numerical integration - [stats](./stdlib_stats.html) - Descriptive Statistics + - [stats_distribution_PRNG](./stdlib_stats_distribution_PRNG.html) - Probability Distributions random number generator + - [stats_distribution_uniform](./stdlib_stats_distribution_uniform.html) - Uniform probability distribution + - [stats_distribution_normal](./stdlib_stats_distribution_normal.html) - Normal probability distribution + ## Missing specs - [ascii](https://github.com/fortran-lang/stdlib/blob/master/src/stdlib_ascii.f90) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md new file mode 100644 index 000000000..cefbaef87 --- /dev/null +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -0,0 +1,247 @@ +--- +title: stats_distribution +--- + +# Statistical Distributions Normal Module + +[TOC] + +## `normal_distribution_rvs` - normal distribution random variates + +### Status + +Experimental + +### Description + +A normal continuous random variate distribution, also known as Gaussian, or Gauss or Laplace-Gauss distribution. The location `loc` specifies the mean or expectation. The `scale` specifies the standard deviation. + +Without augument the function returns a standard normal distributed random variate N(0,1). The function is elemental. + +With two arguments, the function returns a normal distributed random variate N(loc, scale^2). For complex auguments, the real and imaginary parts are independent of each other. The function is elemental. + +With three auguments, the function returns a rank one array of normal distributed random variates. + +### Syntax + +`result = [[stdlib_stats_distribution_normal(module):normal_distribution_rvs(interface)]]([loc, scale] [[, array_size]])` + +### Arguments + +`array_size`: optional argument has `intent(in)` and is a scalar of type `integer`. + +`loc`: optional argument has `intent(in)` and is a scalar of type `real` or `complx`. + +`scale`: optional argument has `intent(in)` and is a scalar of type `real` or `complx`. + +`loc` and `scale` arguments must have the same type. + +### Return value + +The result is a scalar or rank one array, with a size of `array_size`, of type `real` or `complx`. + +### Example + +```fortran +program demo_normal_rvs + use stdlib_stats_distribution_normal, only: norm => normal_distribution_rvs + use stdlib_stats_distribution_PRNG, only: random_seed + implicit none + real :: a(2,3,4), b(2,3,4) + complx :: loc, scale + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + print *, norm( ) !single standard normal random variate + +! 0.563655198 + + print *, norm(1.0, 2.0) + !normal random variate \(\mu\)=1.0, \(\sigma\)=2.0 + +! -0.633261681 + + print *, norm(0.,1.0,10) !an array of 10 standard norml random variates + +! [-3.38123664E-02, -0.190365672, 0.220678389, -0.424612164, -0.249541596, +! 0.865260184, 1.11086845, -0.328349441, 1.10873628, 1.27049923] + + a(:,:,:) = 1.0 + b(:,:,:) = 1.0 + print *, norm(a,b) ! a rank 3 random variates array +![0.152776539, -7.51764774E-02, 1.47208166, 0.180561781, 1.32407105, +! 1.20383692, 0.123445868, -0.455737948, -0.469808221, 1.60750175, +! 1.05748117, 0.720934749, 0.407810807, 1.48165631, 2.31749439, +! 0.414566994, 3.06084275, 1.86505437, 1.36338580, 7.26878643E-02, +! 0.178585172, 1.39557445, 0.828021586, 0.872084975] + + loc = (-1.0, 2.0) + scale = (2.0, 1.0) + print *, norm(loc, scale) + !single complex normal random variate with real part of \(\mu\)=-1, \(\sigma\)=2; imagainary part of \(\mu\)=2.0 and \(\sigma\)=1.0 + +! (1.22566295,2.12518454) + +end program demo_normal_rvs +``` + +## `normal_distribution_pdf` - normal probability density function + +### Status + +Experimental + +### Description + +The probability density function of the continuous normal distribution. + +$$f(x)=\frac{1}{\sigma&space;\sqrt{2&space;\pi}}&space;e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^{2}})$$ + +### Syntax + +`result = [[stdlib_stats_distribution_normal(module):normal_distribution_pdf(interface)]](x, loc, scale)` + +### Arguments + +`x`: has `intent(in)` and is a scalar of type `real` or `complx`. + +`loc`: has `intent(in)` and is a scalar of type `real` or `complx`. + +`scale`: has `intent(in)` and is a scalar of type `real` or `complx`. + +The function is elemental, i.e., all three auguments could be arrays conformable to each other. All three arguments must have the same type. + +### Return value + +The result is a scalar or an array, with a shape conformable to auguments, of type `real`. + +### Example + +```fortran +program demo_normal_pdf + use stdlib_stats_distribution_PRNG, only : random_seed + use stdlib_stats_distribution_normal, only : & + norm_pdf=>normal_distribution_pdf,& + norm => normal_distribution_rvs + + implicit none + real :: x(3,4,5),a(3,4,5),b(3,4,5) + complx :: loc, scale + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + print *, norm_pdf(1.0,0.,1.) !a probability density at 1.0 in standard normal + +! 0.241970733 + + print *, norm_pdf(2.0,-1.0, 2.0) + !a probability density at 2.0 with \(\mu\)=-1.0 \(\sigma\)=2.0 + +!6.47588000E-02 + + x = reshape(norm(0.0, 1.0, 60),[3,4,5]) + ! standard normal random variates array + + a(:,:,:) = 0.0 + b(:,:,:) = 1.0 + print *, norm_pdf(x, a, b) ! standard normal probability density array + +! [0.340346158, 0.285823315, 0.398714304, 0.391778737, 0.389345556, +! 0.364551932, 0.386712372, 0.274370432, 0.215250477, 0.378006011, +! 0.215760440, 0.177990928, 0.278640658, 0.223813817, 0.356875211, +! 0.285167664, 0.378533930, 0.390739858, 0.271684974, 0.138273031, +! 0.135456234, 0.331718773, 0.398283750, 0.383706540] + + loc = (1.0, -0.5) + scale = (1.0, 2.) + print *, norm_pdf((1.5,1.0), loc, scale) + ! a complex normal probability density function at (1.5,1.0) with real part of \(\mu\)=1.0, \(\sigma\)=1.0 and imaginary part of \(\mu\)=-0.5, \(\sigma\)=2.0 + +! 5.30100204E-02 + +end program demo_normal_pdf +``` + +## `normal_distribution_cdf` - normal cumulative distribution function + +### Status + +Experimental + +### Description + +Cumulative distribution function of the normal continuous distribution + +$$F(X)=\frac{1}{2}\left&space;[&space;1&space;+&space;erf(\frac{x-\mu}{\sqrt{2}&space;\sigma})&space;\right&space;])$$ + +### Syntax + +`result = [[stdlib_stats_distribution_normal(module):normal_distribution_cdf(interface)]](x, loc, scale)` + +### Arguments + +`x`: has `intent(in)` and is a scalar of type `real` or `complx`. + +`loc`: has `intent(in)` and is a scalar of type `real` or `complx`. + +`scale`: has `intent(in)` and is a scalar of type `real` or `complx`. + +The function is elemental, i.e., all three auguments could be arrays conformable to each other. All three arguments must have the same type. + +### Return value + +The result is a scalar or an array, with a shape conformable to auguments, of type `real`. + +### Example + +```fortran +program demo_norm_cdf + use stdlib_stats_distribution_PRNG, only : random_seed + use stdlib_stats_distribution_normal, only : & + norm_cdf => normal_distribution_cdf, & + norm => normal_distribution_rvs + + implicit none + real :: x(2,3,4),a(2,3,4),b(2,3,4) + complx :: loc, scale + integer :: seed_put, seed_get + + seed_put = 1234567 + call random_seed(seed_put, seed_get) + + print *, norm_cdf(1.0, 0.0, 1.0) ! a standard normal cumulative at 1.0 + +! 0.841344714 + + print *, norm_cdf(2.0, -1.0, 2.0) + ! a cumulative at 2.0 with \(\mu\)=-1 \(\sigma\)=2 + +! 0.933192849 + + x = reshape(norm(0.0, 1.0, 24),[2,3,4]) + ! standard normal random variates array + + a(:,:,:) = 0.0 + b(:,:,:) = 1.0 + print *, norm_cdf(x, a, b) ! standard normal cumulative array + +! [0.713505626, 0.207069695, 0.486513376, 0.424511284, 0.587328553, +! 0.335559726, 0.401470929, 0.806552052, 0.866687536, 0.371323735, +! 0.866228044, 0.898046613, 0.198435277, 0.141147852, 0.681565762, +! 0.206268221, 0.627057910, 0.580759525, 0.190364420, 7.27325380E-02, +! 7.08068311E-02, 0.728241026, 0.522919059, 0.390097380] + + loc = (1.0,0.0) + scale = (0.5,1.0) + print *, norm_cdf((0.5,-0.5),loc,scale) + !complex normal cumulative distribution at (0.5,-0.5) with real part of \(\mu\)=1.0, \(\sigma\)=0.5 and imaginary part of \(\mu\)=0.0, \(\sigma\)=1.0 + +!4.89511043E-02 + +end program demo_norm_cdf + +``` From 86aa8ce1e415c8e6c112eb4541fc4ff71604a8e5 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sun, 20 Dec 2020 21:17:07 -0500 Subject: [PATCH 05/49] initial commit --- src/tests/stats/CMakeLists.txt | 1 + src/tests/stats/Makefile.manual | 6 +- src/tests/stats/test_distribution_normal.f90 | 557 +++++++++++++++++++ 3 files changed, 562 insertions(+), 2 deletions(-) create mode 100644 src/tests/stats/test_distribution_normal.f90 diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 36ffc7aeb..3fbf70055 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -5,6 +5,7 @@ ADDTEST(moment) ADDTEST(rawmoment) ADDTEST(var) ADDTEST(varn) +ADDTEST(distribution_normal) if(DEFINED CMAKE_MAXIMUM_RANK) if(${CMAKE_MAXIMUM_RANK} GREATER 7) diff --git a/src/tests/stats/Makefile.manual b/src/tests/stats/Makefile.manual index aacaf98ab..9a82079e0 100644 --- a/src/tests/stats/Makefile.manual +++ b/src/tests/stats/Makefile.manual @@ -1,3 +1,5 @@ -PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90 -include ../Makefile.manual.test.mk +PROGS_SRC = test_mean.f90 test_moment.f90 test_var.f90 \ + test_distribution_normal.f90 + +include ../Makefile.manual.test.mk \ No newline at end of file diff --git a/src/tests/stats/test_distribution_normal.f90 b/src/tests/stats/test_distribution_normal.f90 new file mode 100644 index 000000000..fbb0bdaf7 --- /dev/null +++ b/src/tests/stats/test_distribution_normal.f90 @@ -0,0 +1,557 @@ +program test_distribution_normal + use stdlib_kinds + use stdlib_error, only : check + use stdlib_stats_distribution_PRNG, only : random_seed + use stdlib_stats_distribution_normal, nor_rvs => normal_distribution_rvs, & + nor_pdf => normal_distribution_pdf, & + nor_cdf => normal_distribution_cdf + + implicit none + real(sp), parameter :: sptol = 1000 * epsilon(1.0_sp) + real(dp), parameter :: dptol = 1000 * epsilon(1.0_dp) + real(qp), parameter :: qptol = 1000 * epsilon(1.0_qp) + logical :: warn = .true. + integer :: put, get + real :: x(2,3,4),a(2,3,4), b(2,3,4) + complex :: loc, scale + + put = 12345678 + call random_seed(put, get) + + call test_normal_random_generator + + call test_nor_rvs_rsp + call test_nor_rvs_rdp + call test_nor_rvs_rqp + call test_nor_rvs_csp + call test_nor_rvs_cdp + call test_nor_rvs_cqp + + call test_nor_pdf_rsp + call test_nor_pdf_rdp + call test_nor_pdf_rqp + call test_nor_pdf_csp + call test_nor_pdf_cdp + call test_nor_pdf_cqp + + call test_nor_cdf_rsp + call test_nor_cdf_rdp + call test_nor_cdf_rqp + call test_nor_cdf_csp + call test_nor_cdf_cdp + call test_nor_cdf_cqp + stop + + + contains + + + subroutine test_normal_random_generator + integer :: i, j, freq(0:1000), num=10000000 + real(dp) :: chisq, expct + + print *, "" + print *, "Test normal random generator with chi-squared" + freq = 0 + do i = 1, num + j = 1000 * (1 + erf(nor_rvs(0.0, 1.0) / sqrt(2.0))) / 2.0 + freq(j) = freq(j) + 1 + end do + chisq = 0.0_dp + expct = num / 1000 + do i = 0, 999 + chisq = chisq + (freq(i) - expct) ** 2 / expct + end do + write(*,*) "The critical values for chi-squared with 1000 d. of f. is" & + //" 1143.92" + write(*,*) "Chi-squared for normal random generator is : ", chisq + call check((chisq < 1143.9), & + msg="normal randomness failed chi-squared test", warn=warn) + end subroutine test_normal_random_generator + + + subroutine test_nor_rvs_rsp + real(sp) :: res(10), loc, scale + integer :: i, n, k = 5 + integer :: seed, get + real(sp) :: ans(10) = [2.66708039318040679432897377409972250_sp, & + 2.36030794936128329730706809641560540_sp, & + 1.27712218793084242296487218482070602_sp, & + -2.39132544130814794769435138732660562_sp, & + 1.72566595106028652928387145948363468_sp, & + -1.50621775537767632613395107910037041_sp, & + 2.13518835158352082714827702147886157_sp, & + -0.636788253742142318358787633769679815_sp, & + 2.48600787778845799813609573902795091_sp, & + -3.03711473837981227319460231228731573_sp] + + print *, "Test normal_distribution_rvs_rsp" + seed = 25836914 + call random_seed(seed, get) + + loc = 0.5_sp; scale = 2.0_sp + do i = 1, k + res(i) = nor_rvs(loc, scale) ! 2 dummies + end do + res(6:10) = nor_rvs(loc, scale, k) ! 3 dummies + call check(all(abs(res - ans) < sptol), & + msg="normal_distribution_rvs_rsp failed", warn=warn) + end subroutine test_nor_rvs_rsp + + subroutine test_nor_rvs_rdp + real(dp) :: res(10), loc, scale + integer :: i, n, k = 5 + integer :: seed, get + real(dp) :: ans(10) = [2.66708039318040679432897377409972250_dp, & + 2.36030794936128329730706809641560540_dp, & + 1.27712218793084242296487218482070602_dp, & + -2.39132544130814794769435138732660562_dp, & + 1.72566595106028652928387145948363468_dp, & + -1.50621775537767632613395107910037041_dp, & + 2.13518835158352082714827702147886157_dp, & + -0.636788253742142318358787633769679815_dp, & + 2.48600787778845799813609573902795091_dp, & + -3.03711473837981227319460231228731573_dp] + + print *, "Test normal_distribution_rvs_rdp" + seed = 25836914 + call random_seed(seed, get) + + loc = 0.5_dp; scale = 2.0_dp + do i = 1, k + res(i) = nor_rvs(loc, scale) ! 2 dummies + end do + res(6:10) = nor_rvs(loc, scale, k) ! 3 dummies + call check(all(abs(res - ans) < dptol), & + msg="normal_distribution_rvs_rdp failed", warn=warn) + end subroutine test_nor_rvs_rdp + + subroutine test_nor_rvs_rqp + real(qp) :: res(10), loc, scale + integer :: i, n, k = 5 + integer :: seed, get + real(qp) :: ans(10) = [2.66708039318040679432897377409972250_qp, & + 2.36030794936128329730706809641560540_qp, & + 1.27712218793084242296487218482070602_qp, & + -2.39132544130814794769435138732660562_qp, & + 1.72566595106028652928387145948363468_qp, & + -1.50621775537767632613395107910037041_qp, & + 2.13518835158352082714827702147886157_qp, & + -0.636788253742142318358787633769679815_qp, & + 2.48600787778845799813609573902795091_qp, & + -3.03711473837981227319460231228731573_qp] + + print *, "Test normal_distribution_rvs_rqp" + seed = 25836914 + call random_seed(seed, get) + + loc = 0.5_qp; scale = 2.0_qp + do i = 1, k + res(i) = nor_rvs(loc, scale) ! 2 dummies + end do + res(6:10) = nor_rvs(loc, scale, k) ! 3 dummies + call check(all(abs(res - ans) < qptol), & + msg="normal_distribution_rvs_rqp failed", warn=warn) + end subroutine test_nor_rvs_rqp + + subroutine test_nor_rvs_csp + complex(sp) :: res(10), loc, scale + integer :: i, n, k = 5 + integer :: seed, get + complex(sp) :: ans(10) = [(2.12531018257141113281250000000000000_sp, & + 1.46507704257965087890625000000000000_sp), & + (1.08284163475036621093750000000000000_sp, & + 0.277168631553649902343750000000000000_sp), & + (1.41924941539764404296875000000000000_sp, & + 0.498445570468902587890625000000000000_sp), & + (1.72639131546020507812500000000000000_sp, & + 0.715802907943725585937500000000000000_sp), & + (1.98950588703155517578125000000000000_sp, & + 0.115721315145492553710937500000000000_sp), & + (-1.16929018497467041015625000000000000_sp, & + 0.250744730234146118164062500000000000_sp), & + (1.57160544395446777343750000000000000_sp, & + 0.638282597064971923828125000000000000_sp), & + (-1.36106109619140625000000000000000000_sp, & + 0.166259199380874633789062500000000000_sp), & + (1.13403093814849853515625000000000000_sp, & + 1.04232621192932128906250000000000000_sp), & + (-1.68220531940460205078125000000000000_sp, & + 1.63361442089080810546875000000000000_sp)] + + print *, "Test normal_distribution_rvs_csp" + seed = 25836914 + call random_seed(seed, get) + + loc = (0.5_sp, 1.0_sp); scale = (1.5_sp, 0.5_sp) + do i = 1, k + res(i) = nor_rvs(loc, scale) ! 2 dummies + end do + res(6:10) = nor_rvs(loc, scale, k) ! 3 dummies + call check(all(abs(real(res) - real(ans)) < sptol) .and. & + all(abs(aimag(res) - aimag(ans)) < sptol), & + msg="normal_distribution_rvs_csp failed", warn=warn) + end subroutine test_nor_rvs_csp + + subroutine test_nor_rvs_cdp + complex(dp) :: res(10), loc, scale + integer :: i, n, k = 5 + integer :: seed, get + complex(dp) :: ans(10) = [(2.12531018257141113281250000000000000_dp, & + 1.46507704257965087890625000000000000_dp), & + (1.08284163475036621093750000000000000_dp, & + 0.277168631553649902343750000000000000_dp), & + (1.41924941539764404296875000000000000_dp, & + 0.498445570468902587890625000000000000_dp), & + (1.72639131546020507812500000000000000_dp, & + 0.715802907943725585937500000000000000_dp), & + (1.98950588703155517578125000000000000_dp, & + 0.115721315145492553710937500000000000_dp), & + (-1.16929018497467041015625000000000000_dp, & + 0.250744730234146118164062500000000000_dp), & + (1.57160544395446777343750000000000000_dp, & + 0.638282597064971923828125000000000000_dp), & + (-1.36106109619140625000000000000000000_dp, & + 0.166259199380874633789062500000000000_dp), & + (1.13403093814849853515625000000000000_dp, & + 1.04232621192932128906250000000000000_dp), & + (-1.68220531940460205078125000000000000_dp, & + 1.63361442089080810546875000000000000_dp)] + + print *, "Test normal_distribution_rvs_cdp" + seed = 25836914 + call random_seed(seed, get) + + loc = (0.5_dp, 1.0_dp); scale = (1.5_dp, 0.5_dp) + do i = 1, k + res(i) = nor_rvs(loc, scale) ! 2 dummies + end do + res(6:10) = nor_rvs(loc, scale, k) ! 3 dummies + call check(all(abs(real(res) - real(ans)) < dptol) .and. & + all(abs(aimag(res) - aimag(ans)) < dptol), & + msg="normal_distribution_rvs_cdp failed", warn=warn) + end subroutine test_nor_rvs_cdp + + subroutine test_nor_rvs_cqp + complex(qp) :: res(10), loc, scale + integer :: i, n, k = 5 + integer :: seed, get + complex(qp) :: ans(10) = [(2.12531018257141113281250000000000000_qp, & + 1.46507704257965087890625000000000000_qp), & + (1.08284163475036621093750000000000000_qp, & + 0.277168631553649902343750000000000000_qp), & + (1.41924941539764404296875000000000000_qp, & + 0.498445570468902587890625000000000000_qp), & + (1.72639131546020507812500000000000000_qp, & + 0.715802907943725585937500000000000000_qp), & + (1.98950588703155517578125000000000000_qp, & + 0.115721315145492553710937500000000000_qp), & + (-1.16929018497467041015625000000000000_qp, & + 0.250744730234146118164062500000000000_qp), & + (1.57160544395446777343750000000000000_qp, & + 0.638282597064971923828125000000000000_qp), & + (-1.36106109619140625000000000000000000_qp, & + 0.166259199380874633789062500000000000_qp), & + (1.13403093814849853515625000000000000_qp, & + 1.04232621192932128906250000000000000_qp), & + (-1.68220531940460205078125000000000000_qp, & + 1.63361442089080810546875000000000000_qp)] + + print *, "Test normal_distribution_rvs_cqp" + seed = 25836914 + call random_seed(seed, get) + + loc = (0.5_qp, 1.0_qp); scale = (1.5_qp, 0.5_qp) + do i = 1, k + res(i) = nor_rvs(loc, scale) ! 2 dummies + end do + res(6:10) = nor_rvs(loc, scale, k) ! 3 dummies + call check(all(abs(real(res) - real(ans)) < qptol) .and. & + all(abs(aimag(res) - aimag(ans)) < qptol), & + msg="normal_distribution_rvs_cqp failed", warn=warn) + end subroutine test_nor_rvs_cqp + + + + subroutine test_nor_pdf_rsp + real(sp) :: x1, x2(3,4), loc, scale + integer :: i, n, k = 5 + integer :: seed, get + real :: res(3,5) + real :: ans(15) = [0.215050772, 0.215050772, 0.215050772, 0.200537622, & + 5.66161536E-02, 0.238986954, 0.265935957,0.262147546,& + 0.249866411, 3.98721099E-02, 0.265902370,0.161311597,& + 0.249177739, 0.237427220, 0.155696079] + + print *, "Test normal_distribution_pdf_rsp" + seed = 741852963 + call random_seed(seed, get) + + loc = -0.5_sp; scale = 1.5_sp + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_pdf(x1, loc, scale) + res(:, 2:5) = nor_pdf(x2, loc, scale) + call check(all(abs(res - reshape(ans, [3,5])) < sptol), & + msg="normal_distribution_pdf_rsp failed", warn=warn) + end subroutine test_nor_pdf_rsp + + subroutine test_nor_pdf_rdp + real(dp) :: x1, x2(3,4), loc, scale + integer :: i, n, k = 5 + integer :: seed, get + real :: res(3,5) + real :: ans(15) = [0.215050772, 0.215050772, 0.215050772, 0.200537622, & + 5.66161536E-02, 0.238986954, 0.265935957,0.262147546,& + 0.249866411, 3.98721099E-02, 0.265902370,0.161311597,& + 0.249177739, 0.237427220, 0.155696079] + + print *, "Test normal_distribution_pdf_rdp" + seed = 741852963 + call random_seed(seed, get) + + loc = -0.5_dp; scale = 1.5_dp + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_pdf(x1, loc, scale) + res(:, 2:5) = nor_pdf(x2, loc, scale) + call check(all(abs(res - reshape(ans, [3,5])) < dptol), & + msg="normal_distribution_pdf_rdp failed", warn=warn) + end subroutine test_nor_pdf_rdp + + subroutine test_nor_pdf_rqp + real(qp) :: x1, x2(3,4), loc, scale + integer :: i, n, k = 5 + integer :: seed, get + real :: res(3,5) + real :: ans(15) = [0.215050772, 0.215050772, 0.215050772, 0.200537622, & + 5.66161536E-02, 0.238986954, 0.265935957,0.262147546,& + 0.249866411, 3.98721099E-02, 0.265902370,0.161311597,& + 0.249177739, 0.237427220, 0.155696079] + + print *, "Test normal_distribution_pdf_rqp" + seed = 741852963 + call random_seed(seed, get) + + loc = -0.5_qp; scale = 1.5_qp + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_pdf(x1, loc, scale) + res(:, 2:5) = nor_pdf(x2, loc, scale) + call check(all(abs(res - reshape(ans, [3,5])) < qptol), & + msg="normal_distribution_pdf_rqp failed", warn=warn) + end subroutine test_nor_pdf_rqp + + subroutine test_nor_pdf_csp + complex(sp) :: x1, x2(3,4), loc, scale + integer :: i, n, k = 5 + integer :: seed, get + real :: res(3,5) + real :: ans(15) = [0.129377306, 0.129377306, 0.129377306,4.05915640E-02,& + 0.209143385, 2.98881028E-02,0.128679410, 0.177484736,& + 3.82205285E-02, 7.09915683E-02, 4.56126593E-02, & + 6.57454208E-02,0.165161029,3.86104845E-02,0.196922958] + + print *, "Test normal_distribution_pdf_csp" + seed = 741852963 + call random_seed(seed, get) + + loc = (-0.5_sp, 0.5_sp); scale = (0.5_sp, 1.5_sp) + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_pdf(x1, loc, scale) + res(:, 2:5) = nor_pdf(x2, loc, scale) + call check(all(abs(res - reshape(ans, [3,5])) < sptol), & + msg="normal_distribution_pdf_csp failed", warn=warn) + end subroutine test_nor_pdf_csp + + subroutine test_nor_pdf_cdp + complex(dp) :: x1, x2(3,4), loc, scale + integer :: i, n, k = 5 + integer :: seed, get + real :: res(3,5) + real :: ans(15) = [0.129377306, 0.129377306, 0.129377306,4.05915640E-02,& + 0.209143385, 2.98881028E-02,0.128679410, 0.177484736,& + 3.82205285E-02, 7.09915683E-02, 4.56126593E-02, & + 6.57454208E-02,0.165161029,3.86104845E-02,0.196922958] + + print *, "Test normal_distribution_pdf_cdp" + seed = 741852963 + call random_seed(seed, get) + + loc = (-0.5_dp, 0.5_dp); scale = (0.5_dp, 1.5_dp) + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_pdf(x1, loc, scale) + res(:, 2:5) = nor_pdf(x2, loc, scale) + call check(all(abs(res - reshape(ans, [3,5])) < dptol), & + msg="normal_distribution_pdf_cdp failed", warn=warn) + end subroutine test_nor_pdf_cdp + + subroutine test_nor_pdf_cqp + complex(qp) :: x1, x2(3,4), loc, scale + integer :: i, n, k = 5 + integer :: seed, get + real :: res(3,5) + real :: ans(15) = [0.129377306, 0.129377306, 0.129377306,4.05915640E-02,& + 0.209143385, 2.98881028E-02,0.128679410, 0.177484736,& + 3.82205285E-02, 7.09915683E-02, 4.56126593E-02, & + 6.57454208E-02,0.165161029,3.86104845E-02,0.196922958] + + print *, "Test normal_distribution_pdf_cqp" + seed = 741852963 + call random_seed(seed, get) + + loc = (-0.5_qp, 0.5_qp); scale = (0.5_qp, 1.5_qp) + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_pdf(x1, loc, scale) + res(:, 2:5) = nor_pdf(x2, loc, scale) + call check(all(abs(res - reshape(ans, [3,5])) < qptol), & + msg="normal_distribution_pdf_cqp failed", warn=warn) + end subroutine test_nor_pdf_cqp + + + subroutine test_nor_cdf_rsp + real(sp) :: x1, x2(3,4), loc, scale + integer :: i, n + integer :: seed, get + real :: res(3,5) + real :: ans(15) = [7.50826299E-02, 7.50826299E-02, 7.50826299E-02, & + 0.143119827, 0.241425425, 0.284345865, 0.233239830, & + 0.341059506,0.353156865,6.81066737E-02,4.38792333E-02,& + 0.763679624, 0.363722175, 0.868187129, 0.626506805] + + print *, "Test normal_distribution_cdf_rsp" + seed = 369147582 + call random_seed(seed, get) + + loc = -1.0_sp; scale = 2.0_sp + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_cdf(x1, loc, scale) + res(:, 2:5) = nor_cdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < sptol), & + msg="normal_distribution_cdf_rsp failed", warn=warn) + end subroutine test_nor_cdf_rsp + + subroutine test_nor_cdf_rdp + real(dp) :: x1, x2(3,4), loc, scale + integer :: i, n + integer :: seed, get + real :: res(3,5) + real :: ans(15) = [7.50826299E-02, 7.50826299E-02, 7.50826299E-02, & + 0.143119827, 0.241425425, 0.284345865, 0.233239830, & + 0.341059506,0.353156865,6.81066737E-02,4.38792333E-02,& + 0.763679624, 0.363722175, 0.868187129, 0.626506805] + + print *, "Test normal_distribution_cdf_rdp" + seed = 369147582 + call random_seed(seed, get) + + loc = -1.0_dp; scale = 2.0_dp + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_cdf(x1, loc, scale) + res(:, 2:5) = nor_cdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < dptol), & + msg="normal_distribution_cdf_rdp failed", warn=warn) + end subroutine test_nor_cdf_rdp + + subroutine test_nor_cdf_rqp + real(qp) :: x1, x2(3,4), loc, scale + integer :: i, n + integer :: seed, get + real :: res(3,5) + real :: ans(15) = [7.50826299E-02, 7.50826299E-02, 7.50826299E-02, & + 0.143119827, 0.241425425, 0.284345865, 0.233239830, & + 0.341059506,0.353156865,6.81066737E-02,4.38792333E-02,& + 0.763679624, 0.363722175, 0.868187129, 0.626506805] + + print *, "Test normal_distribution_cdf_rqp" + seed = 369147582 + call random_seed(seed, get) + + loc = -1.0_qp; scale = 2.0_qp + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_cdf(x1, loc, scale) + res(:, 2:5) = nor_cdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < qptol), & + msg="normal_distribution_cdf_rqp failed", warn=warn) + end subroutine test_nor_cdf_rqp + + subroutine test_nor_cdf_csp + complex(sp) :: x1, x2(3,4), loc, scale + integer :: i, n + integer :: seed, get + real :: res(3,5) + real :: ans(15) = [1.07458141E-02, 1.07458141E-02, 1.07458141E-02, & + 6.86483160E-02, 7.95486569E-02, 2.40523368E-02, & + 3.35096717E-02, 0.315778911,0.446311295,0.102010213, & + 7.66918957E-02, 0.564691007, 0.708769500, & + 6.40553832E-02, 5.39999120E-02] + + print *, "Test normal_distribution_cdf_csp" + seed = 369147582 + call random_seed(seed, get) + + loc = (-1.0_sp, 1.0_sp); scale = (1.7_sp, 2.4_sp) + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_cdf(x1, loc, scale) + res(:, 2:5) = nor_cdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < sptol), & + msg="normal_distribution_cdf_csp failed", warn=warn) + end subroutine test_nor_cdf_csp + + subroutine test_nor_cdf_cdp + complex(dp) :: x1, x2(3,4), loc, scale + integer :: i, n + integer :: seed, get + real :: res(3,5) + real :: ans(15) = [1.07458141E-02, 1.07458141E-02, 1.07458141E-02, & + 6.86483160E-02, 7.95486569E-02, 2.40523368E-02, & + 3.35096717E-02, 0.315778911,0.446311295,0.102010213, & + 7.66918957E-02, 0.564691007, 0.708769500, & + 6.40553832E-02, 5.39999120E-02] + + print *, "Test normal_distribution_cdf_cdp" + seed = 369147582 + call random_seed(seed, get) + + loc = (-1.0_dp, 1.0_dp); scale = (1.7_dp, 2.4_dp) + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_cdf(x1, loc, scale) + res(:, 2:5) = nor_cdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < dptol), & + msg="normal_distribution_cdf_cdp failed", warn=warn) + end subroutine test_nor_cdf_cdp + + subroutine test_nor_cdf_cqp + complex(qp) :: x1, x2(3,4), loc, scale + integer :: i, n + integer :: seed, get + real :: res(3,5) + real :: ans(15) = [1.07458141E-02, 1.07458141E-02, 1.07458141E-02, & + 6.86483160E-02, 7.95486569E-02, 2.40523368E-02, & + 3.35096717E-02, 0.315778911,0.446311295,0.102010213, & + 7.66918957E-02, 0.564691007, 0.708769500, & + 6.40553832E-02, 5.39999120E-02] + + print *, "Test normal_distribution_cdf_cqp" + seed = 369147582 + call random_seed(seed, get) + + loc = (-1.0_qp, 1.0_qp); scale = (1.7_qp, 2.4_qp) + x1 = nor_rvs(loc, scale) + x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) + res(:,1) = nor_cdf(x1, loc, scale) + res(:, 2:5) = nor_cdf(x2, loc, scale) + call check(all(abs(res - reshape(ans,[3,5])) < qptol), & + msg="normal_distribution_cdf_cqp failed", warn=warn) + end subroutine test_nor_cdf_cqp + + +end program test_distribution_normal \ No newline at end of file From 849c278246afdad26411ea053fb6d4e4c73b2299 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sun, 20 Dec 2020 21:17:42 -0500 Subject: [PATCH 06/49] Update CMakeLists.txt --- src/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 33144ad07..b1157149e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -20,7 +20,7 @@ set(fppFiles stdlib_quadrature_simps.fypp stdlib_stats_distribution_PRNG.fypp stdlib_stats_distribution_uniform.fypp - stdlib_stats_distribution_normal.fypp + stdlib_stats_distribution_normal.fypp ) From 448fbdbc1461088c7ba0e8ec7529e92c7a3af259 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sun, 20 Dec 2020 21:18:55 -0500 Subject: [PATCH 07/49] Update Makefile.manual --- src/Makefile.manual | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 725b4aec9..bbab36e42 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -18,7 +18,7 @@ SRC = f18estop.f90 \ stdlib_stats_var.f90 \ stdlib_stats_distribution_PRNG.f90\ stdlib_stats_distribution_uniform.f90 \ - stdlib_stats_distribution_normal.f90 + stdlib_stats_distribution_normal.f90 LIB = libstdlib.a @@ -77,9 +77,9 @@ stdlib_stats_distribution_uniform.o: \ stdlib_stats_distribution_PRNG.o stdlib_stats_distribution_normal.o \ stdlib_kinds.o - stdlib_error.o \ - stdlib_stats_distribution.PRNG.o \ - stdlib_stats_distribution.uniform.o + stdlib_error.o \ + stdlib_stats_distribution.PRNG.o \ + stdlib_stats_distribution.uniform.o # Fortran sources that are built from fypp templates stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp @@ -95,4 +95,4 @@ stdlib_stats_moment.f90: stdlib_stats_moment.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp stdlib_stats_distribution_uniform.f90: stdlib_stats_distribution_uniform.fypp -stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp \ No newline at end of file +stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp From 7f2bb28396447c6b2f8cd5af234be377711f4672 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sun, 20 Dec 2020 21:22:17 -0500 Subject: [PATCH 08/49] Update stdlib_stats_distribution_normal.md --- doc/specs/stdlib_stats_distribution_normal.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index cefbaef87..a32d8d85c 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -71,6 +71,7 @@ program demo_normal_rvs a(:,:,:) = 1.0 b(:,:,:) = 1.0 print *, norm(a,b) ! a rank 3 random variates array + ![0.152776539, -7.51764774E-02, 1.47208166, 0.180561781, 1.32407105, ! 1.20383692, 0.123445868, -0.455737948, -0.469808221, 1.60750175, ! 1.05748117, 0.720934749, 0.407810807, 1.48165631, 2.31749439, @@ -127,7 +128,7 @@ program demo_normal_pdf norm => normal_distribution_rvs implicit none - real :: x(3,4,5),a(3,4,5),b(3,4,5) + real :: x(2,3,4),a(2,3,4),b(2,3,4) complx :: loc, scale integer :: seed_put, seed_get @@ -143,7 +144,7 @@ program demo_normal_pdf !6.47588000E-02 - x = reshape(norm(0.0, 1.0, 60),[3,4,5]) + x = reshape(norm(0.0, 1.0, 24),[2,3,4]) ! standard normal random variates array a(:,:,:) = 0.0 From 2e399eaaf5e4a5711e116f4291e83d958d38402b Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sun, 20 Dec 2020 21:26:07 -0500 Subject: [PATCH 09/49] Update stdlib_stats_distribution_normal.md --- doc/specs/stdlib_stats_distribution_normal.md | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index a32d8d85c..e07c0cbb1 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -59,7 +59,7 @@ program demo_normal_rvs ! 0.563655198 print *, norm(1.0, 2.0) - !normal random variate \(\mu\)=1.0, \(\sigma\)=2.0 + !normal random variate miu=1.0, sigma=2.0 ! -0.633261681 @@ -81,7 +81,7 @@ program demo_normal_rvs loc = (-1.0, 2.0) scale = (2.0, 1.0) print *, norm(loc, scale) - !single complex normal random variate with real part of \(\mu\)=-1, \(\sigma\)=2; imagainary part of \(\mu\)=2.0 and \(\sigma\)=1.0 + !single complex normal random variate with real part of mu=-1, sigma=2; imagainary part of mu=2.0 and sigma=1.0 ! (1.22566295,2.12518454) @@ -140,7 +140,7 @@ program demo_normal_pdf ! 0.241970733 print *, norm_pdf(2.0,-1.0, 2.0) - !a probability density at 2.0 with \(\mu\)=-1.0 \(\sigma\)=2.0 + !a probability density at 2.0 with mu=-1.0 sigma=2.0 !6.47588000E-02 @@ -160,7 +160,7 @@ program demo_normal_pdf loc = (1.0, -0.5) scale = (1.0, 2.) print *, norm_pdf((1.5,1.0), loc, scale) - ! a complex normal probability density function at (1.5,1.0) with real part of \(\mu\)=1.0, \(\sigma\)=1.0 and imaginary part of \(\mu\)=-0.5, \(\sigma\)=2.0 + ! a complex normal probability density function at (1.5,1.0) with real part of mu=1.0, sigma=1.0 and imaginary part of mu=-0.5, sigma=2.0 ! 5.30100204E-02 @@ -219,7 +219,7 @@ program demo_norm_cdf ! 0.841344714 print *, norm_cdf(2.0, -1.0, 2.0) - ! a cumulative at 2.0 with \(\mu\)=-1 \(\sigma\)=2 + ! a cumulative at 2.0 with mu=-1 sigma=2 ! 0.933192849 @@ -239,7 +239,7 @@ program demo_norm_cdf loc = (1.0,0.0) scale = (0.5,1.0) print *, norm_cdf((0.5,-0.5),loc,scale) - !complex normal cumulative distribution at (0.5,-0.5) with real part of \(\mu\)=1.0, \(\sigma\)=0.5 and imaginary part of \(\mu\)=0.0, \(\sigma\)=1.0 + !complex normal cumulative distribution at (0.5,-0.5) with real part of mu=1.0, sigma=0.5 and imaginary part of mu=0.0, sigma=1.0 !4.89511043E-02 From 8da306464bf81cf5fb7daa8d4fc31cb3e669adbe Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Mon, 21 Dec 2020 17:10:25 -0500 Subject: [PATCH 10/49] Update CMakeLists.txt --- src/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b1157149e..35dd035a4 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -20,7 +20,7 @@ set(fppFiles stdlib_quadrature_simps.fypp stdlib_stats_distribution_PRNG.fypp stdlib_stats_distribution_uniform.fypp - stdlib_stats_distribution_normal.fypp + stdlib_stats_distribution_normal.fypp ) From 9d3b2b8995c0aa14b0edb3f3e19707694385abdc Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Mon, 21 Dec 2020 17:13:06 -0500 Subject: [PATCH 11/49] Add files via upload --- CMakeLists.txt | 108 ++++++++++++++++++++---------------------------- Makefile.manual | 108 ++++++++---------------------------------------- 2 files changed, 62 insertions(+), 154 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 33144ad07..6143b5ab5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,70 +1,52 @@ -### Pre-process: .fpp -> .f90 via Fypp - -# Create a list of the files to be preprocessed -set(fppFiles - stdlib_bitsets.fypp - stdlib_bitsets_64.fypp - stdlib_bitsets_large.fypp - stdlib_io.fypp - stdlib_linalg.fypp - stdlib_linalg_diag.fypp - stdlib_optval.fypp - stdlib_stats.fypp - stdlib_stats_corr.fypp - stdlib_stats_cov.fypp - stdlib_stats_mean.fypp - stdlib_stats_moment.fypp - stdlib_stats_var.fypp - stdlib_quadrature.fypp - stdlib_quadrature_trapz.fypp - stdlib_quadrature_simps.fypp - stdlib_stats_distribution_PRNG.fypp - stdlib_stats_distribution_uniform.fypp - stdlib_stats_distribution_normal.fypp -) +cmake_minimum_required(VERSION 3.14.0) +project(stdlib Fortran) +enable_testing() + +include(${CMAKE_SOURCE_DIR}/cmake/stdlib.cmake) + +# --- compiler options +if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU) + add_compile_options(-fimplicit-none) + add_compile_options(-ffree-line-length-132) + add_compile_options(-Wall) + add_compile_options(-Wextra) + add_compile_options(-Wimplicit-procedure) + add_compile_options(-Wconversion-extra) + # -pedantic-errors triggers a false positive for optional arguments of elemental functions, + # see test_optval and https://gcc.gnu.org/bugzilla/show_bug.cgi?id=95446 + add_compile_options(-pedantic-errors) + if(CMAKE_Fortran_COMPILER_VERSION VERSION_GREATER_EQUAL 8.0) + add_compile_options(-std=f2018) + else() + add_compile_options(-std=f2008ts) + endif() +elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Intel) + add_compile_options(-warn declarations,general,usage,interfaces,unused) + add_compile_options(-standard-semantics) + if(CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 18.0) + add_compile_options(-stand f15) + else() + add_compile_options(-stand f18) + endif() +elseif(CMAKE_Fortran_COMPILER_ID STREQUAL PGI) + add_compile_options(-Mdclchk) +endif() +# --- compiler feature checks +include(CheckFortranSourceCompiles) +include(CheckFortranSourceRuns) +check_fortran_source_compiles("error stop i; end" f18errorstop SRC_EXT f90) +check_fortran_source_compiles("real, allocatable :: array(:, :, :, :, :, :, :, :, :, :); end" f03rank SRC_EXT f90) +check_fortran_source_runs("use, intrinsic :: iso_fortran_env, only : real128; real(real128) :: x; x = x+1; end" f03real128) -# Custom preprocessor flags if(DEFINED CMAKE_MAXIMUM_RANK) - set(fyppFlags "-DMAXRANK=${CMAKE_MAXIMUM_RANK}") -elseif(f03rank) - set(fyppFlags) -else() - set(fyppFlags "-DVERSION90") + set(CMAKE_MAXIMUM_RANK ${CMAKE_MAXIMUM_RANK}) endif() -fypp_f90("${fyppFlags}" "${fppFiles}" outFiles) - -set(SRC - stdlib_ascii.f90 - stdlib_error.f90 - stdlib_kinds.f90 - stdlib_logger.f90 - stdlib_system.F90 - ${outFiles} -) - -add_library(fortran_stdlib ${SRC}) - -set(LIB_MOD_DIR ${CMAKE_CURRENT_BINARY_DIR}/mod_files/) -set_target_properties(fortran_stdlib PROPERTIES - Fortran_MODULE_DIRECTORY ${LIB_MOD_DIR}) -target_include_directories(fortran_stdlib PUBLIC - $ - $ -) - -if(f18errorstop) - target_sources(fortran_stdlib PRIVATE f18estop.f90) -else() - target_sources(fortran_stdlib PRIVATE f08estop.f90) +# --- find preprocessor +find_program(FYPP fypp) +if(NOT FYPP) + message(FATAL_ERROR "Preprocessor fypp not found!") endif() -add_subdirectory(tests) - -install(TARGETS fortran_stdlib - RUNTIME DESTINATION bin - ARCHIVE DESTINATION lib - LIBRARY DESTINATION lib - ) -install(DIRECTORY ${LIB_MOD_DIR} DESTINATION include) +add_subdirectory(src) diff --git a/Makefile.manual b/Makefile.manual index 725b4aec9..b8280b102 100644 --- a/Makefile.manual +++ b/Makefile.manual @@ -1,98 +1,24 @@ -SRC = f18estop.f90 \ - stdlib_ascii.f90 \ - stdlib_bitsets.f90 \ - stdlib_bitsets_64.f90 \ - stdlib_bitsets_large.f90 \ - stdlib_error.f90 \ - stdlib_io.f90 \ - stdlib_kinds.f90 \ - stdlib_linalg.f90 \ - stdlib_linalg_diag.f90 \ - stdlib_logger.f90 \ - stdlib_optval.f90 \ - stdlib_quadrature.f90 \ - stdlib_quadrature_trapz.f90 \ - stdlib_stats.f90 \ - stdlib_stats_mean.f90 \ - stdlib_stats_moment.f90 \ - stdlib_stats_var.f90 \ - stdlib_stats_distribution_PRNG.f90\ - stdlib_stats_distribution_uniform.f90 \ - stdlib_stats_distribution_normal.f90 - -LIB = libstdlib.a +# Fortran stdlib Makefile +FC = gfortran +FFLAGS = -Wall -Wextra -Wimplicit-interface -fPIC -g -fcheck=all +FYPPFLAGS= +export FC +export FFLAGS +export FYPPFLAGS -OBJS = $(SRC:.f90=.o) -MODS = $(OBJS:.o=.mod) -SMODS = $(OBJS:.o=*.smod) +.PHONY: all clean test -.PHONY: all clean +all: + $(MAKE) -f Makefile.manual --directory=src + $(MAKE) -f Makefile.manual --directory=src/tests -all: $(LIB) - -$(LIB): $(OBJS) - ar rcs $@ $(OBJS) +test: + $(MAKE) -f Makefile.manual --directory=src/tests test + @echo + @echo "All tests passed." clean: - $(RM) $(LIB) $(OBJS) $(MODS) $(SMODS) - -%.o: %.f90 - $(FC) $(FFLAGS) -c $< - -%.f90: %.fypp - fypp $(FYPPFLAGS) $< $@ - -# Fortran module dependencies -f18estop.o: stdlib_error.o -stdlib_bitsets.o: stdlib_kinds.o -stdlib_bitsets_64.o: stdlib_bitsets.o -stdlib_bitsets_large.o: stdlib_bitsets.o -stdlib_error.o: stdlib_optval.o -stdlib_io.o: \ - stdlib_error.o \ - stdlib_optval.o \ - stdlib_kinds.o -stdlib_linalg_diag.o: stdlib_kinds.o -stdlib_logger.o: stdlib_ascii.o stdlib_optval.o -stdlib_optval.o: stdlib_kinds.o -stdlib_quadrature.o: stdlib_kinds.o -stdlib_stats_mean.o: \ - stdlib_optval.o \ - stdlib_kinds.o \ - stdlib_stats.o -stdlib_stats_moment.o: \ - stdlib_optval.o \ - stdlib_kinds.o \ - stdlib_stats.o -stdlib_stats_var.o: \ - stdlib_optval.o \ - stdlib_kinds.o \ - stdlib_stats.o -stdlib_stats_distribution_PRNG.o: stdlib_kinds.o -stdlib_stats_distribution_uniform.o: \ - stdlib_kinds.o \ - stdlib_error.o \ - stdlib_stats_distribution_PRNG.o -stdlib_stats_distribution_normal.o \ - stdlib_kinds.o - stdlib_error.o \ - stdlib_stats_distribution.PRNG.o \ - stdlib_stats_distribution.uniform.o - -# Fortran sources that are built from fypp templates -stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp -stdlib_bitsets_large.f90: stdlib_bitsets_large.fypp -stdlib_bitsets.f90: stdlib_bitsets.fypp -stdlib_io.f90: stdlib_io.fypp -stdlib_linalg.f90: stdlib_linalg.fypp -stdlib_linalg_diag.f90: stdlib_linalg_diag.fypp -stdlib_quadrature.f90: stdlib_quadrature.fypp -stdlib_stats.f90: stdlib_stats.fypp -stdlib_stats_mean.f90: stdlib_stats_mean.fypp -stdlib_stats_moment.f90: stdlib_stats_moment.fypp -stdlib_stats_var.f90: stdlib_stats_var.fypp -stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp -stdlib_stats_distribution_uniform.f90: stdlib_stats_distribution_uniform.fypp -stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp \ No newline at end of file + $(MAKE) -f Makefile.manual clean --directory=src + $(MAKE) -f Makefile.manual clean --directory=src/tests From 6795c346da1f4a2108d77afcbf86e4f314331121 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Mon, 21 Dec 2020 17:38:55 -0500 Subject: [PATCH 12/49] Add files via upload --- CMakeLists.txt | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 6143b5ab5..d62689913 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,7 +22,9 @@ if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU) endif() elseif(CMAKE_Fortran_COMPILER_ID STREQUAL Intel) add_compile_options(-warn declarations,general,usage,interfaces,unused) - add_compile_options(-standard-semantics) + if(NOT CMAKE_Fortran_COMPILER_VERSION VERSION_EQUAL 20.2.1.20200827) + add_compile_options(-standard-semantics) + endif() if(CMAKE_Fortran_COMPILER_VERSION VERSION_LESS 18.0) add_compile_options(-stand f15) else() @@ -35,7 +37,7 @@ endif() # --- compiler feature checks include(CheckFortranSourceCompiles) include(CheckFortranSourceRuns) -check_fortran_source_compiles("error stop i; end" f18errorstop SRC_EXT f90) +check_fortran_source_runs("i=0; error stop i; end" f18errorstop SRC_EXT f90) check_fortran_source_compiles("real, allocatable :: array(:, :, :, :, :, :, :, :, :, :); end" f03rank SRC_EXT f90) check_fortran_source_runs("use, intrinsic :: iso_fortran_env, only : real128; real(real128) :: x; x = x+1; end" f03real128) From 8f9e681efd99dbda550efe685da64787b968dd53 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Mon, 21 Dec 2020 17:44:30 -0500 Subject: [PATCH 13/49] Update stdlib_stats_distribution_normal.fypp --- src/stdlib_stats_distribution_normal.fypp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index de55412ab..da7abfc2f 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -25,7 +25,7 @@ Module stdlib_stats_distribution_normal !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# !! description)) !! - module procedure norm_dist_rvs_0_rsp !0 dummy varaible + module procedure norm_dist_rvs_0_rsp !0 dummy variable #:for k1, t1 in RC_KINDS_TYPES module procedure norm_dist_rvs_${t1[0]}$${k1}$ !2 dummy variables @@ -363,4 +363,4 @@ Module stdlib_stats_distribution_normal #:endfor -end module stdlib_stats_distribution_normal \ No newline at end of file +end module stdlib_stats_distribution_normal From 80bd8acbfaf666c3f707c0a05b558267d0571adc Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Mon, 21 Dec 2020 19:09:10 -0500 Subject: [PATCH 14/49] Add files via upload From 1a1dfe4fe36528ba37d859551f9fa7c32a9dd5f4 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Mon, 21 Dec 2020 19:15:41 -0500 Subject: [PATCH 15/49] Update Makefile.manual --- src/Makefile.manual | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index bbab36e42..a778aa04c 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -16,7 +16,7 @@ SRC = f18estop.f90 \ stdlib_stats_mean.f90 \ stdlib_stats_moment.f90 \ stdlib_stats_var.f90 \ - stdlib_stats_distribution_PRNG.f90\ + stdlib_stats_distribution_PRNG.f90 \ stdlib_stats_distribution_uniform.f90 \ stdlib_stats_distribution_normal.f90 @@ -72,14 +72,14 @@ stdlib_stats_var.o: \ stdlib_stats.o stdlib_stats_distribution_PRNG.o: stdlib_kinds.o stdlib_stats_distribution_uniform.o: \ - stdlib_kinds.o \ + stdlib_kinds.o \ stdlib_error.o \ stdlib_stats_distribution_PRNG.o stdlib_stats_distribution_normal.o \ - stdlib_kinds.o - stdlib_error.o \ - stdlib_stats_distribution.PRNG.o \ - stdlib_stats_distribution.uniform.o + stdlib_kinds.o + stdlib_error.o \ + stdlib_stats_distribution.PRNG.o \ + stdlib_stats_distribution.uniform.o # Fortran sources that are built from fypp templates stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp From 8457bbe4b1b3b649da1db852f10b57ccee253e60 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Mon, 21 Dec 2020 19:28:36 -0500 Subject: [PATCH 16/49] Update Makefile.manual --- src/Makefile.manual | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index a778aa04c..422d9b028 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -76,7 +76,7 @@ stdlib_stats_distribution_uniform.o: \ stdlib_error.o \ stdlib_stats_distribution_PRNG.o stdlib_stats_distribution_normal.o \ - stdlib_kinds.o + stdlib_kinds.o \ stdlib_error.o \ stdlib_stats_distribution.PRNG.o \ stdlib_stats_distribution.uniform.o From f8881a094b5a6b41fffa4533dd00f17296bd5996 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Mon, 21 Dec 2020 19:41:17 -0500 Subject: [PATCH 17/49] Update Makefile.manual --- src/Makefile.manual | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 422d9b028..06d1e02c5 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -75,7 +75,7 @@ stdlib_stats_distribution_uniform.o: \ stdlib_kinds.o \ stdlib_error.o \ stdlib_stats_distribution_PRNG.o -stdlib_stats_distribution_normal.o \ +stdlib_stats_distribution_normal.o: \ stdlib_kinds.o \ stdlib_error.o \ stdlib_stats_distribution.PRNG.o \ From 256b815c5e7658b0e7a9c20e56a07c5b68945f67 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Mon, 21 Dec 2020 20:09:18 -0500 Subject: [PATCH 18/49] Update Makefile.manual --- src/Makefile.manual | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 06d1e02c5..24a176938 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -78,8 +78,8 @@ stdlib_stats_distribution_uniform.o: \ stdlib_stats_distribution_normal.o: \ stdlib_kinds.o \ stdlib_error.o \ - stdlib_stats_distribution.PRNG.o \ - stdlib_stats_distribution.uniform.o + stdlib_stats_distribution_PRNG.o \ + stdlib_stats_distribution_uniform.o # Fortran sources that are built from fypp templates stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp From 75ee05e781cbb28b6a79625aa5c36bc4a90f181c Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Mon, 21 Dec 2020 20:18:04 -0500 Subject: [PATCH 19/49] Update stdlib_stats_distribution_uniform.fypp --- src/stdlib_stats_distribution_uniform.fypp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/stdlib_stats_distribution_uniform.fypp b/src/stdlib_stats_distribution_uniform.fypp index d092129ce..d51680f0f 100644 --- a/src/stdlib_stats_distribution_uniform.fypp +++ b/src/stdlib_stats_distribution_uniform.fypp @@ -185,7 +185,6 @@ Module stdlib_stats_distribution_uniform ! ${t1}$, intent(in) :: scale ${t1}$ :: res - integer(int64) :: tmp real(${k1}$) :: r1, r2, tr, ti if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & @@ -218,7 +217,6 @@ Module stdlib_stats_distribution_uniform ! ${t1}$, intent(in) :: loc, scale ${t1}$ :: res - integer(int64) :: tmp real(${k1}$) :: r1, r2, tr, ti if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & @@ -481,4 +479,4 @@ Module stdlib_stats_distribution_uniform end function shuffle_${t1[0]}$${k1}$ #:endfor -end module stdlib_stats_distribution_uniform \ No newline at end of file +end module stdlib_stats_distribution_uniform From 86f4f3aab72a396e1ed8e42672feb6d119dd072c Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 22 Dec 2020 10:01:29 -0500 Subject: [PATCH 20/49] Update stdlib_stats_distribution_normal.md --- doc/specs/stdlib_stats_distribution_normal.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index e07c0cbb1..2f80cb685 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -2,7 +2,7 @@ title: stats_distribution --- -# Statistical Distributions Normal Module +# Statistical Distributions -- Normal Distribution Module [TOC] From 374c4b939c54db3b1743961da8e3ae668fbb77f2 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sun, 27 Dec 2020 11:30:50 -0500 Subject: [PATCH 21/49] chn. complex number with kinds --- src/stdlib_stats_distribution_normal.fypp | 36 +++++++++++----------- src/stdlib_stats_distribution_uniform.fypp | 16 +++++----- 2 files changed, 26 insertions(+), 26 deletions(-) diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index da7abfc2f..8397ba08c 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -114,7 +114,7 @@ Module stdlib_stats_distribution_normal if( .not. zig_norm_initialized ) call zigset iz = 0 ! original algorithm use 32bit - hz = dist_rand(iz) + hz = dist_rand(1_int32) iz = iand( hz, 127 ) if( abs( hz ) < kn(iz) ) then @@ -123,8 +123,8 @@ Module stdlib_stats_distribution_normal L1: do L2: if( iz == 0 ) then do - x = -log( uni( ) ) * rr - y = -log( uni( ) ) + x = -log( uni(1.0_${k1}$) ) * rr + y = -log( uni(1.0_${k1}$) ) if( y + y >= x * x ) exit end do res = r + x @@ -132,14 +132,14 @@ Module stdlib_stats_distribution_normal exit L1 end if L2 x = hz * wn(iz) - if( fn(iz) + uni( ) * (fn(iz-1) - fn(iz)) < & + if( fn(iz) + uni(1.0_${k1}$) * (fn(iz-1) - fn(iz)) < & exp(-HALF * x * x) ) then res = x exit L1 end if !original algorithm use 32bit - hz = dist_rand(iz) + hz = dist_rand(1_int32) iz = iand( hz, 127 ) if( abs( hz ) < kn(iz) ) then res = hz * wn(iz) @@ -168,7 +168,7 @@ Module stdlib_stats_distribution_normal if( .not. zig_norm_initialized ) call zigset iz = 0 ! original algorithm use 32bit - hz = dist_rand(iz) + hz = dist_rand(1_int32) iz = iand( hz, 127 ) if( abs( hz ) < kn(iz) ) then @@ -177,8 +177,8 @@ Module stdlib_stats_distribution_normal L1: do L2: if( iz == 0 ) then do - x = -log( uni( ) ) * rr - y = -log( uni( ) ) + x = -log( uni(1.0_${k1}$) ) * rr + y = -log( uni(1.0_${k1}$) ) if( y + y >= x * x ) exit end do res = r + x @@ -186,14 +186,14 @@ Module stdlib_stats_distribution_normal exit L1 end if L2 x = hz * wn(iz) - if( fn(iz) + uni( ) * (fn(iz-1) - fn(iz)) < & + if( fn(iz) + uni(1.0_${k1}$) * (fn(iz-1) - fn(iz)) < & exp(-HALF * x * x) ) then res = x exit L1 end if !original algorithm use 32bit - hz = dist_rand(iz) + hz = dist_rand(1_int32) iz = iand( hz, 127 ) if( abs( hz ) < kn(iz) ) then res = hz * wn(iz) @@ -219,7 +219,7 @@ Module stdlib_stats_distribution_normal tr = norm_dist_rvs_r${k1}$(real(loc), real(scale)) ti = norm_dist_rvs_r${k1}$(aimag(loc), aimag(scale)) - res = cmplx(tr, ti) + res = cmplx(tr, ti, kind=${k1}$) return end function norm_dist_rvs_${t1[0]}$${k1}$ @@ -242,7 +242,7 @@ Module stdlib_stats_distribution_normal do i = 1, array_size iz = 0 ! original algorithm use 32bit - hz = dist_rand(iz) + hz = dist_rand(1_int32) iz = iand( hz, 127 ) if( abs( hz ) < kn(iz) ) then @@ -251,8 +251,8 @@ Module stdlib_stats_distribution_normal L1: do L2: if( iz == 0 ) then do - x = -log( uni( ) ) * rr - y = -log( uni( ) ) + x = -log( uni(1.0_${k1}$) ) * rr + y = -log( uni(1.0_${k1}$) ) if( y + y >= x * x ) exit end do re = r + x @@ -260,14 +260,14 @@ Module stdlib_stats_distribution_normal exit L1 end if L2 x = hz * wn(iz) - if( fn(iz) + uni( ) * (fn(iz-1) - fn(iz)) < & + if( fn(iz) + uni(1.0_${k1}$) * (fn(iz-1) - fn(iz)) < & exp(-HALF * x * x) ) then re = x exit L1 end if !original algorithm use 32bit - hz = dist_rand(iz) + hz = dist_rand(1_int32) iz = iand( hz, 127 ) if( abs( hz ) < kn(iz) ) then re = hz * wn(iz) @@ -295,7 +295,7 @@ Module stdlib_stats_distribution_normal do i = 1, array_size tr = norm_dist_rvs_r${k1}$(real(loc), real(scale)) ti = norm_dist_rvs_r${k1}$(aimag(loc), aimag(scale)) - res(i) = cmplx(tr, ti) + res(i) = cmplx(tr, ti, kind=${k1}$) end do return end function norm_dist_rvs_array_${t1[0]}$${k1}$ @@ -363,4 +363,4 @@ Module stdlib_stats_distribution_normal #:endfor -end module stdlib_stats_distribution_normal +end module stdlib_stats_distribution_normal \ No newline at end of file diff --git a/src/stdlib_stats_distribution_uniform.fypp b/src/stdlib_stats_distribution_uniform.fypp index d51680f0f..572caaf37 100644 --- a/src/stdlib_stats_distribution_uniform.fypp +++ b/src/stdlib_stats_distribution_uniform.fypp @@ -134,7 +134,7 @@ Module stdlib_stats_distribution_uniform ! Uniformly distributed float in [0,1] ! Based on the paper by Frederic Goualard, "Generating Random Floating- ! Point Numbers By Dividing Integers: a Case Study", Proceedings of - ! ICCS 2020, June 20202, Amsterdam, Netherlands + ! ICCS 2020, June 2020, Amsterdam, Netherlands ! ${t1}$ :: res integer(int64) :: tmp @@ -201,7 +201,7 @@ Module stdlib_stats_distribution_uniform tr = real(scale) * r1 ti = aimag(scale) * r2 endif - res = cmplx(tr, ti) + res = cmplx(tr, ti, kind=${k1}$) return end function unif_dist_rvs_1_${t1[0]}$${k1}$ @@ -233,7 +233,7 @@ Module stdlib_stats_distribution_uniform tr = real(loc) + real(scale) * r1 ti = aimag(loc) + aimag(scale) * r2 endif - res = cmplx(tr, ti) + res = cmplx(tr, ti, kind=${k1}$) return end function unif_dist_rvs_${t1[0]}$${k1}$ @@ -329,7 +329,7 @@ Module stdlib_stats_distribution_uniform tr = real(loc) + real(scale) * r1 ti = aimag(loc) + aimag(scale) * r2 endif - res(i) = cmplx(tr, ti) + res(i) = cmplx(tr, ti, kind=${k1}$) enddo return end function unif_dist_rvs_array_${t1[0]}$${k1}$ @@ -343,10 +343,10 @@ Module stdlib_stats_distribution_uniform if(scale == 0) then res = 0.0 - elseif(x < loc .or. x >loc + scale) then + elseif(x < loc .or. x > (loc + scale)) then res = 0.0 else - res = 1. / (scale + 1) + res = 1. / (scale + 1_${k1}$) end if return end function unif_dist_pdf_${t1[0]}$${k1}$ @@ -400,7 +400,7 @@ Module stdlib_stats_distribution_uniform elseif(x < loc) then res = 0.0 elseif(x >= loc .and. x <= (loc + scale)) then - res = real((x - loc + 1)) / real((scale + 1)) + res = real((x - loc + 1_${k1}$)) / real((scale + 1_${k1}$)) else res = 1.0 end if @@ -479,4 +479,4 @@ Module stdlib_stats_distribution_uniform end function shuffle_${t1[0]}$${k1}$ #:endfor -end module stdlib_stats_distribution_uniform +end module stdlib_stats_distribution_uniform \ No newline at end of file From a605a896dc9061ef2ab8a2556de496fc8467145d Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sun, 27 Dec 2020 11:31:27 -0500 Subject: [PATCH 22/49] chg. complex number with kinds --- src/tests/stats/test_distribution_normal.f90 | 168 +++++++++---------- 1 file changed, 84 insertions(+), 84 deletions(-) diff --git a/src/tests/stats/test_distribution_normal.f90 b/src/tests/stats/test_distribution_normal.f90 index fbb0bdaf7..261aa5e85 100644 --- a/src/tests/stats/test_distribution_normal.f90 +++ b/src/tests/stats/test_distribution_normal.f90 @@ -158,26 +158,26 @@ subroutine test_nor_rvs_csp complex(sp) :: res(10), loc, scale integer :: i, n, k = 5 integer :: seed, get - complex(sp) :: ans(10) = [(2.12531018257141113281250000000000000_sp, & - 1.46507704257965087890625000000000000_sp), & - (1.08284163475036621093750000000000000_sp, & - 0.277168631553649902343750000000000000_sp), & - (1.41924941539764404296875000000000000_sp, & - 0.498445570468902587890625000000000000_sp), & - (1.72639131546020507812500000000000000_sp, & - 0.715802907943725585937500000000000000_sp), & - (1.98950588703155517578125000000000000_sp, & - 0.115721315145492553710937500000000000_sp), & - (-1.16929018497467041015625000000000000_sp, & - 0.250744730234146118164062500000000000_sp), & - (1.57160544395446777343750000000000000_sp, & - 0.638282597064971923828125000000000000_sp), & - (-1.36106109619140625000000000000000000_sp, & - 0.166259199380874633789062500000000000_sp), & - (1.13403093814849853515625000000000000_sp, & - 1.04232621192932128906250000000000000_sp), & - (-1.68220531940460205078125000000000000_sp, & - 1.63361442089080810546875000000000000_sp)] + complex(sp) :: ans(10) = [(2.12531029488530509574673033057479188_sp, & + 1.46507698734032082432676702410390135_sp), & + (1.08284164094813181722365413861552952_sp, & + 0.277168639672963013076412153168348595_sp), & + (1.41924946329521489696290359461272601_sp, & + 0.498445561155580918466512230224907398_sp), & + (1.72639126368764062036120776610914618_sp, & + 0.715802936564464420410303091557580046_sp), & + (1.98950590834134349860207180427096318_sp, & + 0.115721315405046931701349421928171068_sp), & + (-1.16929014824793620075382705181255005_sp, & + 0.250744737486995217246033007540972903_sp), & + (1.57160542831869509683428987045772374_sp, & + 0.638282596371312238581197107123443857_sp), & + (-1.36106107654239116833139178197598085_sp, & + 0.166259201494369124318950525776017457_sp), & + (1.13403096805387920698038328737311531_sp, & + 1.04232618148691447146347854868508875_sp), & + (-1.68220535920475811053620418533682823_sp, & + 1.63361446685040256898702182297711261_sp)] print *, "Test normal_distribution_rvs_csp" seed = 25836914 @@ -197,26 +197,26 @@ subroutine test_nor_rvs_cdp complex(dp) :: res(10), loc, scale integer :: i, n, k = 5 integer :: seed, get - complex(dp) :: ans(10) = [(2.12531018257141113281250000000000000_dp, & - 1.46507704257965087890625000000000000_dp), & - (1.08284163475036621093750000000000000_dp, & - 0.277168631553649902343750000000000000_dp), & - (1.41924941539764404296875000000000000_dp, & - 0.498445570468902587890625000000000000_dp), & - (1.72639131546020507812500000000000000_dp, & - 0.715802907943725585937500000000000000_dp), & - (1.98950588703155517578125000000000000_dp, & - 0.115721315145492553710937500000000000_dp), & - (-1.16929018497467041015625000000000000_dp, & - 0.250744730234146118164062500000000000_dp), & - (1.57160544395446777343750000000000000_dp, & - 0.638282597064971923828125000000000000_dp), & - (-1.36106109619140625000000000000000000_dp, & - 0.166259199380874633789062500000000000_dp), & - (1.13403093814849853515625000000000000_dp, & - 1.04232621192932128906250000000000000_dp), & - (-1.68220531940460205078125000000000000_dp, & - 1.63361442089080810546875000000000000_dp)] + complex(dp) :: ans(10) = [(2.12531029488530509574673033057479188_dp, & + 1.46507698734032082432676702410390135_dp), & + (1.08284164094813181722365413861552952_dp, & + 0.277168639672963013076412153168348595_dp), & + (1.41924946329521489696290359461272601_dp, & + 0.498445561155580918466512230224907398_dp), & + (1.72639126368764062036120776610914618_dp, & + 0.715802936564464420410303091557580046_dp), & + (1.98950590834134349860207180427096318_dp, & + 0.115721315405046931701349421928171068_dp), & + (-1.16929014824793620075382705181255005_dp, & + 0.250744737486995217246033007540972903_dp), & + (1.57160542831869509683428987045772374_dp, & + 0.638282596371312238581197107123443857_dp), & + (-1.36106107654239116833139178197598085_dp, & + 0.166259201494369124318950525776017457_dp), & + (1.13403096805387920698038328737311531_dp, & + 1.04232618148691447146347854868508875_dp), & + (-1.68220535920475811053620418533682823_dp, & + 1.63361446685040256898702182297711261_dp)] print *, "Test normal_distribution_rvs_cdp" seed = 25836914 @@ -236,26 +236,26 @@ subroutine test_nor_rvs_cqp complex(qp) :: res(10), loc, scale integer :: i, n, k = 5 integer :: seed, get - complex(qp) :: ans(10) = [(2.12531018257141113281250000000000000_qp, & - 1.46507704257965087890625000000000000_qp), & - (1.08284163475036621093750000000000000_qp, & - 0.277168631553649902343750000000000000_qp), & - (1.41924941539764404296875000000000000_qp, & - 0.498445570468902587890625000000000000_qp), & - (1.72639131546020507812500000000000000_qp, & - 0.715802907943725585937500000000000000_qp), & - (1.98950588703155517578125000000000000_qp, & - 0.115721315145492553710937500000000000_qp), & - (-1.16929018497467041015625000000000000_qp, & - 0.250744730234146118164062500000000000_qp), & - (1.57160544395446777343750000000000000_qp, & - 0.638282597064971923828125000000000000_qp), & - (-1.36106109619140625000000000000000000_qp, & - 0.166259199380874633789062500000000000_qp), & - (1.13403093814849853515625000000000000_qp, & - 1.04232621192932128906250000000000000_qp), & - (-1.68220531940460205078125000000000000_qp, & - 1.63361442089080810546875000000000000_qp)] + complex(qp) :: ans(10) = [(2.12531029488530509574673033057479188_qp, & + 1.46507698734032082432676702410390135_qp), & + (1.08284164094813181722365413861552952_qp, & + 0.277168639672963013076412153168348595_qp), & + (1.41924946329521489696290359461272601_qp, & + 0.498445561155580918466512230224907398_qp), & + (1.72639126368764062036120776610914618_qp, & + 0.715802936564464420410303091557580046_qp), & + (1.98950590834134349860207180427096318_qp, & + 0.115721315405046931701349421928171068_qp), & + (-1.16929014824793620075382705181255005_qp, & + 0.250744737486995217246033007540972903_qp), & + (1.57160542831869509683428987045772374_qp, & + 0.638282596371312238581197107123443857_qp), & + (-1.36106107654239116833139178197598085_qp, & + 0.166259201494369124318950525776017457_qp), & + (1.13403096805387920698038328737311531_qp, & + 1.04232618148691447146347854868508875_qp), & + (-1.68220535920475811053620418533682823_qp, & + 1.63361446685040256898702182297711261_qp)] print *, "Test normal_distribution_rvs_cqp" seed = 25836914 @@ -347,10 +347,10 @@ subroutine test_nor_pdf_csp integer :: i, n, k = 5 integer :: seed, get real :: res(3,5) - real :: ans(15) = [0.129377306, 0.129377306, 0.129377306,4.05915640E-02,& - 0.209143385, 2.98881028E-02,0.128679410, 0.177484736,& - 3.82205285E-02, 7.09915683E-02, 4.56126593E-02, & - 6.57454208E-02,0.165161029,3.86104845E-02,0.196922958] + real :: ans(15) = [0.129377320, 0.129377320,0.129377320,4.05915640E-02, & + 0.209143385,2.98881028E-02, 0.128679410, 0.177484736,& + 3.82205322E-02, 7.09915683E-02, 4.56126593E-02, & + 6.57454133E-02,0.165161043,3.86104807E-02,0.196922958] print *, "Test normal_distribution_pdf_csp" seed = 741852963 @@ -370,10 +370,10 @@ subroutine test_nor_pdf_cdp integer :: i, n, k = 5 integer :: seed, get real :: res(3,5) - real :: ans(15) = [0.129377306, 0.129377306, 0.129377306,4.05915640E-02,& - 0.209143385, 2.98881028E-02,0.128679410, 0.177484736,& - 3.82205285E-02, 7.09915683E-02, 4.56126593E-02, & - 6.57454208E-02,0.165161029,3.86104845E-02,0.196922958] + real :: ans(15) = [0.129377320, 0.129377320,0.129377320,4.05915640E-02, & + 0.209143385,2.98881028E-02, 0.128679410, 0.177484736,& + 3.82205322E-02, 7.09915683E-02, 4.56126593E-02, & + 6.57454133E-02,0.165161043,3.86104807E-02,0.196922958] print *, "Test normal_distribution_pdf_cdp" seed = 741852963 @@ -393,10 +393,10 @@ subroutine test_nor_pdf_cqp integer :: i, n, k = 5 integer :: seed, get real :: res(3,5) - real :: ans(15) = [0.129377306, 0.129377306, 0.129377306,4.05915640E-02,& - 0.209143385, 2.98881028E-02,0.128679410, 0.177484736,& - 3.82205285E-02, 7.09915683E-02, 4.56126593E-02, & - 6.57454208E-02,0.165161029,3.86104845E-02,0.196922958] + real :: ans(15) = [0.129377320, 0.129377320,0.129377320,4.05915640E-02, & + 0.209143385,2.98881028E-02, 0.128679410, 0.177484736,& + 3.82205322E-02, 7.09915683E-02, 4.56126593E-02, & + 6.57454133E-02,0.165161043,3.86104807E-02,0.196922958] print *, "Test normal_distribution_pdf_cqp" seed = 741852963 @@ -486,11 +486,11 @@ subroutine test_nor_cdf_csp integer :: i, n integer :: seed, get real :: res(3,5) - real :: ans(15) = [1.07458141E-02, 1.07458141E-02, 1.07458141E-02, & - 6.86483160E-02, 7.95486569E-02, 2.40523368E-02, & - 3.35096717E-02, 0.315778911,0.446311295,0.102010213, & + real :: ans(15) = [1.07458131E-02, 1.07458131E-02, 1.07458131E-02, & + 6.86483234E-02, 7.95486644E-02, 2.40523387E-02, & + 3.35096754E-02,0.315778911,0.446311295, 0.102010213, & 7.66918957E-02, 0.564691007, 0.708769500, & - 6.40553832E-02, 5.39999120E-02] + 6.40553832E-02, 5.39999157E-02] print *, "Test normal_distribution_cdf_csp" seed = 369147582 @@ -510,11 +510,11 @@ subroutine test_nor_cdf_cdp integer :: i, n integer :: seed, get real :: res(3,5) - real :: ans(15) = [1.07458141E-02, 1.07458141E-02, 1.07458141E-02, & - 6.86483160E-02, 7.95486569E-02, 2.40523368E-02, & - 3.35096717E-02, 0.315778911,0.446311295,0.102010213, & + real :: ans(15) = [1.07458131E-02, 1.07458131E-02, 1.07458131E-02, & + 6.86483234E-02, 7.95486644E-02, 2.40523387E-02, & + 3.35096754E-02,0.315778911,0.446311295, 0.102010213, & 7.66918957E-02, 0.564691007, 0.708769500, & - 6.40553832E-02, 5.39999120E-02] + 6.40553832E-02, 5.39999157E-02] print *, "Test normal_distribution_cdf_cdp" seed = 369147582 @@ -534,11 +534,11 @@ subroutine test_nor_cdf_cqp integer :: i, n integer :: seed, get real :: res(3,5) - real :: ans(15) = [1.07458141E-02, 1.07458141E-02, 1.07458141E-02, & - 6.86483160E-02, 7.95486569E-02, 2.40523368E-02, & - 3.35096717E-02, 0.315778911,0.446311295,0.102010213, & + real :: ans(15) = [1.07458131E-02, 1.07458131E-02, 1.07458131E-02, & + 6.86483234E-02, 7.95486644E-02, 2.40523387E-02, & + 3.35096754E-02,0.315778911,0.446311295, 0.102010213, & 7.66918957E-02, 0.564691007, 0.708769500, & - 6.40553832E-02, 5.39999120E-02] + 6.40553832E-02, 5.39999157E-02] print *, "Test normal_distribution_cdf_cqp" seed = 369147582 From 9e0f2b9b2017bcda9af059d8a6ad226f69421ec8 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:43:44 -0500 Subject: [PATCH 23/49] Update Makefile.manual --- src/Makefile.manual | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 24a176938..4d57b7417 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -14,7 +14,9 @@ SRC = f18estop.f90 \ stdlib_quadrature_trapz.f90 \ stdlib_stats.f90 \ stdlib_stats_mean.f90 \ - stdlib_stats_moment.f90 \ + stdlib_stats_moment_all.f90 \ + stdlib_stats_moment_mask.f90 \ + stdlib_stats_moment_scalar.f90 \ stdlib_stats_var.f90 \ stdlib_stats_distribution_PRNG.f90 \ stdlib_stats_distribution_uniform.f90 \ From e662e6409c655fd96fcdbe399c91366c17b0de81 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Tue, 29 Dec 2020 17:47:33 -0500 Subject: [PATCH 24/49] Update Makefile.manual --- src/Makefile.manual | 1 + 1 file changed, 1 insertion(+) diff --git a/src/Makefile.manual b/src/Makefile.manual index 4d57b7417..ea6f587a0 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -14,6 +14,7 @@ SRC = f18estop.f90 \ stdlib_quadrature_trapz.f90 \ stdlib_stats.f90 \ stdlib_stats_mean.f90 \ + stdlib_stats_moment.f90 \ stdlib_stats_moment_all.f90 \ stdlib_stats_moment_mask.f90 \ stdlib_stats_moment_scalar.f90 \ From d83fe37919d055274f8601626f2bd319095f825b Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 31 Dec 2020 15:10:04 -0500 Subject: [PATCH 25/49] Update CMakeLists.txt --- src/CMakeLists.txt | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 35dd035a4..f17389d56 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -18,9 +18,7 @@ set(fppFiles stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp stdlib_quadrature_simps.fypp - stdlib_stats_distribution_PRNG.fypp - stdlib_stats_distribution_uniform.fypp - stdlib_stats_distribution_normal.fypp + ) From d47d5327d261de80585366658ff2e355621dc2b9 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 31 Dec 2020 15:10:58 -0500 Subject: [PATCH 26/49] Update Makefile.manual --- src/Makefile.manual | 22 ++++------------------ 1 file changed, 4 insertions(+), 18 deletions(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index ea6f587a0..5a855db0d 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -18,11 +18,8 @@ SRC = f18estop.f90 \ stdlib_stats_moment_all.f90 \ stdlib_stats_moment_mask.f90 \ stdlib_stats_moment_scalar.f90 \ - stdlib_stats_var.f90 \ - stdlib_stats_distribution_PRNG.f90 \ - stdlib_stats_distribution_uniform.f90 \ - stdlib_stats_distribution_normal.f90 - + stdlib_stats_var.f90 + LIB = libstdlib.a @@ -73,16 +70,7 @@ stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o -stdlib_stats_distribution_PRNG.o: stdlib_kinds.o -stdlib_stats_distribution_uniform.o: \ - stdlib_kinds.o \ - stdlib_error.o \ - stdlib_stats_distribution_PRNG.o -stdlib_stats_distribution_normal.o: \ - stdlib_kinds.o \ - stdlib_error.o \ - stdlib_stats_distribution_PRNG.o \ - stdlib_stats_distribution_uniform.o + # Fortran sources that are built from fypp templates stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp @@ -96,6 +84,4 @@ stdlib_stats.f90: stdlib_stats.fypp stdlib_stats_mean.f90: stdlib_stats_mean.fypp stdlib_stats_moment.f90: stdlib_stats_moment.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp -stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp -stdlib_stats_distribution_uniform.f90: stdlib_stats_distribution_uniform.fypp -stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp + From 7d24bf643503f3366a95ddd8148f9de0964b711f Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 31 Dec 2020 15:11:47 -0500 Subject: [PATCH 27/49] Update Makefile.manual --- src/Makefile.manual | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 5a855db0d..8b54c2144 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -71,7 +71,6 @@ stdlib_stats_var.o: \ stdlib_kinds.o \ stdlib_stats.o - # Fortran sources that are built from fypp templates stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp stdlib_bitsets_large.f90: stdlib_bitsets_large.fypp @@ -84,4 +83,3 @@ stdlib_stats.f90: stdlib_stats.fypp stdlib_stats_mean.f90: stdlib_stats_mean.fypp stdlib_stats_moment.f90: stdlib_stats_moment.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp - From 9b2c7c16c369c6f88caafbd0997db86c9f45ab49 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 31 Dec 2020 15:15:36 -0500 Subject: [PATCH 28/49] Update Makefile.manual --- src/Makefile.manual | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index dd6f12708..188d85c73 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -18,8 +18,11 @@ SRC = f18estop.f90 \ stdlib_stats_moment_all.f90 \ stdlib_stats_moment_mask.f90 \ stdlib_stats_moment_scalar.f90 \ - stdlib_stats_var.f90 - + stdlib_stats_var.f90 \ + stdlib_stats_distribution_PRNG.f90 \ + stdlib_stats_distribution_uniform.f90 \ + stdlib_stats_distribution_normal.f90 + LIB = libstdlib.a @@ -70,6 +73,16 @@ stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o +stdlib_stats_distribution_PRNG.o: stdlib_kinds.o +stdlib_stats_distribution_uniform.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_PRNG.o +stdlib_stats_distribution_normal.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_PRNG.o \ + stdlib_stats_distribution_uniform.o # Fortran sources that are built from fypp templates stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp @@ -86,3 +99,6 @@ stdlib_stats_moment_all.f90: stdlib_stats_moment_all.fypp stdlib_stats_moment_mask.f90: stdlib_stats_moment_mask.fypp stdlib_stats_moment_scalar.f90: stdlib_stats_moment_scalar.fypp stdlib_stats_var.f90: stdlib_stats_var.fypp +stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp +stdlib_stats_distribution_uniform.f90: stdlib_stats_distribution_uniform.fypp +stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp From 7b8fbc24c3bed83a2c5d4cb27183d7b9717a94e3 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 31 Dec 2020 15:16:40 -0500 Subject: [PATCH 29/49] Update CMakeLists.txt --- src/CMakeLists.txt | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 1704e12ab..3b0836eb7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -21,7 +21,9 @@ set(fppFiles stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp stdlib_quadrature_simps.fypp - + stdlib_stats_distribution_PRNG.fypp + stdlib_stats_distribution_uniform.fypp + stdlib_stats_distribution_normal.fypp ) From be8f7fbd1aaf58e0ac8e56b51f027d01ca505eca Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 14:49:54 -0500 Subject: [PATCH 30/49] Update Makefile.manual --- src/Makefile.manual | 108 ++++++++++++++++++++++---------------------- 1 file changed, 55 insertions(+), 53 deletions(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 188d85c73..9351a374a 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -1,32 +1,35 @@ +SRCFYPP =\ + stdlib_bitsets_64.fypp \ + stdlib_bitsets_large.fypp \ + stdlib_bitsets.fypp \ + stdlib_io.fypp \ + stdlib_linalg.fypp \ + stdlib_linalg_diag.fypp \ + stdlib_optval.fypp \ + stdlib_quadrature.fypp \ + stdlib_quadrature_trapz.fypp \ + stdlib_quadrature_simps.fypp \ + stdlib_stats.fypp \ + stdlib_stats_corr.fypp \ + stdlib_stats_cov.fypp \ + stdlib_stats_mean.fypp \ + stdlib_stats_moment.fypp \ + stdlib_stats_moment_all.fypp \ + stdlib_stats_moment_mask.fypp \ + stdlib_stats_moment_scalar.fypp \ + stdlib_stats_var.fypp + SRC = f18estop.f90 \ stdlib_ascii.f90 \ - stdlib_bitsets.f90 \ - stdlib_bitsets_64.f90 \ - stdlib_bitsets_large.f90 \ stdlib_error.f90 \ - stdlib_io.f90 \ stdlib_kinds.f90 \ - stdlib_linalg.f90 \ - stdlib_linalg_diag.f90 \ stdlib_logger.f90 \ - stdlib_optval.f90 \ - stdlib_quadrature.f90 \ - stdlib_quadrature_trapz.f90 \ - stdlib_stats.f90 \ - stdlib_stats_mean.f90 \ - stdlib_stats_moment.f90 \ - stdlib_stats_moment_all.f90 \ - stdlib_stats_moment_mask.f90 \ - stdlib_stats_moment_scalar.f90 \ - stdlib_stats_var.f90 \ - stdlib_stats_distribution_PRNG.f90 \ - stdlib_stats_distribution_uniform.f90 \ - stdlib_stats_distribution_normal.f90 - -LIB = libstdlib.a + $(SRCGEN) +LIB = libstdlib.a +SRCGEN = $(SRCFYPP:.fypp=.f90) OBJS = $(SRC:.f90=.o) MODS = $(OBJS:.o=.mod) SMODS = $(OBJS:.o=*.smod) @@ -39,12 +42,12 @@ $(LIB): $(OBJS) ar rcs $@ $(OBJS) clean: - $(RM) $(LIB) $(OBJS) $(MODS) $(SMODS) + $(RM) $(LIB) $(OBJS) $(MODS) $(SMODS) $(SRCGEN) %.o: %.f90 $(FC) $(FFLAGS) -c $< -%.f90: %.fypp +$(SRCGEN): %.f90: %.fypp common.fypp fypp $(FYPPFLAGS) $< $@ # Fortran module dependencies @@ -57,10 +60,32 @@ stdlib_io.o: \ stdlib_error.o \ stdlib_optval.o \ stdlib_kinds.o -stdlib_linalg_diag.o: stdlib_kinds.o +stdlib_linalg.o: \ + stdlib_kinds.o +stdlib_linalg_diag.o: \ + stdlib_linalg.o \ + stdlib_kinds.o stdlib_logger.o: stdlib_ascii.o stdlib_optval.o stdlib_optval.o: stdlib_kinds.o stdlib_quadrature.o: stdlib_kinds.o +stdlib_quadrature_simps.o: \ + stdlib_quadrature.o \ + stdlib_error.o \ + stdlib_kinds.o +stdlib_quadrature_trapz.o: \ + stdlib_quadrature.o \ + stdlib_error.o \ + stdlib_kinds.o +stdlib_stats.o: \ + stdlib_kinds.o +stdlib_stats_corr.o: \ + stdlib_optval.o \ + stdlib_kinds.o \ + stdlib_stats.o +stdlib_stats_cov.o: \ + stdlib_optval.o \ + stdlib_kinds.o \ + stdlib_stats.o stdlib_stats_mean.o: \ stdlib_optval.o \ stdlib_kinds.o \ @@ -69,36 +94,13 @@ stdlib_stats_moment.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o +stdlib_stats_moment_all.o: \ + stdlib_stats_moment.o +stdlib_stats_moment_mask.o: \ + stdlib_stats_moment.o +stdlib_stats_moment_scalar.o: \ + stdlib_stats_moment.o stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o -stdlib_stats_distribution_PRNG.o: stdlib_kinds.o -stdlib_stats_distribution_uniform.o: \ - stdlib_kinds.o \ - stdlib_error.o \ - stdlib_stats_distribution_PRNG.o -stdlib_stats_distribution_normal.o: \ - stdlib_kinds.o \ - stdlib_error.o \ - stdlib_stats_distribution_PRNG.o \ - stdlib_stats_distribution_uniform.o - -# Fortran sources that are built from fypp templates -stdlib_bitsets_64.f90: stdlib_bitsets_64.fypp -stdlib_bitsets_large.f90: stdlib_bitsets_large.fypp -stdlib_bitsets.f90: stdlib_bitsets.fypp -stdlib_io.f90: stdlib_io.fypp -stdlib_linalg.f90: stdlib_linalg.fypp -stdlib_linalg_diag.f90: stdlib_linalg_diag.fypp -stdlib_quadrature.f90: stdlib_quadrature.fypp -stdlib_stats.f90: stdlib_stats.fypp -stdlib_stats_mean.f90: stdlib_stats_mean.fypp -stdlib_stats_moment.f90: stdlib_stats_moment.fypp -stdlib_stats_moment_all.f90: stdlib_stats_moment_all.fypp -stdlib_stats_moment_mask.f90: stdlib_stats_moment_mask.fypp -stdlib_stats_moment_scalar.f90: stdlib_stats_moment_scalar.fypp -stdlib_stats_var.f90: stdlib_stats_var.fypp -stdlib_stats_distribution_PRNG.f90: stdlib_stats_distribution_PRNG.fypp -stdlib_stats_distribution_uniform.f90: stdlib_stats_distribution_uniform.fypp -stdlib_stats_distribution_normal.f90: stdlib_stats_distribution_normal.fypp From 0c548896062f99b1ac0b3a01e835ae8c4898640f Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 16:40:12 -0500 Subject: [PATCH 31/49] Update Makefile.manual --- src/Makefile.manual | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 9351a374a..63f43a5b7 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -17,7 +17,10 @@ SRCFYPP =\ stdlib_stats_moment_all.fypp \ stdlib_stats_moment_mask.fypp \ stdlib_stats_moment_scalar.fypp \ - stdlib_stats_var.fypp + stdlib_stats_var.fypp \ + stdlib_stats_distribution_PRNG.fypp \ + stdlib_stats_distribution_uniform.fypp \ + stdlib_stats_distribution_normal.fypp SRC = f18estop.f90 \ stdlib_ascii.f90 \ @@ -104,3 +107,15 @@ stdlib_stats_var.o: \ stdlib_optval.o \ stdlib_kinds.o \ stdlib_stats.o +stdlib_stats_distribution_PRNG.o:\ + stdlib_kinds.o \ + stdlib_error.o +stdlib_stats_distribution_uniform.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_PRNG.o +stdlib_stats_distribution_normal.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_stats_distribution_PRNG.o + stdlib_stats_distribution_uniform.o From da63c88b67df19e03cd82263ef40ed7d6c6ecbe2 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 16:41:00 -0500 Subject: [PATCH 32/49] Add files via upload --- src/stdlib_stats_distribution_PRNG.fypp | 84 ++++------------------ src/stdlib_stats_distribution_normal.fypp | 16 ++--- src/stdlib_stats_distribution_uniform.fypp | 46 ++++++------ 3 files changed, 45 insertions(+), 101 deletions(-) diff --git a/src/stdlib_stats_distribution_PRNG.fypp b/src/stdlib_stats_distribution_PRNG.fypp index 90d9f367e..3fdbf0438 100644 --- a/src/stdlib_stats_distribution_PRNG.fypp +++ b/src/stdlib_stats_distribution_PRNG.fypp @@ -1,16 +1,16 @@ #:include "common.fypp" module stdlib_stats_distribution_PRNG - use stdlib_kinds + use stdlib_kinds, only: int8, int16, int32, int64 + use stdlib_error implicit none private - integer, parameter :: MAX_INT_BIT_SIZE = 64 - integer(int64), save :: st(4), si = 614872703977525537_int64 + integer, parameter :: MAX_INT_BIT_SIZE = bit_size(1_int64) + integer(int64), save :: st(4) ! internal states for xoshiro256ss function + integer(int64), save :: si = 614872703977525537_int64 ! default seed value logical, save :: seed_initialized = .false. public :: random_seed public :: dist_rand - public :: jump - public :: long_jump interface dist_rand @@ -43,13 +43,16 @@ module stdlib_stats_distribution_PRNG function dist_rand_${t1[0]}$${k1}$(n) result(res) !! Random integer generation for various kinds !! result = [-2^k, 2^k - 1], k = 7, 15, 31, 63, depending on input kind - !! Result is used as bit model number instead of regular arithmetic number + !! Result will be operated by bitwise operators to generate desired integer + !! and real pseudorandom numbers !! ${t1}$, intent(in) :: n ${t1}$ :: res integer :: k k = MAX_INT_BIT_SIZE - bit_size(n) + if(k < 0) call error_stop("Error(dist_rand): Integer bit size is" & + //" greater than 64bit") res = shiftr(xoshiro256ss( ), k) end function dist_rand_${t1[0]}$${k1}$ @@ -58,6 +61,7 @@ module stdlib_stats_distribution_PRNG function xoshiro256ss( ) result (res) ! Generate random 64-bit integers ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) + ! http://prng.di.unimi.it/xoshiro256starstar.c ! ! This is xoshiro256** 1.0, one of our all-purpose, rock-solid ! generators. It has excellent (sub-ns) speed, a state (256 bits) that is @@ -81,7 +85,6 @@ module stdlib_stats_distribution_PRNG st(1) = ieor(st(1), st(4)) st(3) = ieor(st(3), t) st(4) = rol64(st(4), 45) - return end function xoshiro256ss function rol64(x, k) result(res) @@ -92,73 +95,9 @@ module stdlib_stats_distribution_PRNG t1 = shiftr(x, (64 - k)) t2 = shiftl(x, k) res = ior(t1, t2) - return end function rol64 - subroutine jump - ! This is the jump function for the xoshiro256ss generator. It is equivalent - ! to 2^128 calls to xoshiro256ss(); it can be used to generate 2^128 - ! non-overlapping subsequences for parallel computations. - ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) - ! - ! Fortran 90 version translated from C by Jim-215-Fisher - integer(int64) :: jp(4) = [1733541517147835066_int64, & - -3051731464161248980_int64, & - -6244198995065845334_int64, & - 4155657270789760540_int64] - integer(int64) :: s1 = 0, s2 = 0, s3 = 0, s4 = 0, c = 1_int64 - integer :: i, j, k - - do i = 1, 4 - do j = 1, 64 - if(iand(jp(i), shiftl(c, j - 1)) /= 0) then - s1 = ieor(s1, st(1)) - s2 = ieor(s2, st(2)) - s3 = ieor(s3, st(3)) - s4 = ieor(s4, st(4)) - end if - k = xoshiro256ss( ) - end do - end do - st(1) = s1 - st(2) = s2 - st(3) = s3 - st(4) = s4 - end subroutine jump - - subroutine long_jump - ! This is the long-jump function for the xoshiro256ss generator. It is - ! equivalent to 2^192 calls to xoshiro256ss(); it can be used to generate - ! 2^64 starting points, from each of which jump() will generate 2^64 - ! non-overlapping subsequences for parallel distributed computations - ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) - ! - ! Fortran 90 version translated from C by Jim-215-Fisher - integer(int64) :: jp(4) = [8566230491382795199_int64, & - -4251311993797857357_int64, & - 8606660816089834049_int64, & - 4111957640723818037_int64] - integer(int64) :: s1 = 0, s2 = 0, s3 = 0, s4 = 0, c = 1_int64 - integer(int32) :: i, j, k - - do i = 1, 4 - do j = 1, 64 - if(iand(jp(i), shiftl(c, j - 1)) /= 0) then - s1 = ieor(s1, st(1)) - s2 = ieor(s2, st(2)) - s3 = ieor(s3, st(3)) - s4 = ieor(s4, st(4)) - end if - k = xoshiro256ss() - end do - end do - st(1) = s1 - st(2) = s2 - st(3) = s3 - st(4) = s4 - end subroutine long_jump - function splitmix64(s) result(res) ! Written in 2015 by Sebastiano Vigna (vigna@acm.org) ! This is a fixed-increment version of Java 8's SplittableRandom @@ -176,6 +115,8 @@ module stdlib_stats_distribution_PRNG data int01, int02, int03/-7046029254386353131_int64, & -4658895280553007687_int64, & -7723592293110705685_int64/ + ! Values are converted from C unsigned integer of 0x9e3779b97f4a7c15, + ! 0xbf58476d1ce4e5b9, 0x94d049bb133111eb if(present(s)) si = s res = si @@ -183,7 +124,6 @@ module stdlib_stats_distribution_PRNG res = ieor(res, shiftr(res, 30)) * int02 res = ieor(res, shiftr(res, 27)) * int03 res = ieor(res, shiftr(res, 31)) - return end function splitmix64 #:for k1, t1 in INT_KINDS_TYPES diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index 8397ba08c..fecef90e8 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -163,8 +163,8 @@ Module stdlib_stats_distribution_normal ${t1}$ :: x, y integer :: hz, iz - if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & - //" parameter must be non-zero") + if(scale==0._${k1}$) call error_stop("Error(norm_dist_rvs): Normal" & + //" distribution scale parameter must be non-zero") if( .not. zig_norm_initialized ) call zigset iz = 0 ! original algorithm use 32bit @@ -235,8 +235,8 @@ Module stdlib_stats_distribution_normal ${t1}$ :: x, y, re integer :: hz, iz, i - if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & - //" parameter must be non-zero") + if(scale==0._${k1}$) call error_stop("Error(norm_dist_rvs_array):" & + //" Normal distribution scale parameter must be non-zero") if( .not. zig_norm_initialized ) call zigset allocate(res(array_size)) do i = 1, array_size @@ -311,8 +311,8 @@ Module stdlib_stats_distribution_normal real :: res ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$)) - if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & - //" parameter must be non-zero") + if(scale==0._${k1}$) call error_stop("Error(norm_dist_pdf):" & + //" Normal distribution scale parameter must be non-zero") res = exp(- 0.5_${k1}$ * (x - loc) * (x - loc) / (scale * scale)) / & (sqrt_2_Pi * scale) return @@ -342,8 +342,8 @@ Module stdlib_stats_distribution_normal real :: res ${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$) - if(scale==0._${k1}$) call error_stop("Error: Normal distribution scale" & - //" parameter must be non-zero") + if(scale==0._${k1}$) call error_stop("Error(norm_dist_cdf):" & + //" Normal distribution scale parameter must be non-zero") res = (1.0_${k1}$ + erf((x - loc) / (scale * sqrt_2))) / 2.0_${k1}$ return end function norm_dist_cdf_${t1[0]}$${k1}$ diff --git a/src/stdlib_stats_distribution_uniform.fypp b/src/stdlib_stats_distribution_uniform.fypp index 572caaf37..b538144f8 100644 --- a/src/stdlib_stats_distribution_uniform.fypp +++ b/src/stdlib_stats_distribution_uniform.fypp @@ -42,7 +42,8 @@ Module stdlib_stats_distribution_uniform interface uniform_distribution_pdf !! Version experiment !! - !! Get uniform distribution probability density (pdf) for integer, real and complex variables + !! Get uniform distribution probability density (pdf) for integer, real and + !! complex variables !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# !! description)) @@ -54,7 +55,8 @@ Module stdlib_stats_distribution_uniform interface uniform_distribution_cdf !! Version experimental !! - !! Get uniform distribution cumulative distribution function (cdf) for integer, real and complex variables + !! Get uniform distribution cumulative distribution function (cdf) for + !! integer, real and complex variables !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# !! description)) !! @@ -66,7 +68,8 @@ Module stdlib_stats_distribution_uniform interface shuffle !! Version experimental !! - !! Fisher-Yates shuffle algorithm for a rank one array of integer, real and complex variables + !! Fisher-Yates shuffle algorithm for a rank one array of integer, real and + !! complex variables !! ([Specification](../page/specs/stdlib_stats_distribution_uniform.html# !! description)) !! @@ -85,13 +88,14 @@ Module stdlib_stats_distribution_uniform ! https://www.pcg-random.org/posts/bounded-rands.html ! ! Fortran 90 translated from c by Jim-215-fisher + ! ${t1}$, intent(in) :: scale ${t1}$ :: res, u, mask, n integer :: zeros, bits_left, bits n = scale - if(n <= 0_${k1}$) call error_stop("Error: Uniform distribution scale" & - //" parameter must be positive") + if(n <= 0_${k1}$) call error_stop("Error(unif_dist_rvs_1): Uniform" & + //" distribution scale parameter must be positive") zeros = leadz(n) bits = bit_size(n) - zeros mask = shiftr(not(0_${k1}$), zeros) @@ -121,8 +125,8 @@ Module stdlib_stats_distribution_uniform ${t1}$, intent(in) :: loc, scale ${t1}$ :: res - if(scale == 0_${k1}$) call error_stop("Error: Uniform distribution" & - //" scale parameter must be non-zero") + if(scale <= 0_${k1}$) call error_stop("Error(unif_dist_rvs): Uniform" & + //" distribution scale parameter must be positive") res = loc + unif_dist_rvs_1_${t1[0]}$${k1}$(scale) return end function unif_dist_rvs_${t1[0]}$${k1}$ @@ -153,8 +157,8 @@ Module stdlib_stats_distribution_uniform ${t1}$, intent(in) :: scale ${t1}$ :: res - if(scale == 0._${k1}$) call error_stop("Error: Uniform distribution" & - //" scale parameter must be non-zero") + if(scale == 0._${k1}$) call error_stop("Error(unif_dist_rvs_1): " & + //"Uniform distribution scale parameter must be non-zero") res = scale * unif_dist_rvs_0_${t1[0]}$${k1}$( ) return end function unif_dist_rvs_1_${t1[0]}$${k1}$ @@ -169,8 +173,8 @@ Module stdlib_stats_distribution_uniform ${t1}$, intent(in) :: loc, scale ${t1}$ :: res - if(scale == 0._${k1}$) call error_stop("Error: Uniform distribution" & - //" scale parameter must be non-zero") + if(scale == 0._${k1}$) call error_stop("Error(unif_dist_rvs): " & + //"Uniform distribution scale parameter must be non-zero") res = loc + scale * unif_dist_rvs_0_${t1[0]}$${k1}$( ) return end function unif_dist_rvs_${t1[0]}$${k1}$ @@ -187,8 +191,8 @@ Module stdlib_stats_distribution_uniform ${t1}$ :: res real(${k1}$) :: r1, r2, tr, ti - if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & - //" distribution scale parameter must be non-zero") + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(uni_dist_" & + //"rvs_1): Uniform distribution scale parameter must be non-zero") r1 = unif_dist_rvs_0_r${k1}$( ) if(real(scale) == 0.0_${k1}$) then ti = aimag(scale) * r1 @@ -219,8 +223,8 @@ Module stdlib_stats_distribution_uniform ${t1}$ :: res real(${k1}$) :: r1, r2, tr, ti - if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & - //" distribution scale parameter must be non-zero") + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(uni_dist_" & + //"rvs): Uniform distribution scale parameter must be non-zero") r1 = unif_dist_rvs_0_r${k1}$( ) if(real(scale) == 0.0_${k1}$) then tr = real(loc) @@ -249,8 +253,8 @@ Module stdlib_stats_distribution_uniform integer :: i, zeros, bits_left, bits n = scale - if(n == 0_${k1}$) call error_stop("Error: Uniform distribution" & - //" scale parameter must be non-zero") + if(n == 0_${k1}$) call error_stop("Error(unif_dist_rvs_array): Uniform" & + //" distribution scale parameter must be non-zero") allocate(res(array_size)) zeros = leadz(n) bits = bit_size(n) - zeros @@ -287,8 +291,8 @@ Module stdlib_stats_distribution_uniform integer :: i - if(scale == 0._${k1}$) call error_stop("Error: Uniform distribution" & - //" scale parameter must be non-zero") + if(scale == 0._${k1}$) call error_stop("Error(unif_dist_rvs_array):" & + //" Uniform distribution scale parameter must be non-zero") allocate(res(array_size)) do i = 1, array_size tmp = shiftr(dist_rand(INT_ONE), 11) @@ -311,8 +315,8 @@ Module stdlib_stats_distribution_uniform integer :: i - if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error: Uniform" & - //" distribution scale parameter must be non-zero") + if(scale == (0.0_${k1}$, 0.0_${k1}$)) call error_stop("Error(unif_dist_"& + //"rvs_array): Uniform distribution scale parameter must be non-zero") allocate(res(array_size)) do i = 1, array_size tmp = shiftr(dist_rand(INT_ONE), 11) From 5158fe19e620c0fc8358fc138f71a896799f3fb7 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 16:41:44 -0500 Subject: [PATCH 33/49] Add files via upload --- doc/specs/stdlib_stats_distribution_normal.md | 88 ++++++++++--------- 1 file changed, 46 insertions(+), 42 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index 2f80cb685..e77626176 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -2,7 +2,7 @@ title: stats_distribution --- -# Statistical Distributions -- Normal Distribution Module +# Statistical Distributions Normal Module [TOC] @@ -30,22 +30,23 @@ With three auguments, the function returns a rank one array of normal distribute `array_size`: optional argument has `intent(in)` and is a scalar of type `integer`. -`loc`: optional argument has `intent(in)` and is a scalar of type `real` or `complx`. +`loc`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. -`scale`: optional argument has `intent(in)` and is a scalar of type `real` or `complx`. +`scale`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. `loc` and `scale` arguments must have the same type. ### Return value -The result is a scalar or rank one array, with a size of `array_size`, of type `real` or `complx`. +The result is a scalar or rank one array, with a size of `array_size`, of type `real` or `complex`. ### Example ```fortran program demo_normal_rvs - use stdlib_stats_distribution_normal, only: norm => normal_distribution_rvs use stdlib_stats_distribution_PRNG, only: random_seed + use stdlib_stats_distribution_normal, only: norm => normal_distribution_rvs + implicit none real :: a(2,3,4), b(2,3,4) complx :: loc, scale @@ -65,23 +66,24 @@ program demo_normal_rvs print *, norm(0.,1.0,10) !an array of 10 standard norml random variates -! [-3.38123664E-02, -0.190365672, 0.220678389, -0.424612164, -0.249541596, -! 0.865260184, 1.11086845, -0.328349441, 1.10873628, 1.27049923] +! -3.38123664E-02 -0.190365672 0.220678389 -0.424612164 -0.249541596 +! 0.865260184 1.11086845 -0.328349441 1.10873628 1.27049923 a(:,:,:) = 1.0 b(:,:,:) = 1.0 print *, norm(a,b) ! a rank 3 random variates array - -![0.152776539, -7.51764774E-02, 1.47208166, 0.180561781, 1.32407105, -! 1.20383692, 0.123445868, -0.455737948, -0.469808221, 1.60750175, -! 1.05748117, 0.720934749, 0.407810807, 1.48165631, 2.31749439, -! 0.414566994, 3.06084275, 1.86505437, 1.36338580, 7.26878643E-02, -! 0.178585172, 1.39557445, 0.828021586, 0.872084975] + +!0.152776539 -7.51764774E-02 1.47208166 0.180561781 1.32407105 +! 1.20383692 0.123445868 -0.455737948 -0.469808221 1.60750175 +! 1.05748117 0.720934749 0.407810807 1.48165631 2.31749439 +! 0.414566994 3.06084275 1.86505437 1.36338580 7.26878643E-02 +! 0.178585172 1.39557445 0.828021586 0.872084975 loc = (-1.0, 2.0) scale = (2.0, 1.0) print *, norm(loc, scale) - !single complex normal random variate with real part of mu=-1, sigma=2; imagainary part of mu=2.0 and sigma=1.0 + !single complex normal random variate with real part of miu=-1, sigma=2; + !imagainary part of miu=2.0 and sigma=1.0 ! (1.22566295,2.12518454) @@ -106,11 +108,11 @@ $$f(x)=\frac{1}{\sigma&space;\sqrt{2&space;\pi}}&space;e^{-\frac{1}{2}(\frac{x-\ ### Arguments -`x`: has `intent(in)` and is a scalar of type `real` or `complx`. +`x`: has `intent(in)` and is a scalar of type `real` or `complex`. -`loc`: has `intent(in)` and is a scalar of type `real` or `complx`. +`loc`: has `intent(in)` and is a scalar of type `real` or `complex`. -`scale`: has `intent(in)` and is a scalar of type `real` or `complx`. +`scale`: has `intent(in)` and is a scalar of type `real` or `complex`. The function is elemental, i.e., all three auguments could be arrays conformable to each other. All three arguments must have the same type. @@ -123,12 +125,12 @@ The result is a scalar or an array, with a shape conformable to auguments, of ty ```fortran program demo_normal_pdf use stdlib_stats_distribution_PRNG, only : random_seed - use stdlib_stats_distribution_normal, only : & - norm_pdf=>normal_distribution_pdf,& - norm => normal_distribution_rvs + use stdlib_stats_distribution_normal, only : & + norm_pdf=>normal_distribution_pdf, & + norm => normal_distribution_rvs implicit none - real :: x(2,3,4),a(2,3,4),b(2,3,4) + real :: x(3,4,5),a(3,4,5),b(3,4,5) complx :: loc, scale integer :: seed_put, seed_get @@ -140,27 +142,28 @@ program demo_normal_pdf ! 0.241970733 print *, norm_pdf(2.0,-1.0, 2.0) - !a probability density at 2.0 with mu=-1.0 sigma=2.0 + !a probability density at 2.0 with miu=-1.0 sigma=2.0 !6.47588000E-02 - x = reshape(norm(0.0, 1.0, 24),[2,3,4]) + x = reshape(norm(0.0, 1.0, 60),[3,4,5]) ! standard normal random variates array a(:,:,:) = 0.0 b(:,:,:) = 1.0 print *, norm_pdf(x, a, b) ! standard normal probability density array -! [0.340346158, 0.285823315, 0.398714304, 0.391778737, 0.389345556, -! 0.364551932, 0.386712372, 0.274370432, 0.215250477, 0.378006011, -! 0.215760440, 0.177990928, 0.278640658, 0.223813817, 0.356875211, -! 0.285167664, 0.378533930, 0.390739858, 0.271684974, 0.138273031, -! 0.135456234, 0.331718773, 0.398283750, 0.383706540] +! 0.340346158 0.285823315 0.398714304 0.391778737 0.389345556 +! 0.364551932 0.386712372 0.274370432 0.215250477 0.378006011 +! 0.215760440 0.177990928 0.278640658 0.223813817 0.356875211 +! 0.285167664 0.378533930 0.390739858 0.271684974 0.138273031 +! 0.135456234 0.331718773 0.398283750 0.383706540 loc = (1.0, -0.5) scale = (1.0, 2.) print *, norm_pdf((1.5,1.0), loc, scale) - ! a complex normal probability density function at (1.5,1.0) with real part of mu=1.0, sigma=1.0 and imaginary part of mu=-0.5, sigma=2.0 + ! a complex normal probability density function at (1.5,1.0) with real part + ! of miu=1.0, sigma=1.0 and imaginary part of miu=-0.5, sigma=2.0 ! 5.30100204E-02 @@ -185,11 +188,11 @@ $$F(X)=\frac{1}{2}\left&space;[&space;1&space;+&space;erf(\frac{x-\mu}{\sqr ### Arguments -`x`: has `intent(in)` and is a scalar of type `real` or `complx`. +`x`: has `intent(in)` and is a scalar of type `real` or `complex`. -`loc`: has `intent(in)` and is a scalar of type `real` or `complx`. +`loc`: has `intent(in)` and is a scalar of type `real` or `complex`. -`scale`: has `intent(in)` and is a scalar of type `real` or `complx`. +`scale`: has `intent(in)` and is a scalar of type `real` or `complex`. The function is elemental, i.e., all three auguments could be arrays conformable to each other. All three arguments must have the same type. @@ -202,9 +205,9 @@ The result is a scalar or an array, with a shape conformable to auguments, of ty ```fortran program demo_norm_cdf use stdlib_stats_distribution_PRNG, only : random_seed - use stdlib_stats_distribution_normal, only : & - norm_cdf => normal_distribution_cdf, & - norm => normal_distribution_rvs + use stdlib_stats_distribution_normal, only : & + norm_cdf => normal_distribution_cdf, & + norm => normal_distribution_rvs implicit none real :: x(2,3,4),a(2,3,4),b(2,3,4) @@ -219,7 +222,7 @@ program demo_norm_cdf ! 0.841344714 print *, norm_cdf(2.0, -1.0, 2.0) - ! a cumulative at 2.0 with mu=-1 sigma=2 + ! a cumulative at 2.0 with miu=-1 sigma=2 ! 0.933192849 @@ -230,16 +233,17 @@ program demo_norm_cdf b(:,:,:) = 1.0 print *, norm_cdf(x, a, b) ! standard normal cumulative array -! [0.713505626, 0.207069695, 0.486513376, 0.424511284, 0.587328553, -! 0.335559726, 0.401470929, 0.806552052, 0.866687536, 0.371323735, -! 0.866228044, 0.898046613, 0.198435277, 0.141147852, 0.681565762, -! 0.206268221, 0.627057910, 0.580759525, 0.190364420, 7.27325380E-02, -! 7.08068311E-02, 0.728241026, 0.522919059, 0.390097380] +! 0.713505626 0.207069695 0.486513376 0.424511284 0.587328553 +! 0.335559726 0.401470929 0.806552052 0.866687536 0.371323735 +! 0.866228044 0.898046613 0.198435277 0.141147852 0.681565762 +! 0.206268221 0.627057910 0.580759525 0.190364420 7.27325380E-02 +! 7.08068311E-02 0.728241026 0.522919059 0.390097380 loc = (1.0,0.0) scale = (0.5,1.0) print *, norm_cdf((0.5,-0.5),loc,scale) - !complex normal cumulative distribution at (0.5,-0.5) with real part of mu=1.0, sigma=0.5 and imaginary part of mu=0.0, sigma=1.0 + !complex normal cumulative distribution at (0.5,-0.5) with real part of + !miu=1.0, sigma=0.5 and imaginary part of miu=0.0, sigma=1.0 !4.89511043E-02 From d6140cd62da0f49c65bb2fe10e7675b398f2c99f Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 20:10:45 -0500 Subject: [PATCH 34/49] Update stdlib_stats_distribution_normal.md --- doc/specs/stdlib_stats_distribution_normal.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index e77626176..8651a1365 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -2,7 +2,7 @@ title: stats_distribution --- -# Statistical Distributions Normal Module +# Statistical Distributions -- Normal Module [TOC] From 7817320585a0e7bb2448ced3c06dc2b719e002a1 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Thu, 21 Jan 2021 22:35:52 -0500 Subject: [PATCH 35/49] Update Makefile.manual --- src/Makefile.manual | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Makefile.manual b/src/Makefile.manual index 63f43a5b7..fefc69509 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -117,5 +117,5 @@ stdlib_stats_distribution_uniform.o: \ stdlib_stats_distribution_normal.o: \ stdlib_kinds.o \ stdlib_error.o \ - stdlib_stats_distribution_PRNG.o + stdlib_stats_distribution_PRNG.o \ stdlib_stats_distribution_uniform.o From e0d6dde28865330ee653f50c5c37fc95e29b697f Mon Sep 17 00:00:00 2001 From: Jim-215-Fisher Date: Wed, 6 Oct 2021 20:34:40 -0400 Subject: [PATCH 36/49] change .f90 to .fypp --- ...test_distribution_normal.f90 => test_distribution_normal.fypp} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/tests/stats/{test_distribution_normal.f90 => test_distribution_normal.fypp} (100%) diff --git a/src/tests/stats/test_distribution_normal.f90 b/src/tests/stats/test_distribution_normal.fypp similarity index 100% rename from src/tests/stats/test_distribution_normal.f90 rename to src/tests/stats/test_distribution_normal.fypp From 7b27e63f149bf1a2a8a3ce0c810122512414f9c3 Mon Sep 17 00:00:00 2001 From: Jim-215-Fisher Date: Thu, 7 Oct 2021 15:32:24 -0400 Subject: [PATCH 37/49] add normal distribution --- doc/specs/index.md | 1 + doc/specs/stdlib_stats_distribution_normal.md | 64 +- src/CMakeLists.txt | 1 + src/Makefile.manual | 6 + src/stdlib_stats_distribution_PRNG.fypp | 151 ----- src/stdlib_stats_distribution_normal.fypp | 243 +++---- src/tests/stats/CMakeLists.txt | 2 + src/tests/stats/Makefile.manual | 3 +- src/tests/stats/test_distribution_normal.fypp | 638 +++++------------- 9 files changed, 312 insertions(+), 797 deletions(-) delete mode 100644 src/stdlib_stats_distribution_PRNG.fypp diff --git a/doc/specs/index.md b/doc/specs/index.md index a3b0a5def..720b232f4 100644 --- a/doc/specs/index.md +++ b/doc/specs/index.md @@ -25,6 +25,7 @@ This is and index/directory of the specifications (specs) for each new module/fe - [sorting](./stdlib_sorting.html) - Sorting of rank one arrays - [stats](./stdlib_stats.html) - Descriptive Statistics - [stats_distributions_uniform](./stdlib_stats_distribution_uniform.html) - Uniform Probability Distribution + - [stats_distributions_normal](./stdlib_stats_distribution_normal.html) - Normal Probability Distribution - [string\_type](./stdlib_string_type.html) - Basic string support - [strings](./stdlib_strings.html) - String handling and manipulation routines - [stringlist_type](./stdlib_stringlist_type.html) - 1-Dimensional list of strings diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index 8651a1365..2dd7d504d 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -1,12 +1,12 @@ --- -title: stats_distribution +title: stats_distribution_normal --- -# Statistical Distributions -- Normal Module +# Statistical Distributions -- Normal Distribution Module [TOC] -## `normal_distribution_rvs` - normal distribution random variates +## `rvs_normal` - normal distribution random variates ### Status @@ -16,15 +16,19 @@ Experimental A normal continuous random variate distribution, also known as Gaussian, or Gauss or Laplace-Gauss distribution. The location `loc` specifies the mean or expectation. The `scale` specifies the standard deviation. -Without augument the function returns a standard normal distributed random variate N(0,1). The function is elemental. +Without argument the function returns a standard normal distributed random variate N(0,1). -With two arguments, the function returns a normal distributed random variate N(loc, scale^2). For complex auguments, the real and imaginary parts are independent of each other. The function is elemental. +With two arguments, the function returns a normal distributed random variate N(loc, scale^2). For complex arguments, the real and imaginary parts are independent of each other. -With three auguments, the function returns a rank one array of normal distributed random variates. +With three arguments, the function returns a rank one array of normal distributed random variates. ### Syntax -`result = [[stdlib_stats_distribution_normal(module):normal_distribution_rvs(interface)]]([loc, scale] [[, array_size]])` +`result = [[stdlib_stats_distribution_normal(module):rvs_normal(interface)]]([loc, scale] [[, array_size]])` + +### Class + +Function ### Arguments @@ -44,8 +48,8 @@ The result is a scalar or rank one array, with a size of `array_size`, of type ` ```fortran program demo_normal_rvs - use stdlib_stats_distribution_PRNG, only: random_seed - use stdlib_stats_distribution_normal, only: norm => normal_distribution_rvs + use stdlib_random, only: random_seed + use stdlib_stats_distribution_normal, only: norm => rvs_normal implicit none real :: a(2,3,4), b(2,3,4) @@ -90,7 +94,7 @@ program demo_normal_rvs end program demo_normal_rvs ``` -## `normal_distribution_pdf` - normal probability density function +## `pdf_normal` - normal distribution probability density function ### Status @@ -100,11 +104,15 @@ Experimental The probability density function of the continuous normal distribution. -$$f(x)=\frac{1}{\sigma&space;\sqrt{2&space;\pi}}&space;e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^{2}})$$ +$$f(x) = \frac{1}{\sigma \sqrt{2}} e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^{2}}$$ ### Syntax -`result = [[stdlib_stats_distribution_normal(module):normal_distribution_pdf(interface)]](x, loc, scale)` +`result = [[stdlib_stats_distribution_normal(module):pdf_normal(interface)]](x, loc, scale)` + +### Class + +Elemental function ### Arguments @@ -114,20 +122,19 @@ $$f(x)=\frac{1}{\sigma&space;\sqrt{2&space;\pi}}&space;e^{-\frac{1}{2}(\frac{x-\ `scale`: has `intent(in)` and is a scalar of type `real` or `complex`. -The function is elemental, i.e., all three auguments could be arrays conformable to each other. All three arguments must have the same type. +All three arguments must have the same type. ### Return value -The result is a scalar or an array, with a shape conformable to auguments, of type `real`. +The result is a scalar or an array, with a shape conformable to arguments, of type `real`. ### Example ```fortran program demo_normal_pdf - use stdlib_stats_distribution_PRNG, only : random_seed - use stdlib_stats_distribution_normal, only : & - norm_pdf=>normal_distribution_pdf, & - norm => normal_distribution_rvs + use stdlib_random, only : random_seed + use stdlib_stats_distribution_normal, only : norm_pdf=>pdf_normal, & + norm => rvs_normal implicit none real :: x(3,4,5),a(3,4,5),b(3,4,5) @@ -170,7 +177,7 @@ program demo_normal_pdf end program demo_normal_pdf ``` -## `normal_distribution_cdf` - normal cumulative distribution function +## `cdf_normal` - normal distribution cumulative distribution function ### Status @@ -180,11 +187,15 @@ Experimental Cumulative distribution function of the normal continuous distribution -$$F(X)=\frac{1}{2}\left&space;[&space;1&space;+&space;erf(\frac{x-\mu}{\sqrt{2}&space;\sigma})&space;\right&space;])$$ +$$F(x)=\frac{1}{2}\left [ 1+erf(\frac{x-\mu}{\sigma \sqrt{2}}) \right ]$$ ### Syntax -`result = [[stdlib_stats_distribution_normal(module):normal_distribution_cdf(interface)]](x, loc, scale)` +`result = [[stdlib_stats_distribution_normal(module):cdf_normal(interface)]](x, loc, scale)` + +### Class + +Elemental function ### Arguments @@ -194,20 +205,19 @@ $$F(X)=\frac{1}{2}\left&space;[&space;1&space;+&space;erf(\frac{x-\mu}{\sqr `scale`: has `intent(in)` and is a scalar of type `real` or `complex`. -The function is elemental, i.e., all three auguments could be arrays conformable to each other. All three arguments must have the same type. +All three arguments must have the same type. ### Return value -The result is a scalar or an array, with a shape conformable to auguments, of type `real`. +The result is a scalar or an array, with a shape conformable to arguments, of type `real`. ### Example ```fortran program demo_norm_cdf - use stdlib_stats_distribution_PRNG, only : random_seed - use stdlib_stats_distribution_normal, only : & - norm_cdf => normal_distribution_cdf, & - norm => normal_distribution_rvs + use stdlib_random, only : random_seed + use stdlib_stats_distribution_normal, only : norm_cdf => cdf_normal, & + norm => rvs_normal implicit none real :: x(2,3,4),a(2,3,4),b(2,3,4) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index bb9fb4fd8..250467ca5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -25,6 +25,7 @@ set(fppFiles stdlib_stats_moment_mask.fypp stdlib_stats_moment_scalar.fypp stdlib_stats_distribution_uniform.fypp + stdlib_stats_distribution_normal.fypp stdlib_stats_var.fypp stdlib_quadrature.fypp stdlib_quadrature_trapz.fypp diff --git a/src/Makefile.manual b/src/Makefile.manual index 179fc600f..ec286ff45 100644 --- a/src/Makefile.manual +++ b/src/Makefile.manual @@ -27,6 +27,7 @@ SRCFYPP = \ stdlib_stats_moment_mask.fypp \ stdlib_stats_moment_scalar.fypp \ stdlib_stats_distribution_uniform.fypp \ + stdlib_stats_distribution_normal.fypp \ stdlib_stats_var.fypp \ stdlib_math.fypp \ stdlib_math_linspace.fypp \ @@ -150,6 +151,11 @@ stdlib_stats_distribution_uniform.o: \ stdlib_kinds.o \ stdlib_error.o \ stdlib_random.o +stdlib_stats_distribution_normal.o: \ + stdlib_kinds.o \ + stdlib_error.o \ + stdlib_random.o \ + stdlib_stats_distribution_uniform.o stdlib_random.o: \ stdlib_kinds.o \ stdlib_error.o diff --git a/src/stdlib_stats_distribution_PRNG.fypp b/src/stdlib_stats_distribution_PRNG.fypp deleted file mode 100644 index 3fdbf0438..000000000 --- a/src/stdlib_stats_distribution_PRNG.fypp +++ /dev/null @@ -1,151 +0,0 @@ -#:include "common.fypp" -module stdlib_stats_distribution_PRNG - use stdlib_kinds, only: int8, int16, int32, int64 - use stdlib_error - implicit none - private - integer, parameter :: MAX_INT_BIT_SIZE = bit_size(1_int64) - integer(int64), save :: st(4) ! internal states for xoshiro256ss function - integer(int64), save :: si = 614872703977525537_int64 ! default seed value - logical, save :: seed_initialized = .false. - - public :: random_seed - public :: dist_rand - - - interface dist_rand - !! Version experimental - !! - !! Generation of random integers with different kinds - !! ([Specification](../page/specs/stdlib_stats_distribution_PRNG.html# - !! description)) - #:for k1, t1 in INT_KINDS_TYPES - module procedure dist_rand_${t1[0]}$${k1}$ - #:endfor - end interface dist_rand - - interface random_seed - !! Version experimental - !! - !! Set seed value for random number generator - !! ([Specification](../page/specs/stdlib_stats_distribution_PRNG.html# - !! description)) - !! - #:for k1, t1 in INT_KINDS_TYPES - module procedure random_distribution_seed_${t1[0]}$${k1}$ - #:endfor - end interface random_seed - - - contains - - #:for k1, t1 in INT_KINDS_TYPES - function dist_rand_${t1[0]}$${k1}$(n) result(res) - !! Random integer generation for various kinds - !! result = [-2^k, 2^k - 1], k = 7, 15, 31, 63, depending on input kind - !! Result will be operated by bitwise operators to generate desired integer - !! and real pseudorandom numbers - !! - ${t1}$, intent(in) :: n - ${t1}$ :: res - integer :: k - - k = MAX_INT_BIT_SIZE - bit_size(n) - if(k < 0) call error_stop("Error(dist_rand): Integer bit size is" & - //" greater than 64bit") - res = shiftr(xoshiro256ss( ), k) - end function dist_rand_${t1[0]}$${k1}$ - - #:endfor - - function xoshiro256ss( ) result (res) - ! Generate random 64-bit integers - ! Written in 2018 by David Blackman and Sebastiano Vigna (vigna@acm.org) - ! http://prng.di.unimi.it/xoshiro256starstar.c - ! - ! This is xoshiro256** 1.0, one of our all-purpose, rock-solid - ! generators. It has excellent (sub-ns) speed, a state (256 bits) that is - ! large enough for any parallel application, and it passes all tests we - ! are aware of. - ! - ! The state must be seeded so that it is not everywhere zero. If you have - ! a 64-bit seed, we suggest to seed a splitmix64 generator and use its - ! output to fill st. - ! - ! Fortran 90 version translated from C by Jim-215-Fisher - ! - integer(int64) :: res, t - - if(.not. seed_initialized) call random_distribution_seed_iint64(si,t) - res = rol64(st(2) * 5 , 7) * 9 - t = shiftl(st(2), 17) - st(3) = ieor(st(3), st(1)) - st(4) = ieor(st(4), st(2)) - st(2) = ieor(st(2), st(3)) - st(1) = ieor(st(1), st(4)) - st(3) = ieor(st(3), t) - st(4) = rol64(st(4), 45) - end function xoshiro256ss - - function rol64(x, k) result(res) - integer(int64), intent(in) :: x - integer, intent(in) :: k - integer(int64) :: t1, t2, res - - t1 = shiftr(x, (64 - k)) - t2 = shiftl(x, k) - res = ior(t1, t2) - end function rol64 - - - function splitmix64(s) result(res) - ! Written in 2015 by Sebastiano Vigna (vigna@acm.org) - ! This is a fixed-increment version of Java 8's SplittableRandom - ! generator. - ! See http://dx.doi.org/10.1145/2714064.2660195 and - ! http://docs.oracle.com/javase/8/docs/api/java/util/SplittableRandom.html - ! - ! It is a very fast generator passing BigCrush, and it can be useful if - ! for some reason you absolutely want 64 bits of state. - ! - ! Fortran 90 translated from C by Jim-215-Fisher - ! - integer(int64) :: res, int01, int02, int03 - integer(int64), intent(in), optional :: s - data int01, int02, int03/-7046029254386353131_int64, & - -4658895280553007687_int64, & - -7723592293110705685_int64/ - ! Values are converted from C unsigned integer of 0x9e3779b97f4a7c15, - ! 0xbf58476d1ce4e5b9, 0x94d049bb133111eb - - if(present(s)) si = s - res = si - si = res + int01 - res = ieor(res, shiftr(res, 30)) * int02 - res = ieor(res, shiftr(res, 27)) * int03 - res = ieor(res, shiftr(res, 31)) - end function splitmix64 - - #:for k1, t1 in INT_KINDS_TYPES - subroutine random_distribution_seed_${t1[0]}$${k1}$(put, get) - !! Set seed value for random number generator - !! - ${t1}$, intent(in) :: put - ${t1}$, intent(out) :: get - integer(int64) :: tmp - integer :: i - - tmp = splitmix64(int(put, kind = int64)) - do i = 1, 10 - tmp = splitmix64( ) - end do - do i = 1, 4 - tmp = splitmix64( ) - st(i) = tmp - end do - get = int(tmp, kind = ${k1}$) - seed_initialized = .true. - end subroutine random_distribution_seed_${t1[0]}$${k1}$ - - #:endfor -end module stdlib_stats_distribution_PRNG \ No newline at end of file diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index fecef90e8..1a5930b95 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -1,10 +1,10 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES -Module stdlib_stats_distribution_normal +module stdlib_stats_distribution_normal use stdlib_kinds use stdlib_error, only : error_stop - use stdlib_stats_distribution_PRNG, only : dist_rand - use stdlib_stats_distribution_uniform, only : uni=>uniform_distribution_rvs + use stdlib_random, only : dist_rand + use stdlib_stats_distribution_uniform, only : uni=>rvs_uniform implicit none private @@ -14,54 +14,63 @@ Module stdlib_stats_distribution_normal real(dp), save :: wn(0:127), fn(0:127) logical, save :: zig_norm_initialized = .false. - public :: normal_distribution_rvs - public :: normal_distribution_pdf - public :: normal_distribution_cdf + public :: rvs_normal + public :: pdf_normal + public :: cdf_normal - interface normal_distribution_rvs + + + interface rvs_normal !! Version experimental !! !! Normal Distribution Random Variates - !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# - !! description)) + !! ([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! rvs_normal-normal-distribution-random-variates)) !! - module procedure norm_dist_rvs_0_rsp !0 dummy variable + module procedure rvs_norm_0_rsp !0 dummy variable #:for k1, t1 in RC_KINDS_TYPES - module procedure norm_dist_rvs_${t1[0]}$${k1}$ !2 dummy variables + module procedure rvs_norm_${t1[0]}$${k1}$ !2 dummy variables #:endfor #:for k1, t1 in RC_KINDS_TYPES - module procedure norm_dist_rvs_array_${t1[0]}$${k1}$ !3 dummy variables + module procedure rvs_norm_array_${t1[0]}$${k1}$ !3 dummy variables #:endfor - end interface normal_distribution_rvs + end interface rvs_normal + - interface normal_distribution_pdf + + interface pdf_normal !! Version experimental !! !! Normal Distribution Probability Density Function - !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# - !! description)) + !! ([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! pdf_normal-normal-distribution-probability-density-function)) !! #:for k1, t1 in RC_KINDS_TYPES - module procedure norm_dist_pdf_${t1[0]}$${k1}$ + module procedure pdf_norm_${t1[0]}$${k1}$ #:endfor - end interface normal_distribution_pdf + end interface pdf_normal + - interface normal_distribution_cdf + + interface cdf_normal !! Version experimental !! !! Normal Distribution Cumulative Distribution Function - !!([Specification](../page/specs/stdlib_stats_distribution_normal.html# - !! description)) + !! ([Specification](../page/specs/stdlib_stats_distribution_normal.html# + !! cdf_normal-normal-distribution-cumulative-distribution-function)) !! #:for k1, t1 in RC_KINDS_TYPES - module procedure norm_dist_cdf_${t1[0]}$${k1}$ + module procedure cdf_norm_${t1[0]}$${k1}$ #:endfor - end interface normal_distribution_cdf + end interface cdf_normal + + + - contains +contains subroutine zigset ! Marsaglia & Tsang generator for random normals & random exponentials. @@ -78,12 +87,12 @@ Module stdlib_stats_distribution_normal ! Latest version - 1 January 2001 ! real(dp), parameter :: M1 = 2147483648.0_dp - real(dp) :: dn = 3.442619855899_dp, tn, & - vn = 0.00991256303526217_dp, q + real(dp) :: dn, tn, q integer :: i + dn = 3.442619855899_dp tn = dn - ! tables for random normals + !tables for random normals q = vn * exp(HALF * dn * dn) kn(0) = int((dn / q) * M1, kind = int32) kn(1) = 0 @@ -99,11 +108,13 @@ Module stdlib_stats_distribution_normal wn(i) = dn / M1 end do zig_norm_initialized = .true. - return end subroutine zigset + + #:for k1, t1 in REAL_KINDS_TYPES - impure elemental function norm_dist_rvs_0_${t1[0]}$${k1}$( ) result(res) + function rvs_norm_0_${t1[0]}$${k1}$( ) result(res) + ! ! Standard normal random vairate (0,1) ! ${t1}$ :: res @@ -111,12 +122,10 @@ Module stdlib_stats_distribution_normal ${t1}$ :: x, y integer :: hz, iz - if( .not. zig_norm_initialized ) call zigset + if(.not. zig_norm_initialized) call zigset iz = 0 - ! original algorithm use 32bit - hz = dist_rand(1_int32) - - iz = iand( hz, 127 ) + hz = dist_rand(1_int32) !32bit random integer + iz = iand( hz, 127 ) !random integer in [0, 127] if( abs( hz ) < kn(iz) ) then res = hz * wn(iz) else @@ -137,8 +146,6 @@ Module stdlib_stats_distribution_normal res = x exit L1 end if - - !original algorithm use 32bit hz = dist_rand(1_int32) iz = iand( hz, 127 ) if( abs( hz ) < kn(iz) ) then @@ -147,103 +154,64 @@ Module stdlib_stats_distribution_normal end if end do L1 end if - return - end function norm_dist_rvs_0_${t1[0]}$${k1}$ + end function rvs_norm_0_${t1[0]}$${k1}$ #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES - impure elemental function norm_dist_rvs_${t1[0]}$${k1}$(loc, scale) & - result(res) + function rvs_norm_${t1[0]}$${k1}$(loc, scale) result(res) + ! ! Normal random variate (loc, scale) ! ${t1}$, intent(in) :: loc, scale ${t1}$ :: res - ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r - ${t1}$ :: x, y - integer :: hz, iz - if(scale==0._${k1}$) call error_stop("Error(norm_dist_rvs): Normal" & - //" distribution scale parameter must be non-zero") - if( .not. zig_norm_initialized ) call zigset - iz = 0 - ! original algorithm use 32bit - hz = dist_rand(1_int32) - - iz = iand( hz, 127 ) - if( abs( hz ) < kn(iz) ) then - res = hz * wn(iz) - else - L1: do - L2: if( iz == 0 ) then - do - x = -log( uni(1.0_${k1}$) ) * rr - y = -log( uni(1.0_${k1}$) ) - if( y + y >= x * x ) exit - end do - res = r + x - if( hz <= 0 ) res = -res - exit L1 - end if L2 - x = hz * wn(iz) - if( fn(iz) + uni(1.0_${k1}$) * (fn(iz-1) - fn(iz)) < & - exp(-HALF * x * x) ) then - res = x - exit L1 - end if - - !original algorithm use 32bit - hz = dist_rand(1_int32) - iz = iand( hz, 127 ) - if( abs( hz ) < kn(iz) ) then - res = hz * wn(iz) - exit L1 - end if - end do L1 - end if + if(scale == 0._${k1}$) call error_stop("Error(rvs_norm): Normal" & + //" distribution scale parameter must be non-zero") + res = rvs_norm_0_${t1[0]}$${k1}$( ) res = res * scale + loc - return - end function norm_dist_rvs_${t1[0]}$${k1}$ + end function rvs_norm_${t1[0]}$${k1}$ #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES - impure elemental function norm_dist_rvs_${t1[0]}$${k1}$(loc, scale) & - result(res) - ! Normal distributed complex. The real part and imaginary part are & + function rvs_norm_${t1[0]}$${k1}$(loc, scale) result(res) + ! + ! Normally distributed complex. The real part and imaginary part are & ! independent of each other. ! ${t1}$, intent(in) :: loc, scale ${t1}$ :: res real(${k1}$) :: tr, ti - tr = norm_dist_rvs_r${k1}$(real(loc), real(scale)) - ti = norm_dist_rvs_r${k1}$(aimag(loc), aimag(scale)) + tr = rvs_norm_r${k1}$(loc % re, scale % re) + ti = rvs_norm_r${k1}$(loc % im, scale % im) res = cmplx(tr, ti, kind=${k1}$) - return - end function norm_dist_rvs_${t1[0]}$${k1}$ + end function rvs_norm_${t1[0]}$${k1}$ #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES - function norm_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & - result(res) + function rvs_norm_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res) ${t1}$, intent(in) :: loc, scale integer, intent(in) :: array_size - ${t1}$, allocatable :: res(:) + ${t1}$ :: res(array_size) ${t1}$, parameter :: r = 3.442619855899_${k1}$, rr = 1.0_${k1}$ / r ${t1}$ :: x, y, re integer :: hz, iz, i - if(scale==0._${k1}$) call error_stop("Error(norm_dist_rvs_array):" & - //" Normal distribution scale parameter must be non-zero") - if( .not. zig_norm_initialized ) call zigset - allocate(res(array_size)) + if(scale == 0._${k1}$) call error_stop("Error(rvs_norm_array): Normal" & + //"distribution scale parameter must be non-zero") + if(.not. zig_norm_initialized) call zigset do i = 1, array_size iz = 0 - ! original algorithm use 32bit hz = dist_rand(1_int32) - iz = iand( hz, 127 ) if( abs( hz ) < kn(iz) ) then re = hz * wn(iz) @@ -266,7 +234,6 @@ Module stdlib_stats_distribution_normal exit L1 end if - !original algorithm use 32bit hz = dist_rand(1_int32) iz = iand( hz, 127 ) if( abs( hz ) < kn(iz) ) then @@ -277,90 +244,90 @@ Module stdlib_stats_distribution_normal end if res(i) = re * scale + loc end do - return - end function norm_dist_rvs_array_${t1[0]}$${k1}$ + end function rvs_norm_array_${t1[0]}$${k1}$ #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES - function norm_dist_rvs_array_${t1[0]}$${k1}$(loc, scale, array_size) & - result(res) + function rvs_norm_array_${t1[0]}$${k1}$(loc, scale, array_size) result(res) ${t1}$, intent(in) :: loc, scale integer, intent(in) :: array_size integer :: i - ${t1}$, allocatable :: res(:) + ${t1}$ :: res(array_size) real(${k1}$) :: tr, ti - allocate(res(array_size)) do i = 1, array_size - tr = norm_dist_rvs_r${k1}$(real(loc), real(scale)) - ti = norm_dist_rvs_r${k1}$(aimag(loc), aimag(scale)) + tr = rvs_norm_r${k1}$(loc % re, scale % re) + ti = rvs_norm_r${k1}$(loc % im, scale % im) res(i) = cmplx(tr, ti, kind=${k1}$) end do - return - end function norm_dist_rvs_array_${t1[0]}$${k1}$ + end function rvs_norm_array_${t1[0]}$${k1}$ #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES - impure elemental function norm_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) & - result(res) - ! Normal distributed probability function + impure elemental function pdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res) + ! + ! Normal distribution probability density function ! ${t1}$, intent(in) :: x, loc, scale real :: res ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$)) - if(scale==0._${k1}$) call error_stop("Error(norm_dist_pdf):" & - //" Normal distribution scale parameter must be non-zero") - res = exp(- 0.5_${k1}$ * (x - loc) * (x - loc) / (scale * scale)) / & + if(scale == 0._${k1}$) call error_stop("Error(pdf_norm): Normal" & + //"distribution scale parameter must be non-zero") + res = exp(- 0.5_${k1}$ * (x - loc) * (x - loc) / (scale * scale)) / & (sqrt_2_Pi * scale) - return - end function norm_dist_pdf_${t1[0]}$${k1}$ + end function pdf_norm_${t1[0]}$${k1}$ #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES - impure elemental function norm_dist_pdf_${t1[0]}$${k1}$(x, loc, scale) & - result(res) + impure elemental function pdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res) ${t1}$, intent(in) :: x, loc, scale real :: res - res = norm_dist_pdf_r${k1}$(real(x), real(loc), real(scale)) - res = res * norm_dist_pdf_r${k1}$(aimag(x), aimag(loc), aimag(scale)) - return - end function norm_dist_pdf_${t1[0]}$${k1}$ + res = pdf_norm_r${k1}$(x % re, loc % re, scale % re) + res = res * pdf_norm_r${k1}$(x % im, loc % im, scale % im) + end function pdf_norm_${t1[0]}$${k1}$ #:endfor + + #:for k1, t1 in REAL_KINDS_TYPES - impure elemental function norm_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) & - result(res) - ! Normal random cumulative distribution function + impure elemental function cdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res) + ! + ! Normal distribution cumulative distribution function ! ${t1}$, intent(in) :: x, loc, scale real :: res ${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$) - if(scale==0._${k1}$) call error_stop("Error(norm_dist_cdf):" & - //" Normal distribution scale parameter must be non-zero") + if(scale == 0._${k1}$) call error_stop("Error(cdf_norm): Normal" & + //"distribution scale parameter must be non-zero") res = (1.0_${k1}$ + erf((x - loc) / (scale * sqrt_2))) / 2.0_${k1}$ - return - end function norm_dist_cdf_${t1[0]}$${k1}$ + end function cdf_norm_${t1[0]}$${k1}$ #:endfor + + #:for k1, t1 in CMPLX_KINDS_TYPES - impure elemental function norm_dist_cdf_${t1[0]}$${k1}$(x, loc, scale) & - result(res) + impure elemental function cdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res) ${t1}$, intent(in) :: x, loc, scale real :: res - res = norm_dist_cdf_r${k1}$(real(x), real(loc), real(scale)) - res = res * norm_dist_cdf_r${k1}$(aimag(x), aimag(loc), aimag(scale)) - return - end function norm_dist_cdf_${t1[0]}$${k1}$ + res = cdf_norm_r${k1}$(x % re, loc % re, scale % re) + res = res * cdf_norm_r${k1}$(x % im, loc % im, scale % im) + end function cdf_norm_${t1[0]}$${k1}$ #:endfor -end module stdlib_stats_distribution_normal \ No newline at end of file +end module stdlib_stats_distribution_normal diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt index 14fed4167..16e33dbc9 100644 --- a/src/tests/stats/CMakeLists.txt +++ b/src/tests/stats/CMakeLists.txt @@ -4,6 +4,7 @@ set(fppFiles test_median.fypp test_distribution_uniform.fypp + test_distribution_normal.fypp ) # Custom preprocessor flags @@ -27,6 +28,7 @@ ADDTEST(var) ADDTEST(varn) ADDTEST(distribution_PRNG) ADDTEST(distribution_uniform) +ADDTEST(distribution_normal) if(DEFINED CMAKE_MAXIMUM_RANK) if(${CMAKE_MAXIMUM_RANK} GREATER 7) diff --git a/src/tests/stats/Makefile.manual b/src/tests/stats/Makefile.manual index 344a7639e..9545eac0f 100644 --- a/src/tests/stats/Makefile.manual +++ b/src/tests/stats/Makefile.manual @@ -1,6 +1,7 @@ SRCFYPP =\ test_median.fypp \ - test_distribution_uniform.fypp + test_distribution_uniform.fypp \ + test_distribution_normal.fypp SRCGEN = $(SRCFYPP:.fypp=.f90) diff --git a/src/tests/stats/test_distribution_normal.fypp b/src/tests/stats/test_distribution_normal.fypp index 261aa5e85..37627056c 100644 --- a/src/tests/stats/test_distribution_normal.fypp +++ b/src/tests/stats/test_distribution_normal.fypp @@ -1,557 +1,235 @@ + +#:include "common.fypp" +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES program test_distribution_normal use stdlib_kinds use stdlib_error, only : check - use stdlib_stats_distribution_PRNG, only : random_seed - use stdlib_stats_distribution_normal, nor_rvs => normal_distribution_rvs, & - nor_pdf => normal_distribution_pdf, & - nor_cdf => normal_distribution_cdf + use stdlib_random, only : random_seed + use stdlib_stats_distribution_uniform, only : uni => rvs_uniform + use stdlib_stats_distribution_normal, only : nor_rvs => rvs_normal, & + nor_pdf => pdf_normal, nor_cdf => cdf_normal implicit none - real(sp), parameter :: sptol = 1000 * epsilon(1.0_sp) - real(dp), parameter :: dptol = 1000 * epsilon(1.0_dp) - real(qp), parameter :: qptol = 1000 * epsilon(1.0_qp) + #:for k1, t1 in REAL_KINDS_TYPES + ${t1}$, parameter :: ${k1}$tol = 1000 * epsilon(1.0_${k1}$) + #:endfor logical :: warn = .true. integer :: put, get - real :: x(2,3,4),a(2,3,4), b(2,3,4) - complex :: loc, scale put = 12345678 call random_seed(put, get) call test_normal_random_generator - call test_nor_rvs_rsp - call test_nor_rvs_rdp - call test_nor_rvs_rqp - call test_nor_rvs_csp - call test_nor_rvs_cdp - call test_nor_rvs_cqp - call test_nor_pdf_rsp - call test_nor_pdf_rdp - call test_nor_pdf_rqp - call test_nor_pdf_csp - call test_nor_pdf_cdp - call test_nor_pdf_cqp + #:for k1, t1 in RC_KINDS_TYPES + call test_nor_rvs_${t1[0]}$${k1}$ + #:endfor + + + + #:for k1, t1 in RC_KINDS_TYPES + call test_nor_pdf_${t1[0]}$${k1}$ + #:endfor + + - call test_nor_cdf_rsp - call test_nor_cdf_rdp - call test_nor_cdf_rqp - call test_nor_cdf_csp - call test_nor_cdf_cdp - call test_nor_cdf_cqp - stop + #:for k1, t1 in RC_KINDS_TYPES + call test_nor_cdf_${t1[0]}$${k1}$ + #:endfor - contains + +contains + subroutine test_normal_random_generator - integer :: i, j, freq(0:1000), num=10000000 + + integer, parameter :: num = 10000000, array_size = 1000 + integer :: i, j, freq(0:array_size - 1) real(dp) :: chisq, expct print *, "" print *, "Test normal random generator with chi-squared" freq = 0 do i = 1, num - j = 1000 * (1 + erf(nor_rvs(0.0, 1.0) / sqrt(2.0))) / 2.0 + j = array_size * (1 + erf(nor_rvs(0.0, 1.0) / sqrt(2.0))) / 2.0 freq(j) = freq(j) + 1 end do chisq = 0.0_dp - expct = num / 1000 - do i = 0, 999 + expct = num / array_size + do i = 0, array_size - 1 chisq = chisq + (freq(i) - expct) ** 2 / expct end do - write(*,*) "The critical values for chi-squared with 1000 d. of f. is" & - //" 1143.92" + write(*,*) "The critical values for chi-squared with 1000 d. of f. is" & + //" 1143.92" write(*,*) "Chi-squared for normal random generator is : ", chisq - call check((chisq < 1143.9), & - msg="normal randomness failed chi-squared test", warn=warn) + call check((chisq < 1143.9), & + msg = "normal randomness failed chi-squared test", warn = warn) end subroutine test_normal_random_generator - subroutine test_nor_rvs_rsp - real(sp) :: res(10), loc, scale - integer :: i, n, k = 5 - integer :: seed, get - real(sp) :: ans(10) = [2.66708039318040679432897377409972250_sp, & - 2.36030794936128329730706809641560540_sp, & - 1.27712218793084242296487218482070602_sp, & - -2.39132544130814794769435138732660562_sp, & - 1.72566595106028652928387145948363468_sp, & - -1.50621775537767632613395107910037041_sp, & - 2.13518835158352082714827702147886157_sp, & - -0.636788253742142318358787633769679815_sp, & - 2.48600787778845799813609573902795091_sp, & - -3.03711473837981227319460231228731573_sp] - - print *, "Test normal_distribution_rvs_rsp" - seed = 25836914 - call random_seed(seed, get) - - loc = 0.5_sp; scale = 2.0_sp - do i = 1, k - res(i) = nor_rvs(loc, scale) ! 2 dummies - end do - res(6:10) = nor_rvs(loc, scale, k) ! 3 dummies - call check(all(abs(res - ans) < sptol), & - msg="normal_distribution_rvs_rsp failed", warn=warn) - end subroutine test_nor_rvs_rsp - subroutine test_nor_rvs_rdp - real(dp) :: res(10), loc, scale - integer :: i, n, k = 5 - integer :: seed, get - real(dp) :: ans(10) = [2.66708039318040679432897377409972250_dp, & - 2.36030794936128329730706809641560540_dp, & - 1.27712218793084242296487218482070602_dp, & - -2.39132544130814794769435138732660562_dp, & - 1.72566595106028652928387145948363468_dp, & - -1.50621775537767632613395107910037041_dp, & - 2.13518835158352082714827702147886157_dp, & - -0.636788253742142318358787633769679815_dp, & - 2.48600787778845799813609573902795091_dp, & - -3.03711473837981227319460231228731573_dp] - - print *, "Test normal_distribution_rvs_rdp" - seed = 25836914 - call random_seed(seed, get) - loc = 0.5_dp; scale = 2.0_dp - do i = 1, k - res(i) = nor_rvs(loc, scale) ! 2 dummies - end do - res(6:10) = nor_rvs(loc, scale, k) ! 3 dummies - call check(all(abs(res - ans) < dptol), & - msg="normal_distribution_rvs_rdp failed", warn=warn) - end subroutine test_nor_rvs_rdp - - subroutine test_nor_rvs_rqp - real(qp) :: res(10), loc, scale - integer :: i, n, k = 5 - integer :: seed, get - real(qp) :: ans(10) = [2.66708039318040679432897377409972250_qp, & - 2.36030794936128329730706809641560540_qp, & - 1.27712218793084242296487218482070602_qp, & - -2.39132544130814794769435138732660562_qp, & - 1.72566595106028652928387145948363468_qp, & - -1.50621775537767632613395107910037041_qp, & - 2.13518835158352082714827702147886157_qp, & - -0.636788253742142318358787633769679815_qp, & - 2.48600787778845799813609573902795091_qp, & - -3.03711473837981227319460231228731573_qp] - - print *, "Test normal_distribution_rvs_rqp" - seed = 25836914 - call random_seed(seed, get) - - loc = 0.5_qp; scale = 2.0_qp - do i = 1, k - res(i) = nor_rvs(loc, scale) ! 2 dummies - end do - res(6:10) = nor_rvs(loc, scale, k) ! 3 dummies - call check(all(abs(res - ans) < qptol), & - msg="normal_distribution_rvs_rqp failed", warn=warn) - end subroutine test_nor_rvs_rqp - - subroutine test_nor_rvs_csp - complex(sp) :: res(10), loc, scale - integer :: i, n, k = 5 - integer :: seed, get - complex(sp) :: ans(10) = [(2.12531029488530509574673033057479188_sp, & - 1.46507698734032082432676702410390135_sp), & - (1.08284164094813181722365413861552952_sp, & - 0.277168639672963013076412153168348595_sp), & - (1.41924946329521489696290359461272601_sp, & - 0.498445561155580918466512230224907398_sp), & - (1.72639126368764062036120776610914618_sp, & - 0.715802936564464420410303091557580046_sp), & - (1.98950590834134349860207180427096318_sp, & - 0.115721315405046931701349421928171068_sp), & - (-1.16929014824793620075382705181255005_sp, & - 0.250744737486995217246033007540972903_sp), & - (1.57160542831869509683428987045772374_sp, & - 0.638282596371312238581197107123443857_sp), & - (-1.36106107654239116833139178197598085_sp, & - 0.166259201494369124318950525776017457_sp), & - (1.13403096805387920698038328737311531_sp, & - 1.04232618148691447146347854868508875_sp), & - (-1.68220535920475811053620418533682823_sp, & - 1.63361446685040256898702182297711261_sp)] - - print *, "Test normal_distribution_rvs_csp" - seed = 25836914 - call random_seed(seed, get) - - loc = (0.5_sp, 1.0_sp); scale = (1.5_sp, 0.5_sp) - do i = 1, k - res(i) = nor_rvs(loc, scale) ! 2 dummies - end do - res(6:10) = nor_rvs(loc, scale, k) ! 3 dummies - call check(all(abs(real(res) - real(ans)) < sptol) .and. & - all(abs(aimag(res) - aimag(ans)) < sptol), & - msg="normal_distribution_rvs_csp failed", warn=warn) - end subroutine test_nor_rvs_csp - - subroutine test_nor_rvs_cdp - complex(dp) :: res(10), loc, scale - integer :: i, n, k = 5 - integer :: seed, get - complex(dp) :: ans(10) = [(2.12531029488530509574673033057479188_dp, & - 1.46507698734032082432676702410390135_dp), & - (1.08284164094813181722365413861552952_dp, & - 0.277168639672963013076412153168348595_dp), & - (1.41924946329521489696290359461272601_dp, & - 0.498445561155580918466512230224907398_dp), & - (1.72639126368764062036120776610914618_dp, & - 0.715802936564464420410303091557580046_dp), & - (1.98950590834134349860207180427096318_dp, & - 0.115721315405046931701349421928171068_dp), & - (-1.16929014824793620075382705181255005_dp, & - 0.250744737486995217246033007540972903_dp), & - (1.57160542831869509683428987045772374_dp, & - 0.638282596371312238581197107123443857_dp), & - (-1.36106107654239116833139178197598085_dp, & - 0.166259201494369124318950525776017457_dp), & - (1.13403096805387920698038328737311531_dp, & - 1.04232618148691447146347854868508875_dp), & - (-1.68220535920475811053620418533682823_dp, & - 1.63361446685040256898702182297711261_dp)] - - print *, "Test normal_distribution_rvs_cdp" - seed = 25836914 - call random_seed(seed, get) - - loc = (0.5_dp, 1.0_dp); scale = (1.5_dp, 0.5_dp) - do i = 1, k - res(i) = nor_rvs(loc, scale) ! 2 dummies - end do - res(6:10) = nor_rvs(loc, scale, k) ! 3 dummies - call check(all(abs(real(res) - real(ans)) < dptol) .and. & - all(abs(aimag(res) - aimag(ans)) < dptol), & - msg="normal_distribution_rvs_cdp failed", warn=warn) - end subroutine test_nor_rvs_cdp - - subroutine test_nor_rvs_cqp - complex(qp) :: res(10), loc, scale - integer :: i, n, k = 5 + #:for k1, t1 in RC_KINDS_TYPES + subroutine test_nor_rvs_${t1[0]}$${k1}$ + ${t1}$ :: res(10), loc, scale + integer, parameter :: k = 5 + integer :: i, n integer :: seed, get - complex(qp) :: ans(10) = [(2.12531029488530509574673033057479188_qp, & - 1.46507698734032082432676702410390135_qp), & - (1.08284164094813181722365413861552952_qp, & - 0.277168639672963013076412153168348595_qp), & - (1.41924946329521489696290359461272601_qp, & - 0.498445561155580918466512230224907398_qp), & - (1.72639126368764062036120776610914618_qp, & - 0.715802936564464420410303091557580046_qp), & - (1.98950590834134349860207180427096318_qp, & - 0.115721315405046931701349421928171068_qp), & - (-1.16929014824793620075382705181255005_qp, & - 0.250744737486995217246033007540972903_qp), & - (1.57160542831869509683428987045772374_qp, & - 0.638282596371312238581197107123443857_qp), & - (-1.36106107654239116833139178197598085_qp, & - 0.166259201494369124318950525776017457_qp), & - (1.13403096805387920698038328737311531_qp, & - 1.04232618148691447146347854868508875_qp), & - (-1.68220535920475811053620418533682823_qp, & - 1.63361446685040256898702182297711261_qp)] - - print *, "Test normal_distribution_rvs_cqp" + #:if t1[0] == "r" + #! for real type + ${t1}$, parameter :: ans(10) = & + [2.66708039318040679432897377409972250_${k1}$, & + 2.36030794936128329730706809641560540_${k1}$, & + 1.27712218793084242296487218482070602_${k1}$, & + -2.39132544130814794769435138732660562_${k1}$, & + 1.72566595106028652928387145948363468_${k1}$, & + -1.50621775537767632613395107910037041_${k1}$, & + 2.13518835158352082714827702147886157_${k1}$, & + -0.636788253742142318358787633769679815_${k1}$, & + 2.48600787778845799813609573902795091_${k1}$, & + -3.03711473837981227319460231228731573_${k1}$] + #:else + #! for complex type + ${t1}$, parameter :: ans(10) = & + [(2.12531029488530509574673033057479188_${k1}$, & + 1.46507698734032082432676702410390135_${k1}$), & + (1.08284164094813181722365413861552952_${k1}$, & + 0.277168639672963013076412153168348595_${k1}$), & + (1.41924946329521489696290359461272601_${k1}$, & + 0.498445561155580918466512230224907398_${k1}$), & + (1.72639126368764062036120776610914618_${k1}$, & + 0.715802936564464420410303091557580046_${k1}$), & + (1.98950590834134349860207180427096318_${k1}$, & + 0.115721315405046931701349421928171068_${k1}$), & + (-1.16929014824793620075382705181255005_${k1}$, & + 0.250744737486995217246033007540972903_${k1}$), & + (1.57160542831869509683428987045772374_${k1}$, & + 0.638282596371312238581197107123443857_${k1}$), & + (-1.36106107654239116833139178197598085_${k1}$, & + 0.166259201494369124318950525776017457_${k1}$), & + (1.13403096805387920698038328737311531_${k1}$, & + 1.04232618148691447146347854868508875_${k1}$), & + (-1.68220535920475811053620418533682823_${k1}$, & + 1.63361446685040256898702182297711261_${k1}$)] + #:endif + + print *, "Test normal_distribution_rvs_${t1[0]}$${k1}$" seed = 25836914 call random_seed(seed, get) + #:if t1[0] == "r" + #! for real type + loc = 0.5_${k1}$; scale = 2.0_${k1}$ + #:else + #! for complex type + loc = (0.5_${k1}$, 1.0_${k1}$); scale = (1.5_${k1}$, 0.5_${k1}$) + #:endif - loc = (0.5_qp, 1.0_qp); scale = (1.5_qp, 0.5_qp) do i = 1, k res(i) = nor_rvs(loc, scale) ! 2 dummies end do res(6:10) = nor_rvs(loc, scale, k) ! 3 dummies - call check(all(abs(real(res) - real(ans)) < qptol) .and. & - all(abs(aimag(res) - aimag(ans)) < qptol), & - msg="normal_distribution_rvs_cqp failed", warn=warn) - end subroutine test_nor_rvs_cqp - - - - subroutine test_nor_pdf_rsp - real(sp) :: x1, x2(3,4), loc, scale - integer :: i, n, k = 5 - integer :: seed, get - real :: res(3,5) - real :: ans(15) = [0.215050772, 0.215050772, 0.215050772, 0.200537622, & - 5.66161536E-02, 0.238986954, 0.265935957,0.262147546,& - 0.249866411, 3.98721099E-02, 0.265902370,0.161311597,& - 0.249177739, 0.237427220, 0.155696079] - - print *, "Test normal_distribution_pdf_rsp" - seed = 741852963 - call random_seed(seed, get) - - loc = -0.5_sp; scale = 1.5_sp - x1 = nor_rvs(loc, scale) - x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) - res(:,1) = nor_pdf(x1, loc, scale) - res(:, 2:5) = nor_pdf(x2, loc, scale) - call check(all(abs(res - reshape(ans, [3,5])) < sptol), & - msg="normal_distribution_pdf_rsp failed", warn=warn) - end subroutine test_nor_pdf_rsp - - subroutine test_nor_pdf_rdp - real(dp) :: x1, x2(3,4), loc, scale - integer :: i, n, k = 5 - integer :: seed, get - real :: res(3,5) - real :: ans(15) = [0.215050772, 0.215050772, 0.215050772, 0.200537622, & - 5.66161536E-02, 0.238986954, 0.265935957,0.262147546,& - 0.249866411, 3.98721099E-02, 0.265902370,0.161311597,& - 0.249177739, 0.237427220, 0.155696079] - - print *, "Test normal_distribution_pdf_rdp" - seed = 741852963 - call random_seed(seed, get) - - loc = -0.5_dp; scale = 1.5_dp - x1 = nor_rvs(loc, scale) - x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) - res(:,1) = nor_pdf(x1, loc, scale) - res(:, 2:5) = nor_pdf(x2, loc, scale) - call check(all(abs(res - reshape(ans, [3,5])) < dptol), & - msg="normal_distribution_pdf_rdp failed", warn=warn) - end subroutine test_nor_pdf_rdp - - subroutine test_nor_pdf_rqp - real(qp) :: x1, x2(3,4), loc, scale - integer :: i, n, k = 5 - integer :: seed, get - real :: res(3,5) - real :: ans(15) = [0.215050772, 0.215050772, 0.215050772, 0.200537622, & - 5.66161536E-02, 0.238986954, 0.265935957,0.262147546,& - 0.249866411, 3.98721099E-02, 0.265902370,0.161311597,& - 0.249177739, 0.237427220, 0.155696079] - - print *, "Test normal_distribution_pdf_rqp" - seed = 741852963 - call random_seed(seed, get) + call check(all(abs(res - ans) < ${k1}$tol), & + msg="normal_distribution_rvs_${t1[0]}$${k1}$ failed", warn=warn) + end subroutine test_nor_rvs_${t1[0]}$${k1}$ - loc = -0.5_qp; scale = 1.5_qp - x1 = nor_rvs(loc, scale) - x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) - res(:,1) = nor_pdf(x1, loc, scale) - res(:, 2:5) = nor_pdf(x2, loc, scale) - call check(all(abs(res - reshape(ans, [3,5])) < qptol), & - msg="normal_distribution_pdf_rqp failed", warn=warn) - end subroutine test_nor_pdf_rqp + #:endfor - subroutine test_nor_pdf_csp - complex(sp) :: x1, x2(3,4), loc, scale - integer :: i, n, k = 5 - integer :: seed, get - real :: res(3,5) - real :: ans(15) = [0.129377320, 0.129377320,0.129377320,4.05915640E-02, & - 0.209143385,2.98881028E-02, 0.128679410, 0.177484736,& - 3.82205322E-02, 7.09915683E-02, 4.56126593E-02, & - 6.57454133E-02,0.165161043,3.86104807E-02,0.196922958] - print *, "Test normal_distribution_pdf_csp" - seed = 741852963 - call random_seed(seed, get) - loc = (-0.5_sp, 0.5_sp); scale = (0.5_sp, 1.5_sp) - x1 = nor_rvs(loc, scale) - x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) - res(:,1) = nor_pdf(x1, loc, scale) - res(:, 2:5) = nor_pdf(x2, loc, scale) - call check(all(abs(res - reshape(ans, [3,5])) < sptol), & - msg="normal_distribution_pdf_csp failed", warn=warn) - end subroutine test_nor_pdf_csp - subroutine test_nor_pdf_cdp - complex(dp) :: x1, x2(3,4), loc, scale - integer :: i, n, k = 5 - integer :: seed, get - real :: res(3,5) - real :: ans(15) = [0.129377320, 0.129377320,0.129377320,4.05915640E-02, & - 0.209143385,2.98881028E-02, 0.128679410, 0.177484736,& - 3.82205322E-02, 7.09915683E-02, 4.56126593E-02, & - 6.57454133E-02,0.165161043,3.86104807E-02,0.196922958] - - print *, "Test normal_distribution_pdf_cdp" - seed = 741852963 - call random_seed(seed, get) + #:for k1, t1 in RC_KINDS_TYPES + subroutine test_nor_pdf_${t1[0]}$${k1}$ - loc = (-0.5_dp, 0.5_dp); scale = (0.5_dp, 1.5_dp) - x1 = nor_rvs(loc, scale) - x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) - res(:,1) = nor_pdf(x1, loc, scale) - res(:, 2:5) = nor_pdf(x2, loc, scale) - call check(all(abs(res - reshape(ans, [3,5])) < dptol), & - msg="normal_distribution_pdf_cdp failed", warn=warn) - end subroutine test_nor_pdf_cdp - - subroutine test_nor_pdf_cqp - complex(qp) :: x1, x2(3,4), loc, scale - integer :: i, n, k = 5 - integer :: seed, get + ${t1}$ :: x1, x2(3,4), loc, scale + integer, parameter :: k = 5 + integer :: i, n, seed, get real :: res(3,5) - real :: ans(15) = [0.129377320, 0.129377320,0.129377320,4.05915640E-02, & - 0.209143385,2.98881028E-02, 0.128679410, 0.177484736,& - 3.82205322E-02, 7.09915683E-02, 4.56126593E-02, & - 6.57454133E-02,0.165161043,3.86104807E-02,0.196922958] - - print *, "Test normal_distribution_pdf_cqp" + #:if t1[0] == "r" + #! for real type + real, parameter :: ans(15) = & + [0.215050772, 0.215050772, 0.215050772, 0.200537622, & + 5.66161536E-02, 0.238986954, 0.265935957,0.262147546,& + 0.249866411, 3.98721099E-02, 0.265902370,0.161311597,& + 0.249177739, 0.237427220, 0.155696079] + #:else + #! for complex type + real, parameter :: ans(15) = & + [0.129377320, 0.129377320,0.129377320,4.05915640E-02, & + 0.209143385,2.98881028E-02, 0.128679410, 0.177484736,& + 3.82205322E-02, 7.09915683E-02, 4.56126593E-02, & + 6.57454133E-02,0.165161043,3.86104807E-02,0.196922958] + #:endif + + print *, "Test normal_distribution_pdf_${t1[0]}$${k1}$" seed = 741852963 call random_seed(seed, get) + #:if t1[0] == "r" + #! for real type + loc = -0.5_${k1}$; scale = 1.5_${k1}$ + #:else + #! for complex type + loc = (-0.5_${k1}$, 0.5_${k1}$); scale = (0.5_${k1}$, 1.5_${k1}$) + #:endif - loc = (-0.5_qp, 0.5_qp); scale = (0.5_qp, 1.5_qp) x1 = nor_rvs(loc, scale) x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) res(:,1) = nor_pdf(x1, loc, scale) res(:, 2:5) = nor_pdf(x2, loc, scale) - call check(all(abs(res - reshape(ans, [3,5])) < qptol), & - msg="normal_distribution_pdf_cqp failed", warn=warn) - end subroutine test_nor_pdf_cqp + call check(all(abs(res - reshape(ans, [3,5])) < ${k1}$tol), & + msg="normal_distribution_pdf_${t1[0]}$${k1}$ failed", warn=warn) + end subroutine test_nor_pdf_${t1[0]}$${k1}$ + #:endfor - subroutine test_nor_cdf_rsp - real(sp) :: x1, x2(3,4), loc, scale - integer :: i, n - integer :: seed, get - real :: res(3,5) - real :: ans(15) = [7.50826299E-02, 7.50826299E-02, 7.50826299E-02, & - 0.143119827, 0.241425425, 0.284345865, 0.233239830, & - 0.341059506,0.353156865,6.81066737E-02,4.38792333E-02,& - 0.763679624, 0.363722175, 0.868187129, 0.626506805] - - print *, "Test normal_distribution_cdf_rsp" - seed = 369147582 - call random_seed(seed, get) - loc = -1.0_sp; scale = 2.0_sp - x1 = nor_rvs(loc, scale) - x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) - res(:,1) = nor_cdf(x1, loc, scale) - res(:, 2:5) = nor_cdf(x2, loc, scale) - call check(all(abs(res - reshape(ans,[3,5])) < sptol), & - msg="normal_distribution_cdf_rsp failed", warn=warn) - end subroutine test_nor_cdf_rsp - subroutine test_nor_cdf_rdp - real(dp) :: x1, x2(3,4), loc, scale - integer :: i, n - integer :: seed, get - real :: res(3,5) - real :: ans(15) = [7.50826299E-02, 7.50826299E-02, 7.50826299E-02, & - 0.143119827, 0.241425425, 0.284345865, 0.233239830, & - 0.341059506,0.353156865,6.81066737E-02,4.38792333E-02,& - 0.763679624, 0.363722175, 0.868187129, 0.626506805] - print *, "Test normal_distribution_cdf_rdp" - seed = 369147582 - call random_seed(seed, get) + #:for k1, t1 in RC_KINDS_TYPES + subroutine test_nor_cdf_${t1[0]}$${k1}$ - loc = -1.0_dp; scale = 2.0_dp - x1 = nor_rvs(loc, scale) - x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) - res(:,1) = nor_cdf(x1, loc, scale) - res(:, 2:5) = nor_cdf(x2, loc, scale) - call check(all(abs(res - reshape(ans,[3,5])) < dptol), & - msg="normal_distribution_cdf_rdp failed", warn=warn) - end subroutine test_nor_cdf_rdp - - subroutine test_nor_cdf_rqp - real(qp) :: x1, x2(3,4), loc, scale - integer :: i, n - integer :: seed, get - real :: res(3,5) - real :: ans(15) = [7.50826299E-02, 7.50826299E-02, 7.50826299E-02, & - 0.143119827, 0.241425425, 0.284345865, 0.233239830, & - 0.341059506,0.353156865,6.81066737E-02,4.38792333E-02,& - 0.763679624, 0.363722175, 0.868187129, 0.626506805] - - print *, "Test normal_distribution_cdf_rqp" - seed = 369147582 - call random_seed(seed, get) - - loc = -1.0_qp; scale = 2.0_qp - x1 = nor_rvs(loc, scale) - x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) - res(:,1) = nor_cdf(x1, loc, scale) - res(:, 2:5) = nor_cdf(x2, loc, scale) - call check(all(abs(res - reshape(ans,[3,5])) < qptol), & - msg="normal_distribution_cdf_rqp failed", warn=warn) - end subroutine test_nor_cdf_rqp - - subroutine test_nor_cdf_csp - complex(sp) :: x1, x2(3,4), loc, scale - integer :: i, n - integer :: seed, get - real :: res(3,5) - real :: ans(15) = [1.07458131E-02, 1.07458131E-02, 1.07458131E-02, & - 6.86483234E-02, 7.95486644E-02, 2.40523387E-02, & - 3.35096754E-02,0.315778911,0.446311295, 0.102010213, & - 7.66918957E-02, 0.564691007, 0.708769500, & - 6.40553832E-02, 5.39999157E-02] - - print *, "Test normal_distribution_cdf_csp" - seed = 369147582 - call random_seed(seed, get) - - loc = (-1.0_sp, 1.0_sp); scale = (1.7_sp, 2.4_sp) - x1 = nor_rvs(loc, scale) - x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) - res(:,1) = nor_cdf(x1, loc, scale) - res(:, 2:5) = nor_cdf(x2, loc, scale) - call check(all(abs(res - reshape(ans,[3,5])) < sptol), & - msg="normal_distribution_cdf_csp failed", warn=warn) - end subroutine test_nor_cdf_csp - - subroutine test_nor_cdf_cdp - complex(dp) :: x1, x2(3,4), loc, scale - integer :: i, n - integer :: seed, get - real :: res(3,5) - real :: ans(15) = [1.07458131E-02, 1.07458131E-02, 1.07458131E-02, & - 6.86483234E-02, 7.95486644E-02, 2.40523387E-02, & - 3.35096754E-02,0.315778911,0.446311295, 0.102010213, & - 7.66918957E-02, 0.564691007, 0.708769500, & - 6.40553832E-02, 5.39999157E-02] - - print *, "Test normal_distribution_cdf_cdp" - seed = 369147582 - call random_seed(seed, get) - - loc = (-1.0_dp, 1.0_dp); scale = (1.7_dp, 2.4_dp) - x1 = nor_rvs(loc, scale) - x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) - res(:,1) = nor_cdf(x1, loc, scale) - res(:, 2:5) = nor_cdf(x2, loc, scale) - call check(all(abs(res - reshape(ans,[3,5])) < dptol), & - msg="normal_distribution_cdf_cdp failed", warn=warn) - end subroutine test_nor_cdf_cdp - - subroutine test_nor_cdf_cqp - complex(qp) :: x1, x2(3,4), loc, scale + ${t1}$ :: x1, x2(3,4), loc, scale integer :: i, n integer :: seed, get real :: res(3,5) - real :: ans(15) = [1.07458131E-02, 1.07458131E-02, 1.07458131E-02, & - 6.86483234E-02, 7.95486644E-02, 2.40523387E-02, & - 3.35096754E-02,0.315778911,0.446311295, 0.102010213, & - 7.66918957E-02, 0.564691007, 0.708769500, & - 6.40553832E-02, 5.39999157E-02] - - print *, "Test normal_distribution_cdf_cqp" + #:if t1[0] == "r" + #!for real type + real, parameter :: ans(15) = & + [7.50826299E-02, 7.50826299E-02, 7.50826299E-02, & + 0.143119827, 0.241425425, 0.284345865, 0.233239830, & + 0.341059506,0.353156865,6.81066737E-02,4.38792333E-02,& + 0.763679624, 0.363722175, 0.868187129, 0.626506805] + #:else + #! for complex type + real, parameter :: ans(15) = & + [1.07458131E-02, 1.07458131E-02, 1.07458131E-02, & + 6.86483234E-02, 7.95486644E-02, 2.40523387E-02, & + 3.35096754E-02,0.315778911,0.446311295, 0.102010213, & + 7.66918957E-02, 0.564691007, 0.708769500, & + 6.40553832E-02, 5.39999157E-02] + #:endif + + print *, "Test normal_distribution_cdf_${t1[0]}$${k1}$" seed = 369147582 call random_seed(seed, get) + #:if t1[0] == "r" + #! for real type + loc = -1.0_${k1}$; scale = 2.0_${k1}$ + #:else + #! for complex type + loc = (-1.0_${k1}$, 1.0_${k1}$); scale = (1.7_${k1}$, 2.4_${k1}$) + #:endif - loc = (-1.0_qp, 1.0_qp); scale = (1.7_qp, 2.4_qp) x1 = nor_rvs(loc, scale) x2 = reshape(nor_rvs(loc, scale, 12), [3,4]) res(:,1) = nor_cdf(x1, loc, scale) res(:, 2:5) = nor_cdf(x2, loc, scale) - call check(all(abs(res - reshape(ans,[3,5])) < qptol), & - msg="normal_distribution_cdf_cqp failed", warn=warn) - end subroutine test_nor_cdf_cqp + call check(all(abs(res - reshape(ans,[3,5])) < ${k1}$tol), & + msg="normal_distribution_cdf_${t1[0]}$${k1}$ failed", warn=warn) + end subroutine test_nor_cdf_${t1[0]}$${k1}$ + #:endfor -end program test_distribution_normal \ No newline at end of file +end program test_distribution_normal From 180da5e5e3837199167b260368abcad2e0d44329 Mon Sep 17 00:00:00 2001 From: Jim-215-Fisher Date: Thu, 7 Oct 2021 15:53:28 -0400 Subject: [PATCH 38/49] add vn variable in zigset subroutine --- src/stdlib_stats_distribution_normal.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index 1a5930b95..a6004bd60 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -86,7 +86,7 @@ contains ! ! Latest version - 1 January 2001 ! - real(dp), parameter :: M1 = 2147483648.0_dp + real(dp), parameter :: M1 = 2147483648.0_dp, vn = 0.00991256303526217_dp real(dp) :: dn, tn, q integer :: i From e786bd501fb6cdfac80fc48512081550eb93ae85 Mon Sep 17 00:00:00 2001 From: Jim-215-Fisher Date: Thu, 7 Oct 2021 16:05:34 -0400 Subject: [PATCH 39/49] change freq array size to array_size in test file --- src/tests/stats/test_distribution_normal.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tests/stats/test_distribution_normal.fypp b/src/tests/stats/test_distribution_normal.fypp index 37627056c..b37aeea28 100644 --- a/src/tests/stats/test_distribution_normal.fypp +++ b/src/tests/stats/test_distribution_normal.fypp @@ -47,7 +47,7 @@ contains subroutine test_normal_random_generator integer, parameter :: num = 10000000, array_size = 1000 - integer :: i, j, freq(0:array_size - 1) + integer :: i, j, freq(0:array_size) real(dp) :: chisq, expct print *, "" From 11bef2fae961477f20eeb31bba9e7552c266fb94 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sat, 9 Oct 2021 20:21:55 -0400 Subject: [PATCH 40/49] Update doc/specs/stdlib_stats_distribution_normal.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_stats_distribution_normal.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index 2dd7d504d..5844af449 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -38,7 +38,7 @@ Function `scale`: optional argument has `intent(in)` and is a scalar of type `real` or `complex`. -`loc` and `scale` arguments must have the same type. +`loc` and `scale` arguments must be of the same type. ### Return value From 62f02a862346573df165c963e9ea530e09199650 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sat, 9 Oct 2021 20:23:05 -0400 Subject: [PATCH 41/49] Update doc/specs/stdlib_stats_distribution_normal.md Co-authored-by: Jeremie Vandenplas --- doc/specs/stdlib_stats_distribution_normal.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index 5844af449..ecfe8f20c 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -68,7 +68,7 @@ program demo_normal_rvs ! -0.633261681 - print *, norm(0.,1.0,10) !an array of 10 standard norml random variates + print *, norm(0.0, 1.0, 10) !an array of 10 standard norml random variates ! -3.38123664E-02 -0.190365672 0.220678389 -0.424612164 -0.249541596 ! 0.865260184 1.11086845 -0.328349441 1.10873628 1.27049923 From 314492d3bab2119f0fe811d42d2ef74258685e28 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Sat, 9 Oct 2021 20:23:16 -0400 Subject: [PATCH 42/49] Update src/stdlib_stats_distribution_normal.fypp Co-authored-by: Jeremie Vandenplas --- src/stdlib_stats_distribution_normal.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index a6004bd60..c72e3276a 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -55,7 +55,7 @@ module stdlib_stats_distribution_normal interface cdf_normal - !! Version experimental + !! version: experimental !! !! Normal Distribution Cumulative Distribution Function !! ([Specification](../page/specs/stdlib_stats_distribution_normal.html# From c4d1391e1b5c4a57cfd38c95e4e014779a98d4e3 Mon Sep 17 00:00:00 2001 From: Jim-215-Fisher Date: Sat, 9 Oct 2021 23:11:25 -0400 Subject: [PATCH 43/49] some modification --- doc/specs/stdlib_stats_distribution_normal.md | 14 +++++++------- src/stdlib_stats_distribution_normal.fypp | 12 +++++------- 2 files changed, 12 insertions(+), 14 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index ecfe8f20c..4b0deb27d 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -64,7 +64,7 @@ program demo_normal_rvs ! 0.563655198 print *, norm(1.0, 2.0) - !normal random variate miu=1.0, sigma=2.0 + !normal random variate mu=1.0, sigma=2.0 ! -0.633261681 @@ -86,8 +86,8 @@ program demo_normal_rvs loc = (-1.0, 2.0) scale = (2.0, 1.0) print *, norm(loc, scale) - !single complex normal random variate with real part of miu=-1, sigma=2; - !imagainary part of miu=2.0 and sigma=1.0 + !single complex normal random variate with real part of mu=-1, sigma=2; + !imagainary part of mu=2.0 and sigma=1.0 ! (1.22566295,2.12518454) @@ -149,7 +149,7 @@ program demo_normal_pdf ! 0.241970733 print *, norm_pdf(2.0,-1.0, 2.0) - !a probability density at 2.0 with miu=-1.0 sigma=2.0 + !a probability density at 2.0 with mu=-1.0 sigma=2.0 !6.47588000E-02 @@ -170,7 +170,7 @@ program demo_normal_pdf scale = (1.0, 2.) print *, norm_pdf((1.5,1.0), loc, scale) ! a complex normal probability density function at (1.5,1.0) with real part - ! of miu=1.0, sigma=1.0 and imaginary part of miu=-0.5, sigma=2.0 + ! of mu=1.0, sigma=1.0 and imaginary part of mu=-0.5, sigma=2.0 ! 5.30100204E-02 @@ -232,7 +232,7 @@ program demo_norm_cdf ! 0.841344714 print *, norm_cdf(2.0, -1.0, 2.0) - ! a cumulative at 2.0 with miu=-1 sigma=2 + ! a cumulative at 2.0 with mu=-1 sigma=2 ! 0.933192849 @@ -253,7 +253,7 @@ program demo_norm_cdf scale = (0.5,1.0) print *, norm_cdf((0.5,-0.5),loc,scale) !complex normal cumulative distribution at (0.5,-0.5) with real part of - !miu=1.0, sigma=0.5 and imaginary part of miu=0.0, sigma=1.0 + !mu=1.0, sigma=0.5 and imaginary part of mu=0.0, sigma=1.0 !4.89511043E-02 diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index c72e3276a..00b7cc92c 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -10,9 +10,9 @@ module stdlib_stats_distribution_normal private real(dp), parameter :: HALF = 0.5_dp, ONE = 1.0_dp, TWO = 2.0_dp - integer, save :: kn(0:127) - real(dp), save :: wn(0:127), fn(0:127) - logical, save :: zig_norm_initialized = .false. + integer :: kn(0:127) + real(dp) :: wn(0:127), fn(0:127) + logical :: zig_norm_initialized = .false. public :: rvs_normal public :: pdf_normal @@ -21,7 +21,7 @@ module stdlib_stats_distribution_normal interface rvs_normal - !! Version experimental + !! version: experimental !! !! Normal Distribution Random Variates !! ([Specification](../page/specs/stdlib_stats_distribution_normal.html# @@ -41,7 +41,7 @@ module stdlib_stats_distribution_normal interface pdf_normal - !! Version experimental + !! version: experimental !! !! Normal Distribution Probability Density Function !! ([Specification](../page/specs/stdlib_stats_distribution_normal.html# @@ -82,8 +82,6 @@ contains ! This is an electronic journal which can be downloaded from: ! http://www.jstatsoft.org/v05/i08 ! - ! N.B. It is assumed that all integers are 32-bit. - ! ! Latest version - 1 January 2001 ! real(dp), parameter :: M1 = 2147483648.0_dp, vn = 0.00991256303526217_dp From ed7c2cbd3fdc9823167be447925f4bd00ee2ffe2 Mon Sep 17 00:00:00 2001 From: Jim-215-Fisher Date: Fri, 3 Dec 2021 21:13:26 -0500 Subject: [PATCH 44/49] modify for complex case --- doc/specs/stdlib_stats_distribution_normal.md | 12 ++++++++++-- src/stdlib_stats_distribution_normal.fypp | 3 ++- src/tests/stats/test_distribution_normal.fypp | 4 ++-- 3 files changed, 14 insertions(+), 5 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index 4b0deb27d..4edf6cd6b 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -102,10 +102,14 @@ Experimental ### Description -The probability density function of the continuous normal distribution. +The probability density function (pdf) of the single real variable normal distribution: $$f(x) = \frac{1}{\sigma \sqrt{2}} e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^{2}}$$ +For complex varible (x + y i) with independent real x and imaginary y parts, the joint probability density function is the product of corresponding marginal pdf of real and imaginary pdf (ref. "Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197): + +$$f(x + y \mathit{i}) = f(x) f(y) = \frac{1}{2\sigma_{x}\sigma_{y}} e^{-\frac{1}{2}[(\frac{x-\mu}{\sigma_{x}})^{2}+(\frac{y-\nu}{\sigma_{y}})^{2}]}$$ + ### Syntax `result = [[stdlib_stats_distribution_normal(module):pdf_normal(interface)]](x, loc, scale)` @@ -185,10 +189,14 @@ Experimental ### Description -Cumulative distribution function of the normal continuous distribution +Cumulative distribution function of the single real variable normal distribution: $$F(x)=\frac{1}{2}\left [ 1+erf(\frac{x-\mu}{\sigma \sqrt{2}}) \right ]$$ +For the complex variable (x + y i) with independent real x and imaginary y parts, the joint cumulative distribution function is the product of corresponding marginal cdf of real and imaginary cdf (ref. "Probability and Random Processes with Applications to Signal Processing and Communications", 2nd ed., Scott L. Miller and Donald Childers, 2012, p.197): + +$$F(x+y\mathit{i})=F(x)F(y)=\frac{1}{4} [1+erf(\frac{x-\mu}{\sigma_{x} \sqrt{2}})] [1+erf(\frac{y-\nu}{\sigma_{y} \sqrt{2}})]$$ + ### Syntax `result = [[stdlib_stats_distribution_normal(module):cdf_normal(interface)]](x, loc, scale)` diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index 00b7cc92c..bf3234162 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -74,7 +74,8 @@ contains subroutine zigset ! Marsaglia & Tsang generator for random normals & random exponentials. - ! Translated from C by Alan Miller (amiller@bigpond.net.au) + ! Translated from C by Alan Miller (amiller@bigpond.net.au), released as public + ! domain (https://jblevins.org/mirror/amiller/) ! ! Marsaglia, G. & Tsang, W.W. (2000) `The ziggurat method for generating ! random variables', J. Statist. Software, v5(8). diff --git a/src/tests/stats/test_distribution_normal.fypp b/src/tests/stats/test_distribution_normal.fypp index b37aeea28..4b8a5eba0 100644 --- a/src/tests/stats/test_distribution_normal.fypp +++ b/src/tests/stats/test_distribution_normal.fypp @@ -2,7 +2,7 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES program test_distribution_normal - use stdlib_kinds + use stdlib_kinds, only : sp, dp, qp use stdlib_error, only : check use stdlib_random, only : random_seed use stdlib_stats_distribution_uniform, only : uni => rvs_uniform @@ -63,7 +63,7 @@ contains chisq = chisq + (freq(i) - expct) ** 2 / expct end do write(*,*) "The critical values for chi-squared with 1000 d. of f. is" & - //" 1143.92" + //" 1143.92" write(*,*) "Chi-squared for normal random generator is : ", chisq call check((chisq < 1143.9), & msg = "normal randomness failed chi-squared test", warn = warn) From 7194bfd7be05fa0a435d19590d7c8cd0017fbe81 Mon Sep 17 00:00:00 2001 From: Jim-215-Fisher Date: Fri, 3 Dec 2021 22:08:05 -0500 Subject: [PATCH 45/49] add xdp in test file --- src/tests/stats/Makefile.manual | 3 ++- src/tests/stats/test_distribution_normal.fypp | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/tests/stats/Makefile.manual b/src/tests/stats/Makefile.manual index 77c286473..15797d509 100644 --- a/src/tests/stats/Makefile.manual +++ b/src/tests/stats/Makefile.manual @@ -1,7 +1,8 @@ SRCFYPP = \ test_mean.fypp \ test_median.fypp \ - test_distribution_uniform.fypp + test_distribution_uniform.fypp \ + test_distribution_normal.fypp SRCGEN = $(SRCFYPP:.fypp=.f90) diff --git a/src/tests/stats/test_distribution_normal.fypp b/src/tests/stats/test_distribution_normal.fypp index 4b8a5eba0..bc87c8ac3 100644 --- a/src/tests/stats/test_distribution_normal.fypp +++ b/src/tests/stats/test_distribution_normal.fypp @@ -2,7 +2,7 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES program test_distribution_normal - use stdlib_kinds, only : sp, dp, qp + use stdlib_kinds, only : sp, dp, xdp, qp use stdlib_error, only : check use stdlib_random, only : random_seed use stdlib_stats_distribution_uniform, only : uni => rvs_uniform From 09f8fd2296e6d6a4c4e2160f23c601f2eafd2b1d Mon Sep 17 00:00:00 2001 From: Jim-215-Fisher Date: Sat, 4 Dec 2021 22:18:20 -0500 Subject: [PATCH 46/49] add precisions for pdf and cdf --- src/stdlib_stats_distribution_normal.fypp | 12 +-- src/tests/stats/test_distribution_normal.fypp | 89 ++++++++++++++----- 2 files changed, 72 insertions(+), 29 deletions(-) diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index bf3234162..92be5b87e 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -1,7 +1,7 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES module stdlib_stats_distribution_normal - use stdlib_kinds + use stdlib_kinds, only :: sp, dp, xdp, qp, int32 use stdlib_error, only : error_stop use stdlib_random, only : dist_rand use stdlib_stats_distribution_uniform, only : uni=>rvs_uniform @@ -274,7 +274,7 @@ contains ! Normal distribution probability density function ! ${t1}$, intent(in) :: x, loc, scale - real :: res + ${t1}$ :: res ${t1}$, parameter :: sqrt_2_pi = sqrt(2.0_${k1}$ * acos(-1.0_${k1}$)) if(scale == 0._${k1}$) call error_stop("Error(pdf_norm): Normal" & @@ -290,7 +290,7 @@ contains #:for k1, t1 in CMPLX_KINDS_TYPES impure elemental function pdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res) ${t1}$, intent(in) :: x, loc, scale - real :: res + real(${k1}$) :: res res = pdf_norm_r${k1}$(x % re, loc % re, scale % re) res = res * pdf_norm_r${k1}$(x % im, loc % im, scale % im) @@ -306,12 +306,12 @@ contains ! Normal distribution cumulative distribution function ! ${t1}$, intent(in) :: x, loc, scale - real :: res + ${t1}$ :: res ${t1}$, parameter :: sqrt_2 = sqrt(2.0_${k1}$) if(scale == 0._${k1}$) call error_stop("Error(cdf_norm): Normal" & //"distribution scale parameter must be non-zero") - res = (1.0_${k1}$ + erf((x - loc) / (scale * sqrt_2))) / 2.0_${k1}$ + res = erfc(- (x - loc) / (scale * sqrt_2)) / 2.0_${k1}$ end function cdf_norm_${t1[0]}$${k1}$ #:endfor @@ -321,7 +321,7 @@ contains #:for k1, t1 in CMPLX_KINDS_TYPES impure elemental function cdf_norm_${t1[0]}$${k1}$(x, loc, scale) result(res) ${t1}$, intent(in) :: x, loc, scale - real :: res + real(${k1}$) :: res res = cdf_norm_r${k1}$(x % re, loc % re, scale % re) res = res * cdf_norm_r${k1}$(x % im, loc % im, scale % im) diff --git a/src/tests/stats/test_distribution_normal.fypp b/src/tests/stats/test_distribution_normal.fypp index bc87c8ac3..72e2746a4 100644 --- a/src/tests/stats/test_distribution_normal.fypp +++ b/src/tests/stats/test_distribution_normal.fypp @@ -146,21 +146,43 @@ contains ${t1}$ :: x1, x2(3,4), loc, scale integer, parameter :: k = 5 integer :: i, n, seed, get - real :: res(3,5) + real(${k1}$) :: res(3,5) #:if t1[0] == "r" #! for real type - real, parameter :: ans(15) = & - [0.215050772, 0.215050772, 0.215050772, 0.200537622, & - 5.66161536E-02, 0.238986954, 0.265935957,0.262147546,& - 0.249866411, 3.98721099E-02, 0.265902370,0.161311597,& - 0.249177739, 0.237427220, 0.155696079] + real(${k1}$), parameter :: ans(15) = & + [0.215050766989949083210785218076278553_${k1}$, & + 0.215050766989949083210785218076278553_${k1}$, & + 0.215050766989949083210785218076278553_${k1}$, & + 0.200537626692880596839299818454439236_${k1}$, & + 5.66161527403268434368575024104022496E-0002_${k1}$, & + 0.238986957612021514867582359518579138_${k1}$, & + 0.265935969411942911029638292783132425_${k1}$, & + 0.262147558654079961109031890374902943_${k1}$, & + 0.249866408914952245533320687656894701_${k1}$, & + 3.98721117498705317877792757313696510E-0002_${k1}$, & + 0.265902369803533466897906694845094995_${k1}$, & + 0.161311603170650092038944290133124635_${k1}$, & + 0.249177740354276111998717092437037695_${k1}$, & + 0.237427217242213206474603807278971527_${k1}$, & + 0.155696086384122017518186260628090478_${k1}$] #:else #! for complex type - real, parameter :: ans(15) = & - [0.129377320, 0.129377320,0.129377320,4.05915640E-02, & - 0.209143385,2.98881028E-02, 0.128679410, 0.177484736,& - 3.82205322E-02, 7.09915683E-02, 4.56126593E-02, & - 6.57454133E-02,0.165161043,3.86104807E-02,0.196922958] + real(${k1}$), parameter :: ans(15) = & + [0.129377311291944176372137325120411497_${k1}$, & + 0.129377311291944176372137325120411497_${k1}$, & + 0.129377311291944176372137325120411497_${k1}$, & + 4.05915662853246811934977653001971736E-0002_${k1}$, & + 0.209143395418940756076861773161637924_${k1}$, & + 2.98881041363874672676853084975547667E-0002_${k1}$, & + 0.128679412679649127469385460133445161_${k1}$, & + 0.177484732473055532384223611956177231_${k1}$, & + 3.82205306942578982084957100753849738E-0002_${k1}$, & + 7.09915714309796034515515428785324918E-0002_${k1}$, & + 4.56126582912124629544443072483362352E-0002_${k1}$, & + 6.57454133967021123696499056531595921E-0002_${k1}$, & + 0.165161039915667041643464172210282279_${k1}$, & + 3.86104822953520989775015755966798359E-0002_${k1}$, & + 0.196922947431391188040943672442575686_${k1}$] #:endif print *, "Test normal_distribution_pdf_${t1[0]}$${k1}$" @@ -193,22 +215,43 @@ contains ${t1}$ :: x1, x2(3,4), loc, scale integer :: i, n integer :: seed, get - real :: res(3,5) + real(${k1}$) :: res(3,5) #:if t1[0] == "r" #!for real type - real, parameter :: ans(15) = & - [7.50826299E-02, 7.50826299E-02, 7.50826299E-02, & - 0.143119827, 0.241425425, 0.284345865, 0.233239830, & - 0.341059506,0.353156865,6.81066737E-02,4.38792333E-02,& - 0.763679624, 0.363722175, 0.868187129, 0.626506805] + real(${k1}$), parameter :: ans(15) = & + [7.50826305038441048487991102776953948E-0002_${k1}$, & + 7.50826305038441048487991102776953948E-0002_${k1}$, & + 7.50826305038441048487991102776953948E-0002_${k1}$, & + 0.143119834108717983250834016885129863_${k1}$, & + 0.241425421525703182028420560765471735_${k1}$, & + 0.284345878626039240974266199229875972_${k1}$, & + 0.233239836366015928845367994433532757_${k1}$, & + 0.341059506137219171082517155967522896_${k1}$, & + 0.353156850199835111081038166086606192_${k1}$, & + 6.81066766396638231790017005897813244E-0002_${k1}$, & + 4.38792331441682923984716366123285346E-0002_${k1}$, & + 0.763679637882860826030745070304416929_${k1}$, & + 0.363722187587355040667876190724308059_${k1}$, & + 0.868187114884980488672309198087692444_${k1}$, & + 0.626506799809652872401992867475200722_${k1}$] #:else #! for complex type - real, parameter :: ans(15) = & - [1.07458131E-02, 1.07458131E-02, 1.07458131E-02, & - 6.86483234E-02, 7.95486644E-02, 2.40523387E-02, & - 3.35096754E-02,0.315778911,0.446311295, 0.102010213, & - 7.66918957E-02, 0.564691007, 0.708769500, & - 6.40553832E-02, 5.39999157E-02] + real(${k1}$), parameter :: ans(15) = & + [1.07458136221563368133842063954746170E-0002_${k1}$, & + 1.07458136221563368133842063954746170E-0002_${k1}$, & + 1.07458136221563368133842063954746170E-0002_${k1}$, & + 6.86483236063879585051085536740820057E-0002_${k1}$, & + 7.95486634025192048896990048539218724E-0002_${k1}$, & + 2.40523393996423661445007940057223384E-0002_${k1}$, & + 3.35096768781160662250307446207445131E-0002_${k1}$, & + 0.315778916661119434962814841317323376_${k1}$, & + 0.446311293878359175362094845206410428_${k1}$, & + 0.102010220821382542292905161748120877_${k1}$, & + 7.66919007012121545175655727052974512E-0002_${k1}$, & + 0.564690968410069125818268877247699603_${k1}$, & + 0.708769523556518785240723539383512333_${k1}$, & + 6.40553790808161720088070925562830659E-0002_${k1}$, & + 5.39999153072107729358158443133850711E-0002_${k1}$] #:endif print *, "Test normal_distribution_cdf_${t1[0]}$${k1}$" From 915d1434f985fb90689d09b9497a8bdedf21d9f0 Mon Sep 17 00:00:00 2001 From: Jim-215-Fisher Date: Sat, 4 Dec 2021 22:30:02 -0500 Subject: [PATCH 47/49] add limitation note in specs --- doc/specs/stdlib_stats_distribution_normal.md | 2 ++ src/stdlib_stats_distribution_normal.fypp | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index 4edf6cd6b..725d523c4 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -22,6 +22,8 @@ With two arguments, the function returns a normal distributed random variate N(l With three arguments, the function returns a rank one array of normal distributed random variates. +Note: the algorithm used for generating normal random variates is fundamentally limited to double precision. + ### Syntax `result = [[stdlib_stats_distribution_normal(module):rvs_normal(interface)]]([loc, scale] [[, array_size]])` diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index 92be5b87e..7ead8f72c 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -1,7 +1,7 @@ #:include "common.fypp" #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES module stdlib_stats_distribution_normal - use stdlib_kinds, only :: sp, dp, xdp, qp, int32 + use stdlib_kinds, only : sp, dp, xdp, qp, int32 use stdlib_error, only : error_stop use stdlib_random, only : dist_rand use stdlib_stats_distribution_uniform, only : uni=>rvs_uniform From 6dfe4772b744ba5cea925a741a3df508437dafc6 Mon Sep 17 00:00:00 2001 From: Jim-215-Fisher Date: Sat, 4 Dec 2021 22:58:10 -0500 Subject: [PATCH 48/49] pdf formula change to prevent overflow --- src/stdlib_stats_distribution_normal.fypp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stdlib_stats_distribution_normal.fypp b/src/stdlib_stats_distribution_normal.fypp index 7ead8f72c..e21b7af14 100644 --- a/src/stdlib_stats_distribution_normal.fypp +++ b/src/stdlib_stats_distribution_normal.fypp @@ -279,7 +279,7 @@ contains if(scale == 0._${k1}$) call error_stop("Error(pdf_norm): Normal" & //"distribution scale parameter must be non-zero") - res = exp(- 0.5_${k1}$ * (x - loc) * (x - loc) / (scale * scale)) / & + res = exp(- 0.5_${k1}$ * ((x - loc) / scale) * (x - loc) / scale) / & (sqrt_2_Pi * scale) end function pdf_norm_${t1[0]}$${k1}$ From 078e5a76e29be2315074302ab1373469bf324de2 Mon Sep 17 00:00:00 2001 From: Jing <53905783+Jim-215-Fisher@users.noreply.github.com> Date: Wed, 8 Dec 2021 09:28:48 -0500 Subject: [PATCH 49/49] Update stdlib_stats_distribution_normal.md --- doc/specs/stdlib_stats_distribution_normal.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/specs/stdlib_stats_distribution_normal.md b/doc/specs/stdlib_stats_distribution_normal.md index 725d523c4..860aec2dd 100644 --- a/doc/specs/stdlib_stats_distribution_normal.md +++ b/doc/specs/stdlib_stats_distribution_normal.md @@ -44,7 +44,7 @@ Function ### Return value -The result is a scalar or rank one array, with a size of `array_size`, of type `real` or `complex`. +The result is a scalar or rank one array, with a size of `array_size`, and as the same type of `scale` and `loc`. ### Example @@ -132,7 +132,7 @@ All three arguments must have the same type. ### Return value -The result is a scalar or an array, with a shape conformable to arguments, of type `real`. +The result is a scalar or an array, with a shape conformable to arguments, and as the same type of input arguments. ### Example @@ -219,7 +219,7 @@ All three arguments must have the same type. ### Return value -The result is a scalar or an array, with a shape conformable to arguments, of type `real`. +The result is a scalar or an array, with a shape conformable to arguments, as the same type of input arguments. ### Example