Return the covariance matrix and mean vector corresponding to the input (potentially weighted) input sample of shape (ndim, nsam) or (nsam, ndim) or a pair of (potentially weighted) time series x(1:nsam) and y(1:nsam) where ndim is the number of data dimensions (the number of data attributes) and nsam is the number of data points.
More...

Detailed Description

Return the covariance matrix and mean vector corresponding to the input (potentially weighted) input sample of shape (ndim, nsam) or (nsam, ndim) or a pair of (potentially weighted) time series x(1:nsam) and y(1:nsam) where ndim is the number of data dimensions (the number of data attributes) and nsam is the number of data points.

See the documentation of the parent module pm_sampleCov for algorithmic details and sample covariance matrix definition.

Parameters

[in,out]	cov	: The output positive semi-definite square matrix of shape `(1:ndim, 1:ndim)` of, type `complex` of kind any supported by the processor (e.g., CK, CK32, CK64, or CK128), type `real` of kind any supported by the processor (e.g., RK, RK32, RK64, or RK128), representing the covariance matrix corresponding to the input `sample` or the pair `x` and `y`, whichever is present. On output, If the input arguments `x` and `y` are present, then `ndim == 2` holds by definition and the output covariance is of the form \begin{equation} \ms{cov} = \begin{bmatrix} \sigma_{xx} && \sigma_{yx} \\ \sigma_{xy} && \sigma_{yy} \end{bmatrix} \end{equation} Otherwise, is `subset` is present, then only the specified input `subset` will be overwritten with the covariance matrix. Any elements not in the specified input `subset` remain intact.
[in]	subset	: The input scalar constant argument that can be any of the following: The constant lowDia, implying that only the lower-diagonal subset of the output covariance matrix must be computed. The constant uppDia, implying that only the upper-diagonal subset of the output covariance matrix must be computed. This input argument is merely serves to resolve the different procedures of this generic interface from each other at compile-time. (optional. It must be present if and only if the input arguments `x` and `y` are both missing.)
[out]	mean	: The output `contiguous` vector of shape `(ndim)` of the same type and kind as the output `cov` containing the `sample` mean. If the input `sample` is missing and `x` and `y` are present, then `mean` must be a vector of size `ndim = 2`. If the input `sample` is present and is a 2D array, then `mean` must be a vector of size `ndim = size(sample, 3 - dim)` (i.e., computed along the specified input dimension `dim`).
[in]	x	: The input `contiguous` vector of shape `(nsam)` of the same type and kind as the output `cov`, containing the first attribute `x` of the observational sample, where `nsam` is the number of observations in the sample. (optional. It must be present if and only if the input argument `y` is present and `sample` is missing.)
[in]	y	: The input `contiguous` vector of shape `(nsam)` of the same type and kind as the output `cov`, containing the second attribute `x` of the observational sample, where `nsam` is the number of observations in the sample. (optional. It must be present if and only if the input argument `x` is present and `sample` is missing.)
[in]	sample	: The input `contiguous` array of shape `(ndim, nsam)` or `(nsam, ndim)` of the same type and kind as the output `cov`, containing the sample comprised of `nsam` observations each with `ndim` attributes. If `sample` is a matrix, then the direction along which the covariance matrix is computed is dictated by the input argument `dim`. (optional. It must be present if and only if the input argument `dim` is present and `x` and `y` are missing.)
[in]	dim	: The input scalar `integer` of default kind IK indicating the dimension of `sample` along which the covariance matrix must be computed. If `dim = 1`, the input `sample` is assumed to have the shape `(nsam, ndim)`. If `dim = 2`, the input `sample` is assumed to have the shape `(ndim, nsam)`. (optional. It must be present if and only if the input argument `sample` is present.)
[in]	weight	: The `contiguous` vector of length `nsam` of, type `integer` of default kind IK, or type `real` of the same kind as the kind of the output `cov`, containing the corresponding weights of individual `nsam` observations in `sample`. (optional. default = getFilled(1, nsam). It can be present if and only if the input arguments `sample` or `x` and `y` are present.)
[out]	weisum	: The output scalar of the same type and kind as the input `weight`, containing `sum(weight)`. (optional. It must be present if and only if the `weight` argument is present.)
[in]	meang	: The input vector of the same type, kind, and rank as the output `cov` of shape `(1:ndim)` containing the best guess for the `sample` mean. If no good guess is known a priori, `meang` can be set to any observation in `sample`. See the example usage below.

Possible calling interfaces ⛓

: use pm_sampleCov, only: setCovMean

! XY time series covariance matrix.

call setCovMean(cov(1:2, 1:2), mean(1:2), x(1:nsam), y(1:nsam) , meang(1:ndim))

call setCovMean(cov(1:2, 1:2), mean(1:2), x(1:nsam), y(1:nsam), weight(1:nsam), weisum, meang(1:ndim))

! sample covariance matrix.

call setCovMean(cov(1:ndim, 1:ndim), subset, mean(1:ndim) , sample(:,:), dim , meang(1:ndim))

call setCovMean(cov(1:ndim, 1:ndim), subset, mean(1:ndim) , sample(:,:), dim, weight(1:nsam), weisum , meang(1:ndim))

