Structure and Interpretation of WebPPL
Structure and Interpretation of WebPPL
The language: A functional probabilistic DSL embedded in the "good parts" of JS
WebPPL (pronounced ‘web people’), is a small probabilistic DSL (= domain specific language) embedded in a (pure functional) subset of Javascript. This functional subset is a subset of what Douglas Crockford has identified as "Good Part of JS". WebPPL inherits its syntax from JS and its semantics partly from JS and partly frome SCHEME. The language is developed by the team of Stanford's CoCoLab (Computation & Cognition Lab) with the intention to have a simple to implement and fairly pleasant language to write models in ... . This page documents WebPPL as I understand it with some simple tutorial examples. Misinterpretations and errors are under my responsability. The probability of those is greater than zero due to the fact that I am only an user of the opaque glass box WebPPL but not a developer.
References:
Douglas Crockford, JavaScript: The Good Parts, Sebastapol, CA: O'Reilly, 2008, Print ISBN-13: 978-0-596-51774-8
N. D. Goodman and A. Stuhlmüller (electronic). The Design and Implementation of Probabilistic Programming Languages. Retrieved 2018-3-31 from http://dippl.org
R7RS, the Revised⁷ Report on the Algorithmic Language Scheme
Syntax: My Personal Subset of WebPPL
The official syntax of the WebPPL developer team can be found here. The identifiers of the nonterminals are influenced by the imperative nature of JS. I prefer to describe the WebPPL I used by a personal syntax influenced by the revised 7th Scheme report (R7SR).
'<...>' defines a nonterminal which should be replaced by a nonterminal or terminal (in bold). The meaning of '--->' is that the item left should be replaced by item(s) right. '|' separates alternatives. '*' is the Kleene postfix iterations operator with meaning 'at least zero'. '+' is the Kleene postfix iterations operator with meaning 'at least one'. '/* ... */' denote the start and end of a comment and '//' starts a one-line comment.
- Program - a WebPPL program consists of a sequence of definitions and/or expressions and/or commands. My programs are written in a style most similar to SCHEME. That means as few ';' and return as possible.
- <program> ---> <definition>* | <expression>* | <command>*
- example 1:
- /* a simple functional, nonstochastic WebPPL program */
- var pi = 3.141
- var circumference = function(radius) {2 * pi * radius}
- var circularArea = function(radius) {Math.pow(radius, 2) * pi}
- console.log("circumference(1) === 6.282 --->"); print(circumference(1) === 6.282)
- console.log("circle area(1) === 3.141"); print(circularArea(1) === 3.141)
- /* end of program */
- example 2:
- /* a simple functional, recursive, and stochastic WebPPL program */
- var myGeometricDistribution_2 = function(p) { // recursive sampling function
if (flip(p) === true) {0} // success !
else {1 + myGeometricDistribution_2(p)} // one more failure
} - // end of function
- Definition - a definition is a binding of a symbol with the value of an expression:
- <definition> ---> var <symbol> = <expression>;
- examples: var pi = 3.141, var circumference = function(radius) {2 * pi * radius}
- Expression - an expression evaluates to a value:
- <expression> ---> <literal> | <arithmetic expression> | <boolean expression> | <conditional expression> | function expression> | <call expression> | <array expression> | <object expression>
- examples: 2 * pi * radius, Math.pow(radius, 2) * pi
- Literal - ...
- <literal> ---> <boolean value> | <number> | <character> | <string> | ...
- examples: true, -1, 3.141, 'a', "error"
- Arithmetic Expression - ...
- <arithmetic expression> ---> ...
- examples: 2 * pi * radius, Math.pow(radius, 2) * pi
- Boolean Expression - ...
- <boolean expression> ---> ...
- examples: circumference(1) === 6.282, circularArea(1) === 3.141
- Conditional Expression -
- <conditional expression> ---> if(<boolean expression>) <block> else <block>
- <conditional expression> ---> (<boolean expression>) ? <expression> : <expression>
- example: (x < 0.0) ? -x : x, if(x < 0.0) {-x} else {x}
- Function Expression - ...
- <function expression> ---> function (<parameters>) <block>
- examples: function(radius) {2 * pi * radius}, function(radius) {Math.pow(radius, 2) * pi}
- Call Expression - ...
- <call expression> ---> <symbol>(<arguments>)
- examples: circumference(1), circularArea(1), Math.pow(radius, 2), flip(p)
- Array Expression - ...
- <array expression> ---> [<expressions>]
- Object Expression - ...
- <object expression> ---> {<property-value-pairs>}
- Block - ...
- <block> ---> {<program>}, {0}
- <block> ---> {<program>}, {0}
- Expressions - ...
- <expressions> ---> <expression> | <expression>, <expression>+
- Parameters - ...
- <parameters> ---> <symbol> | <symbol>, <symbol>+
- Arguments - ...
- <arguments> ---> <expression> | <expression>, <expression>+
- Property-Value-Pairs - ...
- <property-value-pairs> ---> <property-value-pair>* | <property-value-pair> , <property-value-pair>+
- Property-Value-Pair - ...
- <property-value-pair> ---> <symbol> : <expression>
- Command - a command is a procedure that does not return useful values to its continuation.
- <command> ---> console.log(...) | print(...) | display(...) | <return command>
- Return Command - ...
- <return command> ---> return <expression>
Syntax: Vector and Matrix Multiplication
We chose an example from linear algebra to demonstrate that functional WebPPL is sufficient expressive to deal even with matrix problems. These belong to the basic tools of a machine learner or data scientist.
Chapter 2 of Goodfellow et al. book Deep Learning is about linear algebra. Subchapter 2.1 deals with Scalars, Vectors, Matrices, and Tensors and subchapter 2.2 about Multiplying Matrices and Vectors.
Definition: The matrix product of matrices A and B is a third matrix C. In order for this product to be defined, A must have the same number of columns as B has rows. If A is of shape m x n and B is of shape n x p, then C is of shape m x p. We can write the matrix product just by placing two or more matrices together, for example
$$\mathbf{_mC_p} =\mathbf{_mA_n} \cdot \mathbf{_nB_p}$$
The product operation is defined by the scalar or dot product of the i-th row vector a'i of A with the j-th column vector bj of B:
$$ c_{ij} = \mathbf{a'}_i \cdot \mathbf{b}_j = a_{i1}\cdot b_{1j} + a_{i2}\cdot b_{2j} + ... + a_{in}\cdot b_{nj} = \sum\limits_{k=1}^{n}a_{ik}b_{kj}$$
Goodfellow et al.(2016, p.32), Kowalsky & Michler (2003, p.48).
Example (Kowalsky & Michler, 2003, Example 3.1.19, p.48):
$$\mathbf{_2A_3} = \left(\begin{array}{rrr} 1&-2&0 \\ 0&1&1 \end{array}\right)$$
$$\mathbf{_3B_3} = \left(\begin{array}{rrr} 2&1&1 \\ 1&-1&0 \\0&0&1 \end{array}\right)$$
$$\mathbf{_2C_3} = \left(\begin{array}{rrr} 0&3&1 \\ 1&-1&1 \end{array}\right)$$
References:
Abelson, H., Sussman, G.J., & Sussman, J., Structure and Interpretation of Computer Programs, MIT Press, 2003, 2/e, ISBN 0-262-01153-0
Brown, Jonathon D., Linear Models in Matrix Form, Cham, Heidelberg, New York, Dordrecht, London: Springer, 2014, ISBN 978-3-319-34569-7
Goodfellow, Bengio & Courville, Deep Learning, Cambridge, Mass.: MIT Press, 2016, ISBN 978-0-262-03561-3
Kowalsky & Michler, Lineare Algebra, Berlin-New York: Walter De Gruyter, 2003, 12/e, ISBN 3-11-017963-6
Syntax: Functional Vector and Matrix Multiplication in Core WebPPL
console.log("==============================================================================")
console.log("file: PCM20180327_FunctionalMatrixMultiplication *** PCM *** *2018/04/05* ")
console.log(" - definition is from: Kowalsky & Michler, Lineare Algebra, Berlin & ")
console.log(" Walter de Gruyter, 2003, 12/e, p.48 ")
console.log(" - example is from: Kowalsky & Michler, Lineare Algebra, 2003, p.48 ")
console.log(" - further inspiration in: Abelson, Sussman & Sussman, Structure and ")
console.log(" Interpretation of Computer Programs, ch. 2.2.3, ")
console.log(" 'Sequences as Conventional Interfaces', MIT Press, ")
console.log(" 1996, 2/e, p.113 - 126 ")
console.log(" - application is in: Goodfellow, Bengio & Courville, Deep Learning, ")
console.log(" Cambridge, Mass.: MIT Press, 2016, p.32 ")
console.log("==============================================================================")
var startTime = Date.now()
var errorMsg01 = 'errorMsg01: vectors have incompatible dimensions !'
var errorMsg02 = "errorMsg02: vectors have incompatible types !"
var errorMsg03 = "errorMsg03: vector and matrix have incompatible dimensions !"
var errorMsg04 = "errorMsg04: vector and matrix have incompatible types !"
var errorMsg05 = "errorMsg05: matrices have incompatible dimensions !"
var errorMsg06 = "errorMsg06: matrices have incompatible types !"
var errorMsg07 = "errorMsg07: objects have incompatible types !"
console.log("==============================================================================")
var runTime = function(startTime) {
var stopTime = Date.now()
var runTimeInMsec = (stopTime - startTime)
var runTimeInSec = runTimeInMsec/1000.
var runTimeInMin = runTimeInSec/60.
console.log("runTimeInMsec ---> ", runTimeInMsec)
console.log("runTimeInSec ---> ", runTimeInSec)
console.log("runTimeInMin ---> ", runTimeInMin)
}
//-------------------------------------------------------------------------------------------
var makeVector = function(id, type, array) {
if (type == 'rowVector')
{return {id:id, type:'rowVector', nrows:1, ncols:array.length, vector: array}}
else if (type == 'colVector')
{return {id:id, type:'colVector', nrows:array.length, ncols:1, vector: array}}
else "makeVector does not understand type"
}
//-------------------------------------------------------------------------------------------
var makeMatrix = function(id, type, array) { // input array is rowwise
if (type == 'rowMatrix') {
var mrows = array.length
var mcols = array[0].vector.length
return {id:id, type:'rowMatrix', nrows:mrows, ncols:mcols, matrix:array}}
else if(type == 'colMatrix') {
var mrows = array.length
var mcols = array[0].vector.length
var inc0 = function(j) {j + 0}
var colIndices = mapN(inc0, mcols) // [0, 1, 2, ...]
var transposedArray = map(function(colIndexJ) {
makeVector(id.concat(colIndexJ), 'colVector',
map(function(rowVector){rowVector.vector[colIndexJ]}, array))
}, colIndices)
return {id:id, type:'colMatrix', nrows:mrows, ncols:mcols, matrix:transposedArray}
}
else return {msg:"makeMatrix does not understand type", type:type}
}
//-------------------------------------------------------------------------------------------
var matrixMultiplication = function(xObject, yObject) {
//-----------------------------------------------------------------------------------
// definitions of local functions
//-----------------------------------------------------------------------------------
var scalarProduct = function(xVectorObject, yVectorObject) {
(xVectorObject.type === 'rowVector' && yVectorObject.type === 'colVector') ?
(xVectorObject.ncols === yVectorObject.nrows) ?
sum(map2(function(xi, yi){ xi * yi}, xVectorObject.vector, yVectorObject.vector)) :
{msg:errorMsg01, xId:xVectorObject.id, yId:yVectorObject.id} :
{msg:errorMsg02, xId:xVectorObject.id, yId:yVectorObject.id}
} //-----------------------------------------------------------------------------------
var rowVec_Times_ColMatrix = function(xVectorObject, yMatrixObject) {
(xVectorObject.type === 'rowVector' && yMatrixObject.type === 'colMatrix') ?
(xVectorObject.ncols === yMatrixObject.nrows) ?
{id:xVectorObject.id.concat('_X_', yMatrixObject.id),
type:'rowVector',
nrows: xVectorObject.nrows,
ncols: yMatrixObject.ncols,
vector:map(function(yVectorObject) {scalarProduct(xVectorObject, yVectorObject)},
yMatrixObject.matrix)} :
{msg:errorMsg03, xId:xVectorObject.id, yId:yMatrixObject.id} :
{msg:errorMsg04,xId:xVectorObject.id, yId:yMatrixObject.id}
} //-----------------------------------------------------------------------------------
var rowMatrix_Times_ColMatrix = function(xMatrixObject, yMatrixObject) {
(xMatrixObject.type === 'rowMatrix' && yMatrixObject.type === 'colMatrix') ?
(xMatrixObject.ncols === yMatrixObject.nrows) ?
{type:'rowMatrix',
matrix: map(function(xVectorObject)
{return rowVec_Times_ColMatrix(xVectorObject, yMatrixObject)},
xMatrixObject.matrix)} :
{msg:errorMsg05, xId:xMatrixObject.id, yId:yMatrixObject.id} :
{msg:errorMsg06, xId:xMatrixObject.id, yId:yMatrixObject.id}
}
//-----------------------------------------------------------------------------------
// function body
//-----------------------------------------------------------------------------------
if(xObject.type === 'rowVector' && yObject.type === 'colVector')
{return scalarProduct(xObject, yObject)}
else if (xObject.type === 'rowVector' && yObject.type === 'colMatrix')
{return rowVec_Times_ColMatrix(xObject, yObject)}
else if (xObject.type === 'rowMatrix' && yObject.type === 'colMatrix')
{return rowMatrix_Times_ColMatrix(xObject, yObject)}
else
{return {msg:errorMsg07, xObject:xObject.xId, yObject:yObject.yId}}
}
console.log("==============================================================================")
var a1_rv = makeVector('a(1)_rv', 'rowVector', [1,-2, 0]); print(a1_rv)
var a2_rv = makeVector('a(2)_rv', 'rowVector', [0, 1, 1]); print(a2_rv)
var b1_cv = makeVector('b(1)_cv', 'colVector', [2, 1, 0]); print(b1_cv)
var b2_cv = makeVector('b(2)_cv', 'colVector', [1,-1, 0]); print(b2_cv)
var b3_cv = makeVector('b(3)_cv', 'colVector', [1, 0, 1]); print(b3_cv)
var b1_rv = makeVector('b(1)_rv', 'rowVector', [2, 1, 1]); print(b1_rv)
var b2_rv = makeVector('b(2)_rv', 'rowVector', [1,-1, 0]); print(b2_rv)
var b3_rv = makeVector('b(3)_rv', 'rowVector', [0, 0, 1]); print(b3_rv)
var d_rv = makeVector('d_rv', 'rowVector', [1, 0, 1, 2]); print(d_rv)
var e_cv = makeVector('e_cv', 'colVector', [1, 0, 1, 2]); print(e_cv)
var A_rm = makeMatrix('A_rm', 'rowMatrix', [a1_rv, a2_rv]); print(A_rm)
var B_cm = makeMatrix('B_cm', 'colMatrix', [b1_rv, b2_rv, b3_rv]); print(B_cm)
var C_rm = makeMatrix('C_rm', 'rowMatrix', [d_rv, d_rv]); print(C_rm)
console.log("------------------------------------------------------------------------------")
console.log("matrixMultiplication(a'(1), a'(2)) --->")
print(matrixMultiplication(a1_rv, a2_rv)) // types !
console.log("matrixMultiplication(a'(1), e_cv) --->")
print(matrixMultiplication(a1_rv, e_cv)) //lengths
console.log("------------------------------------------------------")
console.log("matrixMultiplication(a'(1), b(2)) === 3 --->");
console.log(matrixMultiplication(a1_rv, b2_cv) === 3)
console.log("matrixMultiplication(a'(2), b(2)) === -1 --->");
console.log(matrixMultiplication(a2_rv, b2_cv) === -1)
console.log("------------------------------------------------------")
console.log("matrixMultiplication(a'(1), A) --->")
print(matrixMultiplication(a1_rv, A_rm))
console.log("matrixMultiplication(d', B) --->")
print(matrixMultiplication(d_rv, B_cm))
console.log("------------------------------------------------------")
console.log("matrixMultiplication(a'(1), B) --->")
print(matrixMultiplication(a1_rv, B_cm)) //ok
console.log("matrixMultiplication(a'(2), B) --->")
print(matrixMultiplication(a2_rv, B_cm)) //ok
console.log("------------------------------------------------------------------------------")
console.log("matrixMultiplication(B, A) ---> "); print(matrixMultiplication(B_cm, A_rm))
console.log("matrixMultiplication(C, B) ---> "); print(matrixMultiplication(C_rm, B_cm))
console.log("matrixMultiplication(A, B) ---> "); print(matrixMultiplication(A_rm, B_cm))
var matrixC = matrixMultiplication(A_rm, B_cm).matrix
console.log("C = A_Times_B --->");
console.log("[")
print(matrixC[0].vector); print(matrixC[1].vector)
console.log("]")
console.log("------------------------------------------------------------------------------")
runTime(startTime)
console.log("==============================================================================")
Syntax: Output of Functional Vector and Matrix Multiplication in Core WebPPL
==============================================================================
file: PCM20180327_FunctionalMatrixMultiplication *** PCM *** *2018/04/05*
- definition is from: Kowalsky & Michler, Lineare Algebra, Berlin &
Walter de Gruyter, 2003, 12/e, p.48
- example is from: Kowalsky & Michler, Lineare Algebra, 2003, p.48
- further inspiration in: Abelson, Sussman & Sussman, Structure and
Interpretation of Computer Programs, ch. 2.2.3,
'Sequences as Conventional Interfaces', MIT Press,
1996, 2/e, p.113 - 126
- application is in: Goodfellow, Bengio & Courville, Deep Learning,
Cambridge, Mass.: MIT Press, 2016, p.32
(3)
==============================================================================
{"id":"a(1)_rv","type":"rowVector","nrows":1,"ncols":3,"vector":[1,-2,0]}
{"id":"a(2)_rv","type":"rowVector","nrows":1,"ncols":3,"vector":[0,1,1]}
{"id":"b(1)_cv","type":"colVector","nrows":3,"ncols":1,"vector":[2,1,0]}
{"id":"b(2)_cv","type":"colVector","nrows":3,"ncols":1,"vector":[1,-1,0]}
{"id":"b(3)_cv","type":"colVector","nrows":3,"ncols":1,"vector":[1,0,1]}
{"id":"b(1)_rv","type":"rowVector","nrows":1,"ncols":3,"vector":[2,1,1]}
{"id":"b(2)_rv","type":"rowVector","nrows":1,"ncols":3,"vector":[1,-1,0]}
{"id":"b(3)_rv","type":"rowVector","nrows":1,"ncols":3,"vector":[0,0,1]}
{"id":"d_rv","type":"rowVector","nrows":1,"ncols":4,"vector":[1,0,1,2]}
{"id":"e_cv","type":"colVector","nrows":4,"ncols":1,"vector":[1,0,1,2]}
{"id":"A_rm","type":"rowMatrix","nrows":2,"ncols":3,"matrix":[{"id":"a(1)_rv","type":"rowVector","nrows":1,"ncols":3,"vector":[1,-2,0]},{"id":"a(2)_rv","type":"rowVector","nrows":1,"ncols":3,"vector":[0,1,1]}]}
{"id":"B_cm","type":"colMatrix","nrows":3,"ncols":3,"matrix":[{"id":"B_cm0","type":"colVector","nrows":3,"ncols":1,"vector":[2,1,0]},{"id":"B_cm1","type":"colVector","nrows":3,"ncols":1,"vector":[1,-1,0]},{"id":"B_cm2","type":"colVector","nrows":3,"ncols":1,"vector":[1,0,1]}]}
{"id":"C_rm","type":"rowMatrix","nrows":2,"ncols":4,"matrix":[{"id":"d_rv","type":"rowVector","nrows":1,"ncols":4,"vector":[1,0,1,2]},{"id":"d_rv","type":"rowVector","nrows":1,"ncols":4,"vector":[1,0,1,2]}]}
------------------------------------------------------------------------------
matrixMultiplication(a'(1), a'(2)) --->
{"msg":"errorMsg07: objects have incompatible types !"}
matrixMultiplication(a'(1), e_cv) --->
{"msg":"errorMsg01: vectors have incompatible dimensions !","xId":"a(1)_rv","yId":"e_cv"}
------------------------------------------------------
matrixMultiplication(a'(1), b(2)) === 3 --->
true
matrixMultiplication(a'(2), b(2)) === -1 --->
true
------------------------------------------------------
matrixMultiplication(a'(1), A) --->
{"msg":"errorMsg07: objects have incompatible types !"}
matrixMultiplication(d', B) --->
{"msg":"errorMsg03: vector and matrix have incompatible dimensions !","xId":"d_rv","yId":"B_cm"}
------------------------------------------------------
matrixMultiplication(a'(1), B) --->
{"id":"a(1)_rv_X_B_cm","type":"rowVector","nrows":1,"ncols":3,"vector":[0,3,1]}
matrixMultiplication(a'(2), B) --->
{"id":"a(2)_rv_X_B_cm","type":"rowVector","nrows":1,"ncols":3,"vector":[1,-1,1]}
------------------------------------------------------------------------------
matrixMultiplication(B, A) --->
{"msg":"errorMsg07: objects have incompatible types !"}
matrixMultiplication(C, B) --->
{"msg":"errorMsg05: matrices have incompatible dimensions !","xId":"C_rm","yId":"B_cm"}
matrixMultiplication(A, B) --->
{"type":"rowMatrix","matrix":[{"id":"a(1)_rv_X_B_cm","type":"rowVector","nrows":1,"ncols":3,"vector":[0,3,1]},{"id":"a(2)_rv_X_B_cm","type":"rowVector","nrows":1,"ncols":3,"vector":[1,-1,1]}]}
C = A_Times_B --->
[
[0,3,1]
[1,-1,1]
]
------------------------------------------------------------------------------
runTimeInMsec ---> 267
runTimeInSec ---> 0.267
runTimeInMin ---> 0.00445
==============================================================================
Syntax: Functional Vector and Matrix Multiplication in WebPPL Using Package 'ad.tensor'
console.log("==============================================================================")
console.log("file: PCM20180405_MatrixMultiplicationWithAdTensor *** PCM *** *2018/04/05* ")
console.log(" - definition is from: Kowalsky & Michler, Lineare Algebra, Berlin & ")
console.log(" Walter de Gruyter, 2003, 12/e, p.48 ")
console.log(" - example is from: Kowalsky & Michler, Lineare Algebra, 2003, p.48 ")
console.log("==============================================================================")
//-------------------------------------------------------------------------------------------
var runTime = function(startTime) {
var stopTime = Date.now()
var runTimeInMsec = (stopTime - startTime)
var runTimeInSec = runTimeInMsec/1000.
var runTimeInMin = runTimeInSec/60.
console.log("runTimeInMsec ---> ", runTimeInMsec)
console.log("runTimeInSec ---> ", runTimeInSec)
console.log("runTimeInMin ---> ", runTimeInMin)
}
//-------------------------------------------------------------------------------------------
var startTime = Date.now()
var a1_rowVector = [ 1,-2, 0]
var a2_rowVector = [ 0, 1, 1]
var A_rowMatrix = [a1_rowVector, a2_rowVector]
console.log("A_rowMatrix:"); print(A_rowMatrix)
var b1_rowVector = [ 2, 1, 1]
var b2_rowVector = [ 1,-1, 0]
var b3_rowVector = [ 0, 0, 1]
var B_rowMatrix = [b1_rowVector, b2_rowVector, b2_rowVector]
console.log("B_rowMatrix:"); print(B_rowMatrix)
var A_linearMatrix = Matrix(A_rowMatrix)
console.log("A_linearMatrix:"); print(A_linearMatrix)
var B_linearMatrix = Matrix(B_rowMatrix)
console.log("B_linearMatrix:"); print(B_linearMatrix)
var C_linearMatrix = ad.tensor.dot(A_linearMatrix, B_linearMatrix)
console.log("C_linearMatrix:"); print(C_linearMatrix)
var C_rowMatrix = C_linearMatrix.toArray()
console.log("C_rowMatrix = A x B:"); print(C_rowMatrix)
console.log("==============================================================================")
Syntax: Output of Functional Vector and Matrix Multiplication in WebPPL Using Package 'ad.tensor'
==============================================================================
file: PCM20180405_MatrixMultiplicationWithAdTensor *** PCM *** *2018/04/05*
- definition is from: Kowalsky & Michler, Lineare Algebra, Berlin &
Walter de Gruyter, 2003, 12/e, p.48
- example is from: Kowalsky & Michler, Lineare Algebra, 2003, p.48
==============================================================================
A_rowMatrix:
[[1,-2,0],[0,1,1]]
B_rowMatrix:
[[2,1,1],[1,-1,0],[1,-1,0]]
A_linearMatrix:
{"dims":[2,3],"length":6,"data":{"0":1,"1":-2,"2":0,"3":0,"4":1,"5":1}}
B_linearMatrix:
{"dims":[3,3],"length":9,"data":{"0":2,"1":1,"2":1,"3":1,"4":-1,"5":0,"6":1,"7":-1,"8":0}}
C_linearMatrix:
{"dims":[2,3],"length":6,"data":{"0":0,"1":3,"2":1,"3":2,"4":-2,"5":0}}
C_rowMatrix = A x B:
[[0,3,1],[2,-2,0]]
==============================================================================
Geometric Distribution
Random Experiment:
The idealized random experiment which is the basis of the geometric distribution consists of a k of i.i.d. Bernoulli trials. In each trial a coin with P(success) = θ is tossed. Tossing is stopped when a 'success' is observed. Otherwise in the case of 'failure' k-1 coin tosses are repeated. X is defined as the total number of coin tosses necessary to obtain one success.
Definition:
A random variable X is said to be a geometric random variable with parameter θ, shown as
$$X \sim Geometric(θ) ,$$
if its PMF (= probability mass function) is given by
$$P_X(k\:|\:\theta) = P(X = k\:|\:\theta) = \theta(1-\theta)^{k-1}$$
$$\text{where}\:\:\:0<\theta<1\:\:\text{ and }\:\:\Omega=\left\{k\:|\:k=1,2,3,...\right\}.$$
On the next page we demonstrate with exact WebPPL functions and approximate WebPPL sampling methods how to obtain the geometric PMF (probability mass function). Functions and methods are in different styles:
- exact, nonrecursive
- exact, recursive
- exact, tailrecursive
- approximate, recursive
- approximate, tailrecursive
Exact functions compute single probabilities P(X = k) without sampling from the model function.
Approximate methods generate a single sample from the model function. This is called a geometric random trial of an idealized mathematical geometric random experiment. If this experiment is repeated with many times ("replications, samples") we get a population of samples which is the experimental basis of the abstract mathematical distribution. In our examples whe chose a small number of repetitions to demonstrate the difference of the approximations to the exact distributions. The deviations get smaller with increased number of repetitions. This can be controlled by computing the so called Hoeffding bounds (visited 2018/04/11).
References:
Lunn, D., Jackson, Chr., Best, N., Thomas, A. & Spiegelhalter, D., The BUGS Book: A Practical Introduction to Bayesian Analysis, Boca Raton, FL: CRC Press, 2013, ISBN 978-1-58488-849-9
Pishro-Nik, H., Introduction to Probability: Statistics and Random Process, Kappa Research, LLC, 2014, ISBN 978-0-9906372-0-2, www.probabilitycourse.com/ (visited 2018/04/11)
Syntax: Geometric Distribution Generated by Exact and Approximative Functions
console.log("==============================================================================")
console.log("file: PCM20180407_myGeometricDistribution *** PCM *** *2018/04/10* ")
console.log(" - definition is from: H. Pishro-Nik, Introduction to Probability: ")
console.log(" Statistics and Random Processes, ")
console.log(" Kappa Research, LLC, 2014, ISBN 978-0-9906372-0-2")
console.log(" www.probabilitycourse.com/chapter3/3_1_5_special_discrete_distr.php")
console.log(" - example is from: H. Pishro-Nik, Introduction to Probability, p.116f")
console.log("==============================================================================")
var kMax = 25 // upper limit of Bernoulli trials
//-------------------------------------------------------------------------------------------
var runTime = function(startTime) {
var stopTime = Date.now()
var runTimeInMsec = (stopTime - startTime)
var runTimeInSec = runTimeInMsec/1000.
var runTimeInMin = runTimeInSec/60.
console.log("runTimeInMsec ---> ", runTimeInMsec)
console.log("runTimeInSec ---> ", runTimeInSec)
console.log("runTimeInMin ---> ", runTimeInMin)
}
//-------------------------------------------------------------------------------------------
var myGeometricDistribution_1 = function(p) { // exact nonrecursive nonsampling function
var inc = function(x) {x + 1} // increment x
var ks = mapN(inc, kMax) // => ks = [1, 2, 3, ..., kMax]
var probs = map(function (k) { // explicit formula P(X(k)) = ...
Math.pow((1-p),(k-1)) * p // ... [(1-p)^(k-1)] * p
}, ks)
Categorical({ps: probs, vs: ks})
}
//-------------------------------------------------------------------------------------------
var myGeometricDistribution_2 = function(p) { // exact recursive function
var inc = function(x) {x + 1} // increment x
var ks = mapN(inc, kMax) // => ks = [1, 2, 3, ..., kMax]
var geometricRecursive = function(k) {
if (k === 0) {p} // success !
else {(1 - p) * geometricRecursive(k - 1)} // one more failure
}
var probs = map(geometricRecursive, ks)
Categorical({ps: probs, vs: ks})
}
//-------------------------------------------------------------------------------------------
var myGeometricDistribution_3 = function(p) { // exact tailrecursive sampling function
var inc = function(x) {x + 1} // increment x
var ks = mapN(inc, kMax) // => ks = [1, 2, 3, ..., kMax]
var geometricIter = function(px, k) { // tailrecursive helper function
if(k === 0) {px * p} // success!
else {geometricIter((1 - p) * px, k - 1)} // one more failure
}
var probs = map(function(k) {geometricIter(1, k)}, ks) // call of tailrecursive helper
Categorical({ps: probs, vs: ks})
}
//-------------------------------------------------------------------------------------------
var myGeometricDistribution_4 = function(p) { // approximative recursive sampling function
if (flip(p) === true) {1} // success !
else {1 + myGeometricDistribution_4(p)} // one more failure
}
//-------------------------------------------------------------------------------------------
var myGeometricDistribution_5 = function(p) { // approx. tailrecursive sampling function
var geometricIter = function(nFailures) { // tailrecursive helper function
var nSuccesses = 1
if(flip(p) === true) {nFailures + nSuccesses}// success!
else {geometricIter(nFailures + 1)} // one more failure
}
geometricIter(0) // call of tailrecursive helper
}
//===========================================================================================
var startTime = Date.now()
var successP = 0.3 // success probability
console.log("success probability p =", successP)
var nSamples = 500 // #number of model samples
console.log("#model samples =", nSamples)
var xLabel = {xLabel:'k (= success trial after k-1 failures)'}
console.log("==============================================================================")
console.log("myGeometricDistribution_1")
console.log("-------------------------")
console.log("exact distribution generated by nonrecursive nonsampling function for p = 0.3:")
var geometricDistribution_1 = myGeometricDistribution_1(successP)
viz(geometricDistribution_1, xLabel)
console.log("------------------------------------------------------------------------------")
console.log("myGeometricDistribution_2")
console.log("-------------------------")
console.log("exact distribution generated by recursive nonsampling function for p = 0.3:")
var geometricDistribution_2 = myGeometricDistribution_2(successP)
viz.auto(geometricDistribution_2, xLabel)
console.log("------------------------------------------------------------------------------")
console.log("myGeometricDistribution_3")
console.log("-------------------------")
console.log("exact distribution gener. by tailrecursive nonsampling function for p = 0.3:")
var geometricDistribution_3 = myGeometricDistribution_3(successP)
viz.auto(geometricDistribution_3, xLabel)
console.log("------------------------------------------------------------------------------")
console.log("myGeometricDistribution_4")
console.log("-------------------------")
console.log("approx. distribution generated by recursive sampling function for p = 0.3:")
var drawOneSampleFromModel_4 = function() {myGeometricDistribution_4(successP)}
var geometricDistribution_4 = repeat(nSamples,drawOneSampleFromModel_4)
viz.auto(geometricDistribution_4, xLabel)
console.log("------------------------------------------------------------------------------")
console.log("myGeometricDistribution_5")
console.log("-------------------------")
console.log("approx. distribution generated by tailrecursive sampling function for p = 0.3:")
var drawOneSampleFromModel_5 = function() {myGeometricDistribution_5(successP)}
var geometricDistribution_5 = repeat(nSamples,drawOneSampleFromModel_5)
viz.auto(geometricDistribution_5, xLabel)
console.log("------------------------------------------------------------------------------")
runTime(startTime)
console.log("==============================================================================")
Syntax: Output of Geometric Distribution Generated by Exact and Approximative Functions
==============================================================================
file: PCM20180407_myGeometricDistribution *** PCM *** *2018/04/10*
- definition is from: H. Pishro-Nik, Introduction to Probability:
Statistics and Random Processes,
Kappa Research, LLC, 2014, ISBN 978-0-9906372-0-2
https://www.probabilitycourse.com/chapter3/3_1_5_special_discrete_distr.php
- example is from: H. Pishro-Nik, Introduction to Probability, p.116f
==============================================================================
success probability p = 0.3
#model samples = 500
==============================================================================
myGeometricDistribution_1
-------------------------
exact distribution generated by nonrecursive nonsampling function for p = 0.3:
==========
see Fig01 below
===========
------------------------------------------------------------------------------
myGeometricDistribution_2
-------------------------
exact distribution generated by recursive nonsampling function for p = 0.3:
==========
see Fig02 below
===========
------------------------------------------------------------------------------
myGeometricDistribution_3
-------------------------
exact distribution gener. by tailrecursive nonsampling function for p = 0.3:
==========
see Fig03 below
===========
------------------------------------------------------------------------------
myGeometricDistribution_4
-------------------------
approx. distribution generated by recursive sampling function for p = 0.3:
==========
see Fig04 below
===========
------------------------------------------------------------------------------
myGeometricDistribution_5
-------------------------
approx. distribution generated by tailrecursive sampling function for p = 0.3:
==========
see Fig05 below
===========
------------------------------------------------------------------------------
runTimeInMsec ---> 240
runTimeInSec ---> 0.24
runTimeInMin ---> 0.004
==============================================================================