Skip to content

Commit cbd6279

Browse files
committed
Implemented a multivariate normal distribution
1 parent b86fbaa commit cbd6279

File tree

5 files changed

+162
-2
lines changed

5 files changed

+162
-2
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
]
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-Distribution/PMNormalDistribution.class.st

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -40,7 +40,8 @@ PMNormalDistribution class >> new: aNumber1 sigma: aNumber2 [
4040

4141
{ #category : #information }
4242
PMNormalDistribution class >> random [
43-
"Answer a random number distributed according to a (0,1) normal distribution."
43+
"A random value generated from the standard normal distribution N(0,1). Can be used to generate values from any other univariate normal distribution N(mu, sigma): first we generate Z ~ N(0,1) and then, to get X ~ N(mu, sigma), we calculate X = sigma * Z + mu (see #random method)"
44+
4445
| v1 v2 w y |
4546
NextRandom isNil
4647
ifTrue: [ [ v1 := Number random * 2 - 1.
Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
Class {
2+
#name : #PMMultivariateNormalDistributionTest,
3+
#superclass : #TestCase,
4+
#instVars : [
5+
'meanVector',
6+
'covarianceMatrix',
7+
'distribution'
8+
],
9+
#category : #'Math-Tests-Core-Distribution'
10+
}
11+
12+
{ #category : #running }
13+
PMMultivariateNormalDistributionTest >> setUp [
14+
super setUp.
15+
16+
meanVector := #(0 1) asPMVector.
17+
18+
covarianceMatrix := PMMatrix rows: #(
19+
(1 0.8)
20+
(0.8 1)).
21+
22+
distribution := PMMultivariateNormalDistribution
23+
meanVector: meanVector
24+
covarianceMatrix: covarianceMatrix.
25+
]
26+
27+
{ #category : #tests }
28+
PMMultivariateNormalDistributionTest >> testAverage [
29+
self assert: distribution average equals: meanVector
30+
]
31+
32+
{ #category : #tests }
33+
PMMultivariateNormalDistributionTest >> testPower [
34+
self assert: distribution power equals: 2
35+
]
36+
37+
{ #category : #tests }
38+
PMMultivariateNormalDistributionTest >> testRandom [
39+
self assert: distribution random size equals: distribution power
40+
]
41+
42+
{ #category : #tests }
43+
PMMultivariateNormalDistributionTest >> testValue [
44+
| interval points expectedPDF |
45+
46+
interval := -1 to: 1 by: 0.5.
47+
48+
points := interval flatCollect: [ :i |
49+
interval collect: [ :j | { i . j } asPMVector ] ].
50+
51+
"Values calculated using Python's scipy"
52+
expectedPDF := #(0.021773722141141552 0.08146294715905042 0.1521928217104107 0.14198250365871087 0.06614272766298167 0.0066868859054493935 0.043603973467720346 0.14198250365871087 0.23086080368580453 0.18744427741405117 0.0010254671663260122 0.011654633617442983 0.06614272766298167 0.18744427741405117 0.26525823848649227 7.852830360809077e-05 0.001555527859401226 0.015386363253590088 0.07599775773306434 0.18744427741405117 3.0028751901985315e-06 0.00010367250011136624 0.0017872959519927589 0.015386363253590088 0.06614272766298167).
53+
54+
1 to: points size do: [ :i |
55+
self
56+
assert: (distribution value: (points at: i))
57+
closeTo: (expectedPDF at: i) ].
58+
]
Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
Package { #name : #'Math-Tests-Core-Distribution' }

0 commit comments

Comments
 (0)