pm_sampleCov::setCovMean
Return the covariance matrix and mean vector corresponding to the input (potentially weighted) input ...
Definition: pm_sampleCov.F90:5298

pm_sampleCov
This module contains classes and procedures for computing the properties related to the covariance ma...
Definition: pm_sampleCov.F90:170

Warning: The condition any(x(1) /= x) must hold for the corresponding input arguments.
The condition any(y(1) /= x) must hold for the corresponding input arguments.
The input sample must contain at least two unique values per sample attribute.
The condition 0 < sum(weight) must hold for the corresponding input arguments.
The condition all(0. <= weight) must hold for the corresponding input arguments.
The condition 1 < size(sample, dim) must hold for the corresponding input arguments.
The condition 0 < size(sample, 3 - dim) must hold for the corresponding input arguments.
The condition all(shape(cov) == shape(cor)) must hold for the corresponding input arguments.
The condition 1 <= dim .and. dim <= rank(sample) must hold for the corresponding input arguments.
The condition size(sample, 3 - dim) == size(variance) must hold for the corresponding input arguments.
The condition size(sample, 3 - dim) == size(mean) must hold for the corresponding input arguments.
The condition size(sample, dim) == size(weight) must hold for the corresponding input arguments.
The condition size(x) == size(weight) must hold for the corresponding input arguments.
The condition size(x) == size(y) must hold for the corresponding input arguments.
The condition 1 < size(x) must hold for the corresponding input arguments.
The condition 1 < size(y) must hold for the corresponding input arguments.
These conditions are verified only if the library is built with the preprocessor macro CHECK_ENABLED=1.; The pure procedure(s) documented herein become impure when the ParaMonte library is compiled with preprocessor macro CHECK_ENABLED=1.
By default, these procedures are pure in release build and impure in debug and testing builds.

Note: Note the effects of bias-correction in computing the covariance matrix become noticeable only for very sample sample sizes (i.e., when nsam is small).

See also: getCor
setCor
getCov
setCov
getVar
setVar
getMean
setMean
getCovMerged
setCovMerged
setCovMeanMerged
getVarCorrection
getMeanMerged
setMeanMerged
getVarMerged
setVarMerged

Example usage ⛓

: 1program example

2

3 use pm_kind, only: SK, IK, LK

4 use pm_sampleShift, only: getShifted

5 use pm_arraySpace, only: getLinSpace

6 use pm_arrayResize, only: setResized

7 use pm_distUnif, only: getUnifRand

8 use pm_arrayFill, only: getFilled

9 use pm_sampleCov, only: uppDia

10 use pm_sampleCov, only: lowDia

11 use pm_sampleCov, only: setCovMean

12 use pm_io, only: display_type

13 use pm_io, only: getFormat

14

15 implicit none

16

17 integer(IK) :: itry, ntry = 10

18 type(display_type) :: disp

19 character(:), allocatable :: format

20 disp = display_type(file = "main.out.F90")

21

22 call disp%skip()

23 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

24 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

25 call disp%show("!Compute the covariance matrix of a 2-D sample.")

26 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

27 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

28 call disp%skip()

29

30 block

31 use pm_kind, only: TKG => RKS ! All other real types are also supported.

32 real(TKG), allocatable :: sample(:,:), cov(:,:), mean(:), meang(:)

33 integer(IK) :: ndim, nsam, dim

34 call disp%show("ndim = 2; nsam = 10; dim = 2")

35 ndim = 2; nsam = 10; dim = 2

36 call disp%show("sample = reshape(getUnifRand(1, 20, ndim * nsam), shape = [ndim, nsam], order = [2, 1])")

37 sample = reshape(getUnifRand(1, 20, ndim * nsam), shape = [ndim, nsam], order = [2, 1])

38 call disp%show("sample")

39 call disp%show( sample )

40 call disp%show("meang = sample(:,1)")

41 meang = sample(:,1)

42 call disp%show("mean = getFilled(0., ndim)")

43 mean = getFilled(0., ndim)

44 call disp%show("cov = getFilled(0., ndim, ndim)")

45 cov = getFilled(0., ndim, ndim)

46 call disp%show("call setCovMean(cov, uppDia, mean, sample, dim, meang)")

47 call setCovMean(cov, uppDia, mean, sample, dim, meang)

48 call disp%show("mean")

49 call disp%show( mean )

50 call disp%show("cov")

51 call disp%show( cov )

52 call disp%skip()

53 call disp%show("mean = getFilled(0., ndim)")

54 mean = getFilled(0., ndim)

55 call disp%show("cov = getFilled(0., ndim, ndim)")

56 cov = getFilled(0., ndim, ndim)

57 call disp%show("call setCovMean(cov, lowDia, mean, sample, dim, meang)")

58 call setCovMean(cov, lowDia, mean, sample, dim, meang)

59 call disp%show("mean")

60 call disp%show( mean )

61 call disp%show("cov")

62 call disp%show( cov )

63 call disp%skip()

64 call disp%show("Compute the sample covariance along the first dimension.", deliml = SK_'''')

65 call disp%skip()

66 call disp%show("dim = 1")

67 dim = 1

68 call disp%show("mean = getFilled(0., ndim)")

69 mean = getFilled(0., ndim)

70 call disp%show("cov = getFilled(0., ndim, ndim)")

71 cov = getFilled(0., ndim, ndim)

72 call disp%show("call setCovMean(cov, uppDia, mean, transpose(sample), dim, meang)")

73 call setCovMean(cov, uppDia, mean, transpose(sample), dim, meang)

74 call disp%show("mean")

75 call disp%show( mean )

76 call disp%show("cov")

77 call disp%show( cov )

78 call disp%skip()

79 call disp%show("mean = getFilled(0., ndim)")

80 mean = getFilled(0., ndim)

81 call disp%show("cov = getFilled(0., ndim, ndim)")

82 cov = getFilled(0., ndim, ndim)

83 call disp%show("call setCovMean(cov, lowDia, mean, transpose(sample), dim, meang)")

84 call setCovMean(cov, lowDia, mean, transpose(sample), dim, meang)

85 call disp%show("mean")

86 call disp%show( mean )

87 call disp%show("cov")

88 call disp%show( cov )

89 call disp%skip()

90 call disp%show("Compute the full sample covariance for a pair of time series.", deliml = SK_'''')

91 call disp%skip()

92 call disp%show("mean = getFilled(0., ndim)")

93 mean = getFilled(0., ndim)

94 call disp%show("cov = getFilled(0., ndim, ndim)")

95 cov = getFilled(0., ndim, ndim)

96 call disp%show("call setCovMean(cov, mean, sample(1,:), sample(2,:), meang = sample(1:2,1))")

97 call setCovMean(cov, mean, sample(1,:), sample(2,:), meang = sample(1:2,1))

98 call disp%show("mean")

99 call disp%show( mean )

100 call disp%show("cov")

101 call disp%show( cov )

102 call disp%skip()

103 end block

104

105 call disp%skip()

106 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

107 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

108 call disp%show("!Compute the biased covariance matrix of a weighted 2-D sample.")

109 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

110 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

111 call disp%skip()

112

113 block

114 use pm_kind, only: TKG => RKS ! All other real types are also supported.

115 real(TKG) :: rweisum

116 integer(IK) :: iweisum

117 real(TKG), allocatable :: rweight(:)

118 integer(IK), allocatable :: iweight(:)

119 real(TKG), allocatable :: sample(:,:), cov(:,:), mean(:), meang(:)

120 integer(IK) :: ndim, nsam, dim

121 call disp%show("ndim = 2; nsam = 10; dim = 2")

122 ndim = 2; nsam = 10; dim = 2

123 call disp%show("sample = reshape(getUnifRand(1, 20, ndim * nsam), shape = [ndim, nsam], order = [2, 1])")

124 sample = reshape(getUnifRand(1, 20, ndim * nsam), shape = [ndim, nsam], order = [2, 1])

125 call disp%show("sample")

126 call disp%show( sample )

127 call disp%show("meang = sample(:,1)")

128 meang = sample(:,1)

129 call disp%show("call setResized(mean, ndim)")

130 call setResized(mean, ndim)

131 call disp%show("iweight = getUnifRand(1, 10, nsam) ! integer-valued weights.")

132 iweight = getUnifRand(1, 10, nsam) ! integer-valued weights.

133 call disp%show("iweight")

134 call disp%show( iweight )

135 call disp%show("rweight = iweight ! or real-valued weights.")

136 rweight = iweight ! or real-valued weights.

137 call disp%show("iweight")

138 call disp%show( iweight )

139

140 call disp%skip()

141 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

142 call disp%show("!Compute the covariance matrix integer weights.")

143 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

144 call disp%skip()

145

146 call disp%show("mean = getFilled(0., ndim)")

147 mean = getFilled(0., ndim)

148 call disp%show("cov = getFilled(0., ndim, ndim)")

149 cov = getFilled(0., ndim, ndim)

150 call disp%show("call setCovMean(cov, uppDia, mean, sample, dim, iweight, iweisum, meang)")

151 call setCovMean(cov, uppDia, mean, sample, dim, iweight, iweisum, meang)

152 call disp%show("iweisum")

153 call disp%show( iweisum )

154 call disp%show("mean")

155 call disp%show( mean )

156 call disp%show("cov")

157 call disp%show( cov )

158 call disp%skip()

159 call disp%show("mean = getFilled(0., ndim)")

160 mean = getFilled(0., ndim)

161 call disp%show("cov = getFilled(0., ndim, ndim)")

162 cov = getFilled(0., ndim, ndim)

163 call disp%show("call setCovMean(cov, lowDia, mean, sample, dim, iweight, iweisum, meang)")

164 call setCovMean(cov, lowDia, mean, sample, dim, iweight, iweisum, meang)

165 call disp%show("iweisum")

166 call disp%show( iweisum )

167 call disp%show("mean")

168 call disp%show( mean )

169 call disp%show("cov")

170 call disp%show( cov )

171 call disp%skip()

172 call disp%show("Compute the sample covariance along the first dimension.", deliml = SK_'''')

173 call disp%skip()

174 call disp%show("dim = 1")

175 dim = 1

176 call disp%show("mean = getFilled(0., ndim)")

177 mean = getFilled(0., ndim)

178 call disp%show("cov = getFilled(0., ndim, ndim)")

179 cov = getFilled(0., ndim, ndim)

180 call disp%show("call setCovMean(cov, uppDia, mean, transpose(sample), dim, iweight, iweisum, meang)")

181 call setCovMean(cov, uppDia, mean, transpose(sample), dim, iweight, iweisum, meang)

182 call disp%show("iweisum")

183 call disp%show( iweisum )

184 call disp%show("mean")

185 call disp%show( mean )

186 call disp%show("cov")

187 call disp%show( cov )

188 call disp%skip()

189 call disp%show("mean = getFilled(0., ndim)")

190 mean = getFilled(0., ndim)

191 call disp%show("cov = getFilled(0., ndim, ndim)")

192 cov = getFilled(0., ndim, ndim)

193 call disp%show("call setCovMean(cov, lowDia, mean, transpose(sample), dim, iweight, iweisum, meang)")

194 call setCovMean(cov, lowDia, mean, transpose(sample), dim, iweight, iweisum, meang)

195 call disp%show("iweisum")

196 call disp%show( iweisum )

197 call disp%show("mean")

198 call disp%show( mean )

199 call disp%show("cov")

200 call disp%show( cov )

201 call disp%skip()

202 call disp%show("Compute the full sample covariance for a pair of time series.", deliml = SK_'''')

203 call disp%skip()

204 call disp%show("mean = getFilled(0., ndim)")

205 mean = getFilled(0., ndim)

206 call disp%show("cov = getFilled(0., ndim, ndim)")

207 cov = getFilled(0., ndim, ndim)

208 call disp%show("call setCovMean(cov, mean, sample(1,:), sample(2,:), iweight, iweisum, meang = sample(1:2,1))")

209 call setCovMean(cov, mean, sample(1,:), sample(2,:), iweight, iweisum, meang = sample(1:2,1))

210 call disp%show("iweisum")

211 call disp%show( iweisum )

212 call disp%show("mean")

213 call disp%show( mean )

214 call disp%show("cov")

215 call disp%show( cov )

216 call disp%skip()

217

218 call disp%skip()

219 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

220 call disp%show("!Compute the covariance matrix real weights.")

221 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")

222 call disp%skip()

223

224 call disp%show("dim = 2_IK")

225 dim = 2_IK

226 call disp%show("meang = sample(:,1)")

227 meang = sample(:,1)

228 call disp%show("mean = getFilled(0., ndim)")

229 mean = getFilled(0., ndim)

230 call disp%show("cov = getFilled(0., ndim, ndim)")

231 cov = getFilled(0., ndim, ndim)

232 call disp%show("call setCovMean(cov, uppDia, mean, sample, dim, rweight, rweisum, meang)")

233 call setCovMean(cov, uppDia, mean, sample, dim, rweight, rweisum, meang)

234 call disp%show("rweisum")

235 call disp%show( rweisum )

236 call disp%show("mean")

237 call disp%show( mean )

238 call disp%show("cov")

239 call disp%show( cov )

240 call disp%skip()

241 call disp%show("mean = getFilled(0., ndim)")

242 mean = getFilled(0., ndim)

243 call disp%show("cov = getFilled(0., ndim, ndim)")

244 cov = getFilled(0., ndim, ndim)

245 call disp%show("call setCovMean(cov, lowDia, mean, sample, dim, rweight, rweisum, meang)")

246 call setCovMean(cov, lowDia, mean, sample, dim, rweight, rweisum, meang)

247 call disp%show("rweisum")

248 call disp%show( rweisum )

249 call disp%show("mean")

250 call disp%show( mean )

251 call disp%show("cov")

252 call disp%show( cov )

253 call disp%skip()

254 call disp%show("Compute the sample covariance along the first dimension.", deliml = SK_'''')

255 call disp%skip()

256 call disp%show("dim = 1")

257 dim = 1

258 call disp%show("mean = getFilled(0., ndim)")

259 mean = getFilled(0., ndim)

260 call disp%show("cov = getFilled(0., ndim, ndim)")

261 cov = getFilled(0., ndim, ndim)

262 call disp%show("call setCovMean(cov, uppDia, mean, transpose(sample), dim, rweight, rweisum, meang)")

263 call setCovMean(cov, uppDia, mean, transpose(sample), dim, rweight, rweisum, meang)

264 call disp%show("rweisum")

265 call disp%show( rweisum )

266 call disp%show("mean")

267 call disp%show( mean )

268 call disp%show("cov")

269 call disp%show( cov )

270 call disp%skip()

271 call disp%show("mean = getFilled(0., ndim)")

272 mean = getFilled(0., ndim)

273 call disp%show("cov = getFilled(0., ndim, ndim)")

274 cov = getFilled(0., ndim, ndim)

275 call disp%show("call setCovMean(cov, lowDia, mean, transpose(sample), dim, rweight, rweisum, meang)")

276 call setCovMean(cov, lowDia, mean, transpose(sample), dim, rweight, rweisum, meang)

277 call disp%show("rweisum")

278 call disp%show( rweisum )

279 call disp%show("mean")

280 call disp%show( mean )

281 call disp%show("cov")

282 call disp%show( cov )

283 call disp%skip()

284 call disp%show("Compute the full sample covariance for a pair of time series.", deliml = SK_'''')

285 call disp%skip()

286 call disp%show("mean = getFilled(0., ndim)")

287 mean = getFilled(0., ndim)

288 call disp%show("cov = getFilled(0., ndim, ndim)")

289 cov = getFilled(0., ndim, ndim)

290 call disp%show("call setCovMean(cov, mean, sample(1,:), sample(2,:), rweight, rweisum, meang = sample(1:2,1))")

291 call setCovMean(cov, mean, sample(1,:), sample(2,:), rweight, rweisum, meang = sample(1:2,1))

292 call disp%show("rweisum")

293 call disp%show( rweisum )

294 call disp%show("mean")

295 call disp%show( mean )

296 call disp%show("cov")

297 call disp%show( cov )

298 call disp%skip()

299 end block

300

301end program example

pm_arrayFill::getFilled
Generate and return an array of the specified rank and shape of arbitrary intrinsic type and kind wit...
Definition: pm_arrayFill.F90:132

pm_arrayResize::setResized
Allocate or resize (shrink or expand) an input allocatable scalar string or array of rank 1....
Definition: pm_arrayResize.F90:249

pm_arraySpace::getLinSpace
Generate count evenly spaced points over the interval [x1, x2] if x1 < x2, or [x2,...
Definition: pm_arraySpace.F90:108

pm_distUnif::getUnifRand
Generate and return a scalar or a contiguous array of rank 1 of length s1 of randomly uniformly distr...
Definition: pm_distUnif.F90:4150

pm_io::getFormat
Generate and return a generic or type/kind-specific IO format with the requested specifications that ...
Definition: pm_io.F90:18485

pm_io::show
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11726

pm_io::skip
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11508

pm_sampleShift::getShifted
Generate a sample of shape (nsam), or (ndim, nsam) or (nsam, ndim) that is shifted by the specified i...
Definition: pm_sampleShift.F90:197

pm_arrayFill
This module contains procedures and generic interfaces for convenient allocation and filling of array...
Definition: pm_arrayFill.F90:42

pm_arrayResize
This module contains procedures and generic interfaces for resizing allocatable arrays of various typ...
Definition: pm_arrayResize.F90:81

pm_arraySpace
This module contains procedures and generic interfaces for generating arrays with linear or logarithm...
Definition: pm_arraySpace.F90:33

pm_distUnif
This module contains classes and procedures for computing various statistical quantities related to t...
Definition: pm_distUnif.F90:274

pm_io
This module contains classes and procedures for input/output (IO) or generic display operations on st...
Definition: pm_io.F90:252

pm_io::disp
type(display_type) disp
This is a scalar module variable an object of type display_type for general display.
Definition: pm_io.F90:11393

pm_kind
This module defines the relevant Fortran kind type-parameters frequently used in the ParaMonte librar...
Definition: pm_kind.F90:268

pm_kind::LK
integer, parameter LK
The default logical kind in the ParaMonte library: kind(.true.) in Fortran, kind(....
Definition: pm_kind.F90:541

pm_kind::IK
integer, parameter IK
The default integer kind in the ParaMonte library: int32 in Fortran, c_int32_t in C-Fortran Interoper...
Definition: pm_kind.F90:540

pm_kind::SK
integer, parameter SK
The default character kind in the ParaMonte library: kind("a") in Fortran, c_char in C-Fortran Intero...
Definition: pm_kind.F90:539

pm_kind::RKS
integer, parameter RKS
The single-precision real kind in Fortran mode. On most platforms, this is an 32-bit real kind.
Definition: pm_kind.F90:567

pm_sampleShift
This module contains classes and procedures for shifting univariate or multivariate samples by arbitr...
Definition: pm_sampleShift.F90:59

pm_io::display_type
Generate and return an object of type display_type.
Definition: pm_io.F90:10282

Example Unix compile command via Intel ifort compiler ⛓
1#!/usr/bin/env sh

2rm main.exe

3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe

4./main.exe

Example Windows Batch compile command via Intel ifort compiler ⛓
1del main.exe

2set PATH=..\..\..\lib;%PATH%

3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe

4main.exe

Example Unix / MinGW compile command via GNU gfortran compiler ⛓
1#!/usr/bin/env sh

2rm main.exe

3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe

4./main.exe

Example output ⛓
1

2!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

3!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

4!Compute the covariance matrix of a 2-D sample.

5!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

6!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

7

8ndim = 2; nsam = 10; dim = 2

9sample = reshape(getUnifRand(1, 20, ndim * nsam), shape = [ndim, nsam], order = [2, 1])

10sample

11+4.00000000, +13.0000000, +18.0000000, +7.00000000, +14.0000000, +1.00000000, +6.00000000, +15.0000000, +14.0000000, +16.0000000

12+14.0000000, +13.0000000, +6.00000000, +19.0000000, +5.00000000, +1.00000000, +18.0000000, +10.0000000, +10.0000000, +16.0000000

13meang = sample(:,1)

14mean = getFilled(0., ndim)

15cov = getFilled(0., ndim, ndim)

16call setCovMean(cov, uppDia, mean, sample, dim, meang)

17mean

18+10.8000002, +11.1999998

19cov

20+30.1599979, -1.85999906

21+0.00000000, +31.3600006

22

23mean = getFilled(0., ndim)

24cov = getFilled(0., ndim, ndim)

25call setCovMean(cov, lowDia, mean, sample, dim, meang)

26mean

27+10.8000002, +11.1999998

28cov

29+30.1599979, +0.00000000

30-1.85999906, +31.3600006

31

32'Compute the sample covariance along the first dimension.'

33

34dim = 1

35mean = getFilled(0., ndim)

36cov = getFilled(0., ndim, ndim)

37call setCovMean(cov, uppDia, mean, transpose(sample), dim, meang)

38mean

39+10.8000002, +11.1999998

40cov

41+30.1599979, -1.85999906

42+0.00000000, +31.3600006

43

44mean = getFilled(0., ndim)

45cov = getFilled(0., ndim, ndim)

46call setCovMean(cov, lowDia, mean, transpose(sample), dim, meang)

47mean

48+10.8000002, +11.1999998

49cov

50+30.1599979, +0.00000000

51-1.85999906, +31.3600006

52

53'Compute the full sample covariance for a pair of time series.'

54

55mean = getFilled(0., ndim)

56cov = getFilled(0., ndim, ndim)

57call setCovMean(cov, mean, sample(1,:), sample(2,:), meang = sample(1:2,1))

58mean

59+10.8000002, +11.1999998

60cov

61+30.1599979, -1.85999906

62-1.85999906, +31.3600006

63

64

65!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

66!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

67!Compute the biased covariance matrix of a weighted 2-D sample.

68!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

69!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

70

71ndim = 2; nsam = 10; dim = 2

72sample = reshape(getUnifRand(1, 20, ndim * nsam), shape = [ndim, nsam], order = [2, 1])

73sample

74+11.0000000, +4.00000000, +16.0000000, +5.00000000, +10.0000000, +5.00000000, +4.00000000, +16.0000000, +7.00000000, +15.0000000

75+12.0000000, +9.00000000, +18.0000000, +2.00000000, +17.0000000, +12.0000000, +13.0000000, +5.00000000, +14.0000000, +4.00000000

76meang = sample(:,1)

77call setResized(mean, ndim)

78iweight = getUnifRand(1, 10, nsam) ! integer-valued weights.

79iweight

80+7, +8, +7, +3, +4, +10, +8, +8, +6, +4

81rweight = iweight ! or real-valued weights.

82iweight

83+7, +8, +7, +3, +4, +10, +8, +8, +6, +4

84

85!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

86!Compute the covariance matrix integer weights.

87!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

88

89mean = getFilled(0., ndim)

90cov = getFilled(0., ndim, ndim)

91call setCovMean(cov, uppDia, mean, sample, dim, iweight, iweisum, meang)

92iweisum

93+65

94mean

95+9.04615402, +11.0769234

96cov

97+23.7363319, -1.40355027

98+0.00000000, +20.7479286

99

100mean = getFilled(0., ndim)

101cov = getFilled(0., ndim, ndim)

102call setCovMean(cov, lowDia, mean, sample, dim, iweight, iweisum, meang)

103iweisum

104+65

105mean

106+9.04615402, +11.0769234

107cov

108+23.7363319, +0.00000000

109-1.40355027, +20.7479286

110

111'Compute the sample covariance along the first dimension.'

112

113dim = 1

114mean = getFilled(0., ndim)

115cov = getFilled(0., ndim, ndim)

116call setCovMean(cov, uppDia, mean, transpose(sample), dim, iweight, iweisum, meang)

117iweisum

118+65

119mean

120+9.04615402, +11.0769234

121cov

122+23.7363300, -1.40355027

123+0.00000000, +20.7479286

124

125mean = getFilled(0., ndim)

126cov = getFilled(0., ndim, ndim)

127call setCovMean(cov, lowDia, mean, transpose(sample), dim, iweight, iweisum, meang)

128iweisum

129+65

130mean

131+9.04615402, +11.0769234

132cov

133+23.7363319, +0.00000000

134-1.40355027, +20.7479305

135

136'Compute the full sample covariance for a pair of time series.'

137

138mean = getFilled(0., ndim)

139cov = getFilled(0., ndim, ndim)

140call setCovMean(cov, mean, sample(1,:), sample(2,:), iweight, iweisum, meang = sample(1:2,1))

141iweisum

142+65

143mean

144+9.04615402, +11.0769234

145cov

146+23.7363319, -1.40355027

147-1.40355027, +20.7479286

148

149

150!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

151!Compute the covariance matrix real weights.

152!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

153

154dim = 2_IK

155meang = sample(:,1)

156mean = getFilled(0., ndim)

157cov = getFilled(0., ndim, ndim)

158call setCovMean(cov, uppDia, mean, sample, dim, rweight, rweisum, meang)

159rweisum

160+65.0000000

161mean

162+9.04615402, +11.0769234

163cov

164+23.7363319, -1.40355027

165+0.00000000, +20.7479286

166

167mean = getFilled(0., ndim)

168cov = getFilled(0., ndim, ndim)

169call setCovMean(cov, lowDia, mean, sample, dim, rweight, rweisum, meang)

170rweisum

171+65.0000000

172mean

173+9.04615402, +11.0769234

174cov

175+23.7363319, +0.00000000

176-1.40355027, +20.7479286

177

178'Compute the sample covariance along the first dimension.'

179

180dim = 1

181mean = getFilled(0., ndim)

182cov = getFilled(0., ndim, ndim)

183call setCovMean(cov, uppDia, mean, transpose(sample), dim, rweight, rweisum, meang)

184rweisum

185+65.0000000

186mean

187+9.04615402, +11.0769234

188cov

189+23.7363300, -1.40355027

190+0.00000000, +20.7479286

191

192mean = getFilled(0., ndim)

193cov = getFilled(0., ndim, ndim)

194call setCovMean(cov, lowDia, mean, transpose(sample), dim, rweight, rweisum, meang)

195rweisum

196+65.0000000

197mean

198+9.04615402, +11.0769234

199cov

200+23.7363319, +0.00000000

201-1.40355027, +20.7479305

202

203'Compute the full sample covariance for a pair of time series.'

204

205mean = getFilled(0., ndim)

206cov = getFilled(0., ndim, ndim)

207call setCovMean(cov, mean, sample(1,:), sample(2,:), rweight, rweisum, meang = sample(1:2,1))

208rweisum

209+65.0000000

210mean

211+9.04615402, +11.0769234

212cov

213+23.7363319, -1.40355027

214-1.40355027, +20.7479286

215

216

Benchmarks:

Benchmark :: The runtime performance of setCovMean along different sample dimensions. ⛓

! Test the performance of Cholesky factorization computation using an assumed-shape interface vs. explicit-shape interface.
program benchmark
 
    use pm_kind, only: IK, LK, RKG => RKD, SK
    use pm_sampleCov, only: uppDia
    use pm_bench, only: bench_type
 
    implicit none
 
    integer(IK)                         :: itry, ntry
    integer(IK)                         :: i
    integer(IK)                         :: iarr
    integer(IK)                         :: fileUnit
    integer(IK)     , parameter         :: NARR = 18_IK
    real(RKG)       , allocatable       :: samdim1(:,:)
    real(RKG)       , allocatable       :: samdim2(:,:)
    type(bench_type), allocatable       :: bench(:)
    integer(IK)     , parameter         :: nsammax = 2**NARR
    integer(IK)     , parameter         :: ndim = 5_IK
    real(RKG)                           :: cov(ndim, ndim), mean(ndim)
    integer(IK)                         :: isam, nsam
    real(RKG)                           :: dumm
 
    bench = [ bench_type(name = SK_"setCovMeanDIM1", exec = setCovMeanDIM1, overhead = setOverhead) &
            , bench_type(name = SK_"setCovMeanDIM2", exec = setCovMeanDIM2, overhead = setOverhead) &
            ]
 
    write(*,"(*(g0,:,' '))")
    write(*,"(*(g0,:,' '))") "sample covariance benchmarking..."
    write(*,"(*(g0,:,' '))")
 
    open(newunit = fileUnit, file = "main.out", status = "replace")
 
        write(fileUnit, "(*(g0,:,','))") "nsam", (bench(i)%name, i = 1, size(bench))
 
        dumm = 0._RKG
        loopOverMatrixSize: do iarr = 1, NARR - 1
 
            nsam = 2**iarr
            ntry = nsammax / nsam
            allocate(samdim1(nsam, ndim), samdim2(ndim, nsam))
            write(*,"(*(g0,:,' '))") "Benchmarking setCovMeanDIM1() vs. setCovMeanDIM2()", nsam, ntry
 
            do i = 1, size(bench)
                bench(i)%timing = bench(i)%getTiming()
            end do
 
            write(fileUnit,"(*(g0,:,','))") nsam, (bench(i)%timing%mean / ntry, i = 1, size(bench))
            deallocate(samdim1, samdim2)
 
        end do loopOverMatrixSize
        write(*,"(*(g0,:,' '))") dumm
 
    close(fileUnit)
 
contains
 
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    ! procedure wrappers.
    !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
    subroutine setOverhead()
        do itry = 1, ntry
            call random_number(mean)
            dumm = dumm - mean(1)
        end do
    end subroutine
 
    subroutine setSamDIM1()
        do isam = 1, nsam
            call random_number(samdim1(isam, 1 : ndim))
        end do
    end subroutine
 
    subroutine setSamDIM2()
        do isam = 1, nsam
            call random_number(samdim2(1 : ndim, isam))
        end do
    end subroutine
 
    subroutine setCovMeanDIM1()
        block
            use pm_sampleCov, only: setCovMean
            do itry = 1, ntry
                call setSamDIM1()
                call setCovMean(cov, uppDia, mean, samdim1, dim = 1_IK, meang = samdim1(1,:))
                dumm = dumm + cov(1,1) + mean(1)
            end do
        end block
    end subroutine
 
    subroutine setCovMeanDIM2()
        block
            use pm_sampleCov, only: setCovMean
            do itry = 1, ntry
                call setSamDIM2()
                call setCovMean(cov, uppDia, mean, samdim2, dim = 2_IK, meang = samdim2(:,1))
                dumm = dumm + cov(1,1) + mean(1)
            end do
        end block
    end subroutine
 
end program benchmark

Example Unix compile command via Intel ifort compiler ⛓

#!/usr/bin/env sh
rm main.exe
ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Example Windows Batch compile command via Intel ifort compiler ⛓

del main.exe
set PATH=..\..\..\lib;%PATH%
ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
main.exe

Example Unix / MinGW compile command via GNU gfortran compiler ⛓

#!/usr/bin/env sh
rm main.exe
gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
./main.exe

Postprocessing of the benchmark output ⛓

#!/usr/bin/env python
 
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
 
import os
dirname = os.path.basename(os.getcwd()) 
 
fontsize = 14
 
df = pd.read_csv("main.out", delimiter = ",")
colnames = list(df.columns.values)
 
 
 
ax = plt.figure(figsize = 1.25 * np.array([6.4,4.6]), dpi = 200)
ax = plt.subplot()
 
for colname in colnames[1:]:
    plt.plot( df[colnames[0]].values
            , df[colname].values
            , linewidth = 2
            )
 
plt.xticks(fontsize = fontsize)
plt.yticks(fontsize = fontsize)
ax.set_xlabel(colnames[0], fontsize = fontsize)
ax.set_ylabel("Runtime [ seconds ]", fontsize = fontsize)
ax.set_title(" vs. ".join(colnames[1:])+"\nLower is better.", fontsize = fontsize)
ax.set_xscale("log")
ax.set_yscale("log")
plt.minorticks_on()
plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
ax.tick_params(axis = "y", which = "minor")
ax.tick_params(axis = "x", which = "minor")
ax.legend   ( colnames[1:]
           #, loc='center left'
           #, bbox_to_anchor=(1, 0.5)
            , fontsize = fontsize
            )
 
plt.tight_layout()
plt.savefig("benchmark." + dirname + ".runtime.png")
 
 
 
ax = plt.figure(figsize = 1.25 * np.array([6.4,4.6]), dpi = 200)
ax = plt.subplot()
 
plt.plot( df[colnames[0]].values
        , np.ones(len(df[colnames[0]].values))
        , linestyle = "--"
       #, color = "black"
        , linewidth = 2
        )
for colname in colnames[2:]:
    plt.plot( df[colnames[0]].values
            , df[colname].values / df[colnames[1]].values
            , linewidth = 2
            )
 
plt.xticks(fontsize = fontsize)
plt.yticks(fontsize = fontsize)
ax.set_xlabel(colnames[0], fontsize = fontsize)
ax.set_ylabel("Runtime compared to {}".format(colnames[1]), fontsize = fontsize)
ax.set_title("Runtime Ratio Comparison. Lower means faster.\nLower than 1 means faster than {}().".format(colnames[1]), fontsize = fontsize)
ax.set_xscale("log")
ax.set_yscale("log")
plt.minorticks_on()
plt.grid(visible = True, which = "both", axis = "both", color = "0.85", linestyle = "-")
ax.tick_params(axis = "y", which = "minor")
ax.tick_params(axis = "x", which = "minor")
ax.legend   ( colnames[1:]
           #, bbox_to_anchor = (1, 0.5)
           #, loc = "center left"
            , fontsize = fontsize
            )
 
plt.tight_layout()
plt.savefig("benchmark." + dirname + ".runtime.ratio.png")

Visualization of the benchmark output ⛓

Benchmark moral ⛓

The procedures under the generic interface setCovMean can compute the covariance under two different sample axes.
Recall that C is a row-major language while Fortran is a column-major language.
As such, one would expect the computations for a sample whose observations are stored along the second axis would be faster in the Fortran programming language.
However, such an expectation does not appear to hold at all times and appears to depend significantly on the computing architecture and the number of data attributes involved.
The higher the number of data attributes, the more likely the computations along the second axis of sample will be faster.
Note that for small number of data attributes, the computations along the second data axis involve a small loop that has significant computational cost due to the implicit branching involved in the loop.

Test:: test_pm_sampleCov

Todo:: Normal Priority: The examples of this generic interface should be extended to corrected weighted covariance matrices.

Final Remarks ⛓

If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Copyright: Computational Data Science Lab

Author:: Amir Shahmoradi, Monday March 6, 2017, 3:22 pm, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin.
Fatemeh Bagheri, Monday 02:15 AM, September 27, 2021, Dallas, TX

Definition at line 5298 of file pm_sampleCov.F90.

The documentation for this interface was generated from the following file:

src/fortran/main/pm_sampleCov.F90