Skip to content

Commit becefbc

Browse files
Merge pull request #195 from olekscode/master
Implemented Cholesky Decomposition and the Multivariate Normal Distribution
2 parents 482c790 + 4e48b42 commit becefbc

File tree

14 files changed

+485
-10
lines changed

14 files changed

+485
-10
lines changed

src/BaselineOfPolyMath/BaselineOfPolyMath.class.st

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -113,6 +113,8 @@ BaselineOfPolyMath >> baseline: spec [
113113
with: [ spec requires: #('Math-Complex') ];
114114
package: 'Math-Tests-Core'
115115
with: [ spec requires: #('Math-Core') ];
116+
package: 'Math-Tests-Core-Distribution'
117+
with: [ spec requires: #('Math-Core-Distribution') ];
116118
package: 'Math-Tests-Core-Process'
117119
with: [ spec requires: #('Math-Core-Process') ];
118120
package: 'Math-Tests-Numerical'
@@ -163,7 +165,7 @@ BaselineOfPolyMath >> baseline: spec [
163165
#('Math-Clustering' 'Math-Number-Extensions' 'Math-Chromosome' 'Math-PrincipalComponentAnalysis' 'Math-FunctionFit' 'Math-AutomaticDifferenciation' 'Math-KernelSmoothing' 'Math-Permutation' 'Math-KolmogorovSmirnov');
164166
group: 'Tests'
165167
with:
166-
#('Math-Tests-Matrix' 'Math-Tests-Clustering' 'Math-Tests-Numerical' 'Math-Tests-Complex' 'Math-Tests-Quaternion' 'Math-Tests-Random' 'Math-Tests-ODE' 'Math-Tests-KDTree' 'Math-Tests-FunctionFit' 'Math-Tests-AutomaticDifferenciation' 'Math-Tests-FastFourierTransform' 'Math-Tests-Accuracy' 'Math-Tests-ArbitraryPrecisionFloat' 'Math-Tests-KolmogorovSmirnov' 'Math-Tests-Quantile' 'Math-Tests-Polynomials' 'Math-Tests-PrincipalComponentAnalysis' 'Math-Tests-KernelSmoothing' 'Math-Tests-Number-Extensions' 'Math-Tests-Permutation' 'Math-Tests-TSNE' 'Math-Tests-Core-Process' 'Math-Tests-Core');
168+
#('Math-Tests-Matrix' 'Math-Tests-Clustering' 'Math-Tests-Numerical' 'Math-Tests-Complex' 'Math-Tests-Quaternion' 'Math-Tests-Random' 'Math-Tests-ODE' 'Math-Tests-KDTree' 'Math-Tests-FunctionFit' 'Math-Tests-AutomaticDifferenciation' 'Math-Tests-FastFourierTransform' 'Math-Tests-Accuracy' 'Math-Tests-ArbitraryPrecisionFloat' 'Math-Tests-KolmogorovSmirnov' 'Math-Tests-Quantile' 'Math-Tests-Polynomials' 'Math-Tests-PrincipalComponentAnalysis' 'Math-Tests-KernelSmoothing' 'Math-Tests-Number-Extensions' 'Math-Tests-Permutation' 'Math-Tests-TSNE' 'Math-Tests-Core-Process' 'Math-Tests-Core-Distribution' 'Math-Tests-Core');
167169
group: 'default'
168170
with: #('Core' 'Extensions' 'Tests' 'Benchmarks' 'Accuracy') ]
169171
]

src/Math-Complex/Number.extension.st

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,3 +28,21 @@ Number >> i: aNumber [
2828
aNumber isNumber ifFalse: [self error: 'Badly formed complex number'].
2929
^PMComplex real: self imaginary: aNumber
3030
]
31+
32+
{ #category : #'*Math-Complex' }
33+
Number >> isComplexConjugateOf: aNumber [
34+
"A complex conjugate of a real number is same real number"
35+
^ self = aNumber
36+
]
37+
38+
{ #category : #'*Math-Complex' }
39+
Number >> isComplexConjugateOfAComplexNumber: aComplexNumber [
40+
"A complex conjugate of a real number is same real number"
41+
^ self isComplexConjugateOf: aComplexNumber
42+
]
43+
44+
{ #category : #'*Math-Complex' }
45+
Number >> isRealNumber [
46+
"Answer true if receiver is a real number. All instances of Number are real numbers."
47+
^ true
48+
]

src/Math-Complex/Object.extension.st

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,3 +7,10 @@ Object >> isComplexNumber [
77
^ false
88

99
]
10+
11+
{ #category : #'*Math-Complex' }
12+
Object >> isRealNumber [
13+
"Answer true if receiver is a real number. False by default."
14+
^ false
15+
16+
]

src/Math-Complex/PMComplex.class.st

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -522,6 +522,18 @@ PMComplex >> imaginary [
522522
^ imaginary
523523
]
524524

525+
{ #category : #testing }
526+
PMComplex >> isComplexConjugateOf: aNumber [
527+
"Answer true if self and aNumber are complex conjugates of each other. The complex conjugate of a complex number is the number with an equal real part and an imaginary part equal in magnitude but opposite in sign."
528+
^ aNumber isComplexConjugateOfAComplexNumber: self
529+
]
530+
531+
{ #category : #testing }
532+
PMComplex >> isComplexConjugateOfAComplexNumber: aComplexNumber [
533+
"Answer true if self and aComplexNumber are complex conjugates of each other. The complex conjugate of a complex number is the number with an equal real part and an imaginary part equal in magnitude but opposite in sign."
534+
^ (aComplexNumber real = real) and: [ aComplexNumber imaginary = (-1 * imaginary) ]
535+
]
536+
525537
{ #category : #testing }
526538
PMComplex >> isComplexNumber [
527539
^ true

src/Math-Complex/PMComplex.extension.st

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -17,20 +17,20 @@ PMComplex >> productWithVector: aVector [
1717
^ aVector collect: [ :each | each * self ]
1818
]
1919

20-
{ #category : #'*Math-Complex' }
21-
PMComplex class >> random [
22-
"Answers a random number with abs between 0 and 1."
23-
24-
^ self abs: 1.0 random arg: 2 * Float pi random
25-
]
26-
2720
{ #category : #'*Math-Complex' }
2821
PMComplex >> random [
2922
"analog to Number>>random. However, the only bound is that the abs of the produced complex is less than the length of the receive. The receiver effectively defines a disc within which the random element can be produced."
3023
^ self class random * self
3124

3225
]
3326

27+
{ #category : #'*Math-Complex' }
28+
PMComplex class >> random [
29+
"Answers a random number with abs between 0 and 1."
30+
31+
^ self abs: 1.0 random arg: 2 * Float pi random
32+
]
33+
3434
{ #category : #'*Math-Complex' }
3535
PMComplex >> subtractToPolynomial: aPolynomial [
3636
^ aPolynomial addNumber: self negated
Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
Class {
2+
#name : #PMMultivariateNormalDistribution,
3+
#superclass : #PMProbabilityDensity,
4+
#instVars : [
5+
'meanVector',
6+
'covarianceMatrix'
7+
],
8+
#category : #'Math-Core-Distribution'
9+
}
10+
11+
{ #category : #'instance creation' }
12+
PMMultivariateNormalDistribution class >> meanVector: aMeanVector covarianceMatrix: aCovarianceMatrix [
13+
^ self new
14+
initializeMeanVector: aMeanVector
15+
covarianceMatrix: aCovarianceMatrix;
16+
yourself
17+
]
18+
19+
{ #category : #information }
20+
PMMultivariateNormalDistribution >> average [
21+
^ self meanVector
22+
]
23+
24+
{ #category : #transformation }
25+
PMMultivariateNormalDistribution >> changeParametersBy: aVector [
26+
"Modify the parameters of the receiver by aVector."
27+
meanVector := meanVector + aVector first.
28+
covarianceMatrix := covarianceMatrix + aVector second.
29+
]
30+
31+
{ #category : #information }
32+
PMMultivariateNormalDistribution >> covarianceMatrix [
33+
^ covarianceMatrix
34+
]
35+
36+
{ #category : #information }
37+
PMMultivariateNormalDistribution >> distributionValue: aNumber [
38+
"Answers the probability of observing a random variable distributed according to the receiver with a value lower than or equal to aNumber. Also known as the cumulative density function (CDF)."
39+
40+
self shouldBeImplemented
41+
]
42+
43+
{ #category : #initialization }
44+
PMMultivariateNormalDistribution >> initializeMeanVector: aMeanVector covarianceMatrix: aCovarianceMatrix [
45+
meanVector := aMeanVector.
46+
covarianceMatrix := aCovarianceMatrix.
47+
]
48+
49+
{ #category : #information }
50+
PMMultivariateNormalDistribution >> kurtosis [
51+
"Kurtosis is a measure of the 'tailedness' of the probability distribution of a real-valued random variable. Not defined for a multivariate normal distribution"
52+
self shouldNotImplement
53+
]
54+
55+
{ #category : #information }
56+
PMMultivariateNormalDistribution >> meanVector [
57+
^ meanVector
58+
]
59+
60+
{ #category : #information }
61+
PMMultivariateNormalDistribution >> parameters [
62+
"Returns an Array containing the parameters of the distribution.
63+
It is used to print out the distribution and for fitting."
64+
65+
^ { meanVector . covarianceMatrix }
66+
]
67+
68+
{ #category : #information }
69+
PMMultivariateNormalDistribution >> power [
70+
"Number of dimensions of a multivariate normal distribution"
71+
^ meanVector size
72+
]
73+
74+
{ #category : #information }
75+
PMMultivariateNormalDistribution >> random [
76+
"Answer a vector of random numbers distributed accroding to the receiver."
77+
| standardNormalRandom upperTriangular |
78+
79+
standardNormalRandom := (1 to: self power) asPMVector
80+
collect: [ :each | PMNormalDistribution random ].
81+
82+
upperTriangular := covarianceMatrix choleskyDecomposition.
83+
84+
^ upperTriangular transpose * standardNormalRandom + meanVector
85+
]
86+
87+
{ #category : #information }
88+
PMMultivariateNormalDistribution >> skewness [
89+
"Skewness is a measure of the asymmetry of the probability distribution of a real-valued random variable about its mean. Not defined for a multivariate normal distribution"
90+
self shouldNotImplement
91+
]
92+
93+
{ #category : #information }
94+
PMMultivariateNormalDistribution >> value: aVector [
95+
"Answers the probability that a random variable distributed according to the receiver gives a value between aVector and aVector + espilon (infinitesimal interval)."
96+
97+
^ (-1/2 * (aVector - meanVector) * covarianceMatrix inverse * (aVector - meanVector)) exp / (((2 * Float pi) raisedTo: self power) * covarianceMatrix determinant) sqrt
98+
]

src/Math-Core/PMVector.class.st

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -206,6 +206,12 @@ PMVector >> householder [
206206
^Array with: b with: v
207207
]
208208

209+
{ #category : #testing }
210+
PMVector >> isReal [
211+
"Answer true if all values of the vector are real numbers"
212+
^ self allSatisfy: [ :each | each isRealNumber ].
213+
]
214+
209215
{ #category : #operation }
210216
PMVector >> log [
211217
"Apply log function to every element of a vector"

src/Math-Matrix/PMMatrix.class.st

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -335,6 +335,39 @@ PMMatrix >> atRow: rowIndex put: aCollection startingAt: startColumnNumber [
335335

336336
]
337337

338+
{ #category : #'as yet unclassified' }
339+
PMMatrix >> choleskyDecomposition [
340+
| upperTriangular rowSum partialSum diagonalValue nonDiagonalValue factor |
341+
342+
self isPositiveDefinite ifFalse: [
343+
Error signal: 'Choleski decomposition can only be applied to positive-definite matrices' ].
344+
345+
upperTriangular := self class
346+
zerosRows: self numberOfRows
347+
cols: self numberOfColumns.
348+
349+
1 to: self numberOfRows do: [ :i |
350+
1 to: i do: [ :j |
351+
i = j
352+
ifTrue: [
353+
rowSum := (1 to: j - 1) inject: 0 into: [ :sum :k |
354+
sum + (upperTriangular at: k at: j) squared ].
355+
356+
diagonalValue := ((self at: j at: j) - rowSum) sqrt.
357+
358+
upperTriangular at: j at: j put: diagonalValue ]
359+
ifFalse: [
360+
partialSum := (1 to: j - 1) inject: 0 into: [ :sum :k |
361+
sum + (upperTriangular at: k at: i) * (upperTriangular at: k at: j) ].
362+
363+
factor := upperTriangular at: j at: j.
364+
nonDiagonalValue := ((self at: j at: i) - partialSum) / factor.
365+
366+
upperTriangular at: j at: i put: nonDiagonalValue ] ] ].
367+
368+
^ upperTriangular
369+
]
370+
338371
{ #category : #comparing }
339372
PMMatrix >> closeTo: aPMMatrix [
340373
"Tests that we are within the default Float >> #closeTo: precision of aPMMatrix (0.0001)."
@@ -509,6 +542,54 @@ PMMatrix >> inversePureCRL [
509542
^ self squared inversePureCRL * self transpose
510543
]
511544

545+
{ #category : #testing }
546+
PMMatrix >> isHermitian [
547+
"Hermitian matrix (or self-adjoint matrix) is a complex square matrix that is equal to its own conjugate transpose — that is, the element in the i-th row and j-th column is equal to the complex conjugate of the element in the j-th row and i-th column, for all indices i and j"
548+
549+
self isSquare ifFalse: [ ^ false ].
550+
551+
1 to: self numberOfRows do: [ :i |
552+
1 to: (i - 1) do: [ :j |
553+
((self at: i at: j) isComplexConjugateOf: (self at: j at: i))
554+
ifFalse: [ ^ false ] ] ].
555+
556+
^ true
557+
]
558+
559+
{ #category : #testing }
560+
PMMatrix >> isNegativeDefinite [
561+
"A Hermitian matrix is negative definite if and only if all its eigenvalues are strictly negative"
562+
self isHermitian ifFalse: [ ^ false ].
563+
^ self eigen values allSatisfy: [ :each | each < 0 ]
564+
]
565+
566+
{ #category : #testing }
567+
PMMatrix >> isNegativeSemiDefinite [
568+
"A Hermitian matrix is negative semi-definite if and only if all its eigenvalues are non-positive"
569+
self isHermitian ifFalse: [ ^ false ].
570+
^ self eigen values allSatisfy: [ :each | each <= 0 ]
571+
]
572+
573+
{ #category : #testing }
574+
PMMatrix >> isPositiveDefinite [
575+
"A Hermitian matrix is positive definite if and only if all its eigenvalues are strictly positive"
576+
self isHermitian ifFalse: [ ^ false ].
577+
^ self eigen values allSatisfy: [ :each | each > 0 ]
578+
]
579+
580+
{ #category : #testing }
581+
PMMatrix >> isPositiveSemiDefinite [
582+
"A Hermitian matrix is positive semi-definite if and only if all its eigenvalues are non-negative"
583+
self isHermitian ifFalse: [ ^ false ].
584+
^ self eigen values allSatisfy: [ :each | each >= 0 ]
585+
]
586+
587+
{ #category : #testing }
588+
PMMatrix >> isReal [
589+
"Answer true if all values of the matrix are real numbers"
590+
^ rows allSatisfy: [ :vector | vector isReal ].
591+
]
592+
512593
{ #category : #testing }
513594
PMMatrix >> isSquare [
514595
"Answers true if the number of rows is equal to the number of columns."

src/Math-Matrix/PMSymmetricMatrix.class.st

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@ This class can be instantiated like DhbMatrix via #rows:, but the user has to ma
44
Class {
55
#name : #PMSymmetricMatrix,
66
#superclass : #PMMatrix,
7-
#category : 'Math-Matrix'
7+
#category : #'Math-Matrix'
88
}
99

1010
{ #category : #'instance creation' }
@@ -227,6 +227,12 @@ PMSymmetricMatrix >> inversePureLUP [
227227
ifNotNil: [ :i | ^ self class rows: i ]
228228
]
229229

230+
{ #category : #testing }
231+
PMSymmetricMatrix >> isHermitian [
232+
"Every real symmetric matrix is Hermitian"
233+
^ self isReal
234+
]
235+
230236
{ #category : #testing }
231237
PMSymmetricMatrix >> isSquare [
232238
"Answers true because a symmetric matrix is square."

0 commit comments

Comments
 (0)