cowan@marob.masa.com (John Cowan) (07/27/90)
"
Purpose of SEQUENCE.ST:
Sequence is an abstract class which is similar to Stream,
but represents unbounded sequences. It is possible to apply
next to a Sequence object arbitrarily many times. There are
three kinds of Sequences here: GeneratedSequence, Random,
and ProbabilityDistribution.
GeneratedSequences are based on a block which computes new
members given the old members. A buffer containing values
already computed is maintained, and the block is given access
to it.
Random returns sequences of random numbers using a decent
congruential algorithm.
ProbabilityDistribution supports various kinds of continuous
and discrete probability distributions.
"
Object subclass: #Sequence
instanceVariableNames: ''
classVariableNames: ''
poolDictionaries: '' !
!Sequence class methods ! !
!Sequence methods !
atEnd
"Answers false, for Stream compatibility."
^false!
close
"Close the sequence. Do nothing."!
contents
"Signals an error, for Stream compatibility."
^self invalidMessage!
next
"Answer the next value in sequence."
^self subclassResponsibility!
next: anInteger
"Answers an Array with the next anInteger
values of the sequence."
| answer |
answer := Array new: anInteger.
1 to: anInteger do: [:i |
answer at: i put: self next].
^answer!
nextPut: anObject
"Signals an error, for Stream compatibility."
^self invalidMessage!
nextPutAll: anObject
"Signals an error, for Stream compatibility."
^self invalidMessage! !
Sequence subclass: #GeneratedSequence
instanceVariableNames:
'buffer function position limit '
classVariableNames: ''
poolDictionaries: '' !
!GeneratedSequence class methods !
from: aBlock
"Answer a new GeneratedSequence that generates
new elements as needed using aBlock. The
block is passed the GeneratedSequence as its
argument."
^super new initialize: aBlock! !
!GeneratedSequence methods !
at: anIndex ifNone: aBlock
"Friends only - Answer the buffer element whose position is
anIndex. Answer nil if no such position in
the buffer."
anIndex > position ifTrue: [^aBlock value].
^buffer at: anIndex!
contents
"Answer the current contents of the GeneratedSequence,
without any new generation."
^buffer copyFrom: 1 to: limit!
function
"Answer the receiver's generation function."
^function!
grow
"Enlarge the buffer."
| newBuffer |
newBuffer := Array new: buffer size * 4 / 3 + 10.
1 to: buffer size do: [:i |
newBuffer at: i put: (buffer at: i)].
buffer := newBuffer!
initialize: aBlock
"Private - Set up instance variables."
function := aBlock.
buffer := Array new: 12.
position := 1.
limit := 0.!
latest: aBlock
"Friends only - Answer the last buffered element.
If no buffered element, evaluate aBlock."
^self at: limit ifNone: aBlock!
limit
"Friends only - Answer the number of elements actually
present in the buffer."
^limit!
next
"Answer the next element. If needed, expand
the buffer or generate new elements."
| value |
position > buffer size ifTrue:
[self grow].
[position > limit] whileTrue:
[value := function value: self.
limit := limit + 1.
buffer at: limit put: value].
value := buffer at: position.
position := position + 1.
^value!
peek
"Answer the next value without changing the
position."
| value |
value := self next.
position := position - 1.
^value!
position
"Answer the current position of the receiver."
^position!
position: anInteger
"Set the receiver's current position to anInteger."
position := anInteger! !
Sequence subclass: #Random
instanceVariableNames:
'seed increment modulus fmodulus multiplier '
classVariableNames:
'Multipliers Moduli Increments '
poolDictionaries: '' !
!Random class methods !
initialize
"Set the recurrence parameters."
Moduli := #(120050 214326 244944 233280
175000 121500 145800).
Multipliers := #(2311 1807 1597 1861
2661 4081 3661).
Increments := #(25367 45289 51749 49297
36979 25673 30809)!
new
"Answer a new random number generator.
Seed is the millisecond clock value."
^super new generator: 1;
seed: Time millisecondClockValue!
new: g seed: s
"Answer a new random number generator."
^super new generator: g;
seed: s! !
!Random methods !
generator: aSmallInteger
"Private - Chooses parameters."
aSmallInteger < 1 | (aSmallInteger > Increments size)
ifTrue:
[self error: 'No such generator'].
increment := Increments at: aSmallInteger.
modulus := Moduli at: aSmallInteger.
fmodulus := modulus asFloat.
multiplier := Multipliers at: aSmallInteger!
next
"Answer the next random number."
seed := seed * multiplier + increment \\ modulus.
^seed asFloat / fmodulus!
seed: aSmallInteger
"Private - Initialize the first random number."
seed := aSmallInteger \\ modulus! !
Sequence subclass: #ProbabilityDistribution
instanceVariableNames: ''
classVariableNames:
'U '
poolDictionaries: '' !
!ProbabilityDistribution class methods !
initialize
"Uniformly distributed random numbers in the range [o,1]."
U := Random new!
new
^self basicNew! !
!ProbabilityDistribution methods !
computeSample: m outOf: n
m>n ifTrue: [^0.0].
^n factorial / (n-m) factorial!
density: x
self subclassResponsibility!
distribution: aCollection
self subclassResponsibility!
inverseDistribution: x
self subclassResponsibility!
next
"This is a general random number generation
method for any probability law;
use the (0,1) uniformly distributed random
variable U as the value of the law's
distribution function. Obtain the next random
value and then solve for the inverse.
The inverse solution is defined by the subclass."
^self inverseDistribution: U next! !
ProbabilityDistribution subclass: #ContinuousProbability
instanceVariableNames: ''
classVariableNames: ''
poolDictionaries: '' !
!ContinuousProbability class methods ! !
!ContinuousProbability methods !
distribution: aCollection
"This is a slow and dirty trapezoidal integration
to determine the area under the probability
function curve y=density (x) for x in the specified
collection. The method assumes that the
collection contains numerically-ordered elements."
| t aStream x1 x2 y1 y2 |
t := 0.0.
aStream := ReadStream on: aCollection.
x2 := aStream next.
y2 := self density: x2.
[x1 := x2. x2 := aStream next]
whileTrue:
[y1 := y2.
y2 := self density: x2.
t := t + ((x2-x1)*(y2+y1))].
^t*0.5! !
ContinuousProbability subclass: #ExponentialDistribution
instanceVariableNames:
'mu '
classVariableNames: ''
poolDictionaries: '' !
!ExponentialDistribution class methods !
mean: p
^self parameter: 1.0/p!
parameter: p
p > 0.0
ifTrue: [^self new setParameter: p]
ifFalse: [self error: 'The probability parameter must be greater than 0.0']! !
!ExponentialDistribution methods !
density: x
x > 0.0
ifTrue: [^mu * (mu*x) negated exp]
ifFalse: [^0.0]!
distribution: anInterval
anInterval last <= 0.0
ifTrue: [^0.0]
ifFalse: [^1.0 - (mu * anInterval last) negated exp - (anInterval first > 0.0 ifTrue: [self distribution: (0.0 to: anInterval first)] ifFalse: [0.0])]!
inverseDistribution: x
^ x ln negated / mu!
mean
^1.0/mu!
setParameter: p
mu := p!
variance
^1.0/(mu*mu)! !
ExponentialDistribution subclass: #GammaDistribution
instanceVariableNames:
'bigN '
classVariableNames: ''
poolDictionaries: '' !
!GammaDistribution class methods !
events: k mean: p
| events |
events := k truncated.
events > 0
ifTrue: [^(self parameter: events/p) setEvents: events]
ifFalse: [self error: 'the number of events must be greater than 0']! !
!GammaDistribution methods !
density: x
| t |
x > 0.0
ifTrue: [t := mu * x.
^(mu raisedTo: bigN) / (self gamma: bigN) *(x raisedTo: bigN - 1) * t negated exp]
ifFalse: [^0.0]!
gamma: n
| t |
t := n - 1.0.
^self computeSample: t outOf: t!
mean
^super mean*bigN!
setEvents: events
bigN := events!
variance
^super variance*bigN! !
ContinuousProbability subclass: #UniformDistribution
instanceVariableNames:
'startNumber stopNumber '
classVariableNames: ''
poolDictionaries: '' !
!UniformDistribution class methods !
from: begin to: end
begin > end
ifTrue: [self error: 'illegal interval']
ifFalse: [^self new setStart: begin toEnd: end]! !
!UniformDistribution methods !
density: x
(x between: startNumber and: stopNumber)
ifTrue: [^1.0 / (stopNumber - startNumber)]
ifFalse: [^0]!
inverseDistribution: x
"x is a random number between 0 and 1"
^startNumber + (x * (stopNumber - startNumber))!
mean
^ (startNumber + stopNumber)/2!
setStart: begin toEnd: end
startNumber := begin.
stopNumber := end!
variance
^ (stopNumber - stopNumber) squared / 12! !
ContinuousProbability subclass: #NormalDistribution
instanceVariableNames:
'mu sigma '
classVariableNames: ''
poolDictionaries: '' !
!NormalDistribution class methods !
mean: a deviation: b
b > 0.0
ifTrue: [^self new setMean: a standardDeviation: b]
ifFalse: [self error: 'standard deviation must be greater than 0.0']! !
!NormalDistribution methods !
density: x
| twoPi t |
twoPi := 2 * 3.1415926536.
t := x - mu/sigma.
^(-0.5 * t squared) exp / (sigma * twoPi sqrt)!
mean
^mu!
next
| v1 v2 s rand u |
rand := UniformDistribution from: -1.0 to: 1.0.
[v1 := rand next.
v2 := rand next.
s := v1 squared + v2 squared.
s >= 1] whileTrue.
u := (-2.0 * s ln /s) sqrt.
^mu + (sigma * v1 *u)!
setMean: m standardDeviation: s
mu := m.
sigma := s!
variance
^sigma squared! !
ProbabilityDistribution subclass: #DiscreteProbability
instanceVariableNames: ''
classVariableNames: ''
poolDictionaries: '' !
!DiscreteProbability class methods ! !
!DiscreteProbability methods !
distribution: aCollection
"Answer the sum of the discrete values of the density function for each element in the collection."
| t |
t := 0.0.
aCollection do: [:i | t := t + (self density: i)].
^t! !
DiscreteProbability subclass: #BernoulliDistribution
instanceVariableNames:
'prob '
classVariableNames: ''
poolDictionaries: '' !
!BernoulliDistribution class methods !
cardgame
"is the first draw of a card an ace?"
(BernoulliDistribution parameter: 4/52) next
"does a car arrive in the next second?"
"will a machine break down today?"!
parameter: aNumber
(aNumber between: 0.0 and: 1.0)
ifTrue: [^self new setParameter: aNumber]
ifFalse: [^self error: 'The probability must be between 0.0 and 1.0']! !
!BernoulliDistribution methods !
density: x
x=1 ifTrue: [^prob].
x=0 ifTrue: [^1.0 - prob].
self error: ' outcomes of a BernoulliDistribution can only be 1 or 0'!
inverseDistribution: x
x <= prob
ifTrue: [^1]
ifFalse: [^0]!
mean
^prob!
setParameter: aNumber
prob :=aNumber!
variance
^prob * (1.0 - prob)! !
BernoulliDistribution subclass: #BinomialDistribution
instanceVariableNames:
'bigN '
classVariableNames: ''
poolDictionaries: '' !
!BinomialDistribution class methods !
events: n mean: m
n truncated <= 0 ifTrue: [self error: 'number of events must be > 0'].
^self new events: n mean: m!
FlippingCoins
| sampleA sampleB |
sampleA := BernoulliDistribution parameter: 0.5.
"Did I get heads?"
sampleA next.
sampleB := BinomialDistribution events: 5 mean: 2.5.
"How many heads did I get in 5 trials?"
sampleB next! !
!BinomialDistribution methods !
density: x
(x between: 0 and: bigN)
ifTrue: [^((self computeSample: x outOf: bigN) / (self computeSample: x outOf: x)) * (prob raisedTo: x)*((1.0 - prob) raisedTo: bigN - x)]
ifFalse: [^0.0]!
events: n mean: m
bigN := n truncated.
self setParameter: m/bigN!
next
|t|
t := 0.
bigN timesRepeat: [t := t + super next].
^t! !
BernoulliDistribution subclass: #GeometricDistribution
instanceVariableNames: ''
classVariableNames: ''
poolDictionaries: '' !
!GeometricDistribution class methods !
CarsArriving
| sample |
"two cars arrive every minute"
sample := GeometricDistribution mean: 60/2.
"what is the probability that it will take 30 sec before the next car arrives?"
sample density: 30.
"Did the next car arrive in 30 to 40 seconds?"
sample distribution: (30 to: 40)!
mean: m
^self parameter: m! !
!GeometricDistribution methods !
density: x
x>0 ifTrue: [^prob * ((1.0 - prob) raisedTo: x - 1)]
ifFalse: [^0.0]!
inverseDistribution: x
^(x ln / (1.0 - prob) ln) ceiling!
mean
^ 1.0 / prob!
variance
^ (1.0 - prob) / prob squared! !
DiscreteProbability subclass: #PoissonDistribution
instanceVariableNames:
'mu '
classVariableNames: ''
poolDictionaries: '' !
!PoissonDistribution class methods !
mean: p
p > 0.0
ifTrue: [^self new setMean: p]
ifFalse: [self error: 'mean must be greater than 0.0']! !
!PoissonDistribution methods !
density: x
x >= 0
ifTrue: [^ ((mu raisedTo: x) * (mu negated exp)) / x factorial]
ifFalse: [^0.0]!
mean
^mu!
next
| p n q |
p := mu negated exp.
n := 0.
q := 1.0.
[q := q * U next.
q >= p]
whileTrue: [n := n + 1].
^n!
setMean: p
mu := p!
variance
^mu! !
DiscreteProbability subclass: #SampleSpace
instanceVariableNames:
'data '
classVariableNames: ''
poolDictionaries: '' !
!SampleSpace class methods !
data: aCollection
^self new setData: aCollection!
heights
| heights |
heights := SampleSpace data: #(60 60 60 62 62 64 64 64 64 66 66 66 68 68 68 68 68 70 70 70).
"what is the probability of randomly selecting a student with height 64?"
heights density: 64.
"what is the probability of randomly selecting a student whose height is between 60 and 64?"
heights distribution: (60 to: 64 by: 2)! !
!SampleSpace methods !
density: x
"x must be in the sample space; the probability must sum over all occurrences of x in the sample space."
(data includes: x)
ifTrue: [^(data occurrencesOf: x) / data size]
ifFalse: [^0]!
inverseDistribution: x
^data at: (x*data size) truncated + 1!
setData: aCollection
data := aCollection! !
--
cowan@marob.masa.com (aka ...!hombre!marob!cowan)
e'osai ko sarji la lojban