Help! Which language is it?
David Phillip Oster
"random variable" was a first-class data type. My instructor was quick
to point out that a random variable isn't a variable, it is a fuzzy
number. Numbers in this language were represented by curves that were
the probability distributions of the underlying number. For example,
integers were infinitely sharp spikes, with all their energy in one spot
(for example: 2). Floats (2.1) were fairly sharp spikes that accurately
represented the uncertainty in the underlying floating point
representation. The random number generator was just a
predefined-constant with a flat distribution: it is equally probabable
to be any number in the range of the underlying number representation.
The idea was to use it for real-world modelling. You just enter in what
you know, including the uncertainty on all the measured values (for
example, if you measured with a ruler, your numbers are only good to
plus or minus a thirtieth of an inch.), and the language would just
directly calculate the answer by adding and subtracting, multiplying and
dividing random variables. Than answer would itself be a random
variable, which it would print as a number with error bars, or as a 2-d
graph, so you could see how the errors propagated through your formula
and how uncertain the answer was.
Instead of running a monte carlo simulation, casting "dice" millions of
time, it would just directly compute the result.
Under the hood, the language implementation used a constraint
propogation network, and represented a random variable as a curve, and
represented that curve as a polynomial, and represented that polynomial
as a variable length array of coefficients plus an error closure. It
would carry out calculations only to as many digits as necessary to
produce the format of the answer, and generate more coefficients to make
the curves fit the data better, as required. Adding and subtracting,
multiplying and dividing was done the way a symbolic math package would
do it: by applying rules to lists of coefficients.
Does anyone know which language this was? Where can I find out more
about it? and: What has been done with these ideas in the last 20 years?
Ken Shan
It is described in
Daphne Koller, David McAllester and Avi Pfeffer (1997).
Effective Bayesian Inference for Stochastic Programs.
In AAAI-97: Proceedings of the 14th National Conference on
Artificial Intelligence, Providence, RI.
and also in Avi Pfeffer's thesis. Both of these are available online
at http://www.eecs.harvard.edu/~avi/Research/prob.html
If you don't particularly care about numeric *probabilities*, but do
want to keep track of multiple *possibilities*, your interests appear
related to those of Jan van Eijck (http://www.cwi.nl/~jve/dynamo/).
There are probably other related projects, but these two are the ones
I am more familiar with.
--
Edit this signature at http://www.digitas.harvard.edu/cgi-bin/ken/sig
A donut a day keeps the doctor away.
Brian Webb
"David Phillip Oster" <os...@ieee.org> wrote in message
news:oster-C37BB7....@news.sf.sbcglobal.net...
I'm not familiar with the language you're describing, but the techniques
involved sound like "interval analysis". With interval analysis, values
are represented as a range of values, like X = [0.9,1.0]. Operations on
these values propagate errors. For instance [1,2] + [2,4] = [3,6], where
the result covers every possible value of the operands.
These techniques are being used for real-world modeling to implement
global optimizers and for verifying calculations by supplying error bounds
with each operation. To learn more about it, go to the following site:
http://www.cs.utep.edu/interval-comp/
Sun has built-in support for interval variables in it's C++ and FORTRAN
compilers. GNU FORTRAN also has interval types.
An online tool called UniCalc uses intervals and constraint propagation
to implement a global nonlinear equation solver. You can see it at:
http://www.rriai.org.ru/UniCalc/
- Brian
David Phillip Oster
"Brian Webb" <iRo...@swbell.net> wrote:
> I'm not familiar with the language you're describing, but the techniques
> involved sound like "interval analysis". With interval analysis, values
> are represented as a range of values, like X = [0.9,1.0]. Operations on
> these values propagate errors. For instance [1,2] + [2,4] = [3,6], where
> the result covers every possible value of the operands.
Yes, it does sound like interval analysis, but interval analysis is the
simple cousin of the language I'm looking for. In interval analysis, the
assumption is that the actual value has an equal probability of being
anywhere in the interval: always a flat probability curve. In the
language I'd heard about, you could specify the exact shape of the
probability distribution curve.
Thanks for the references, I will check them out.
-- David
Jón Fairbairn
> Yes, it does sound like interval analysis, but interval analysis is the
> probability distribution curve.
This is surely something that one could achieve fairly easily in a
strongly-typed functional language like Haskell? Represent the
number-analogues using functions and compute with those.
That might be a bit slow, but it would get the principle down. There
are probably more efficient ways of representing things, and laziness
might come in useful.
--
Jón Fairbairn Jon.Fa...@cl.cam.ac.uk
George Russell
errors caused by floating-point rounding (rather than measurement error), you
should read
http://HTTP.CS.Berkeley.EDU/~wkahan/improber.pdf
This would hopefully dissuade you from modelling the error distribution as a
smooth curve . . .
Jerzy Karczmarczuk
numbers" :
I am not sure about that last point.
Actually I did it, I implemented a long time ago an arithmetic package
like that, just for fun.
Let's forget about annoying problems I had in order to lift arithmetic
(+,*, etc.) to functional objects in Haskell, enough to say that Num
required Eq, and comparing function was, hm... not so kosher.
If f,g are prob. distributions, you know what you have to do in order to
compute f+g? You have to compute the convolution integral of f and g.
Multiplication is even worse. Division? If the divisor support includes zero,
the result blows up, you end up with an infinite interval. Plenty of rather
costly numerics. Laziness won't help you. It is easier with constant intervals,
and with the conservative approach to interval combination, you take unions...
That's why from the practical point of view, interval arithmetic is not so
popular. The intervals grow up rather fast at every operation performed, and
a few operations suffice to produce results not very significant.
(I won't be categorical; my friends hate intervals, but the organizers of this
conference:
http://eia.udg.es/~misc99/course.html
will tell you that intervals will save the world, and that critics are simply
not too competent.)
You can do better, of course. For example postulate that all distributions
*are* gaussian, and construct the best fit to a gaussian even if you divide
those distributions, neglecting the fact that in principle you might divide
by zero. Then the variances add, the standard deviations which are square roots
of variances grow more slowly.
I can't say much about the usability of such a system, I hadn't even the courage
to propose it to a workshop, although I gave the stuff to my physicist friends,
suggesting that it might be useful for an automatized error analysis...
And, anyway, I found out that delaying lazily
all those convolutions, etc. does no good, it just clogs the memory with
thunks.
==
Oh. George Russel proposes the article by Kahan, I wanted to cite it, and also
a talk given at the last PPDP, but all that concerns rounding errors rather
than a priori fuzzy numbers, with indeterminacy much bigger than the
rounding bits.
There are many articles about probabilistic intervals. The last I found is here:
http://space.iias.spb.su/ai/doc/IntervalProb-98-Rus.zip
the only (yours!) problem is that you should read Cyrillic...
Jerzy Karczmarczuk
Caen, France
Alex Wu
news:3BDFE855...@info.unicaen.fr...
> I can't say much about the usability of such a system, I hadn't even the
courage
> to propose it to a workshop, although I gave the stuff to my physicist
friends,
> suggesting that it might be useful for an automatized error analysis...
I was part of an industrial risk analysis project over five years ago that
used a lot of this -- however, I did the whole thing in obfuscated C++,
i.e., C ;-). However, at the time I don't remember a widely available
language that provided a "natural" way of approaching the problem ...
Alex
Simon Slavin
David Phillip Oster <os...@ieee.org> wrote:
> In 1978 I heard about a Phd thesis describing a computer language where
> "random variable" was a first-class data type.
Try posting that again to alt.folklore.computers. But any old
PhD can invent a computing language and many do. For all I
know only the author and his supervisor might recognise this
language. Also, a facility like the variables you mention
could be added to any existing constraint-language like PROLOG.
Simon.
--
http://www.hearsay.demon.co.uk | I have a hunch that [] the unknown sequences
No junk email please. | of DNA [will decode into] copyright notices
| and patent protections. -- Donald E. Knuth
The French Was There.
ol...@pobox.com
> > are probably more efficient ways of representing things, and laziness
> > might come in useful.
> I am not sure about that last point.
possible. Monads and lazy evaluation come in very handy. This article
shows the example.
The key insight is that we don't need to "symbolically" manipulate
distributions. We represent a random variable by a function, which,
when evaluated, generates one number. If evaluated many times, the
function generates many numbers -- a sample of a given distribution.
We can combine such functions into expression. At the very end, we
evaluate that expression many times -- which gives us the final
sample. From that sample, we can compute the estimate of the
mathematical expectation, and the estimate of the variance (and of
higher moments if needed). This is basically a Monte Carlo
approach. True, Willem Kahan doesn't like it -- when it is applied to
numerical instabilities that is. In physics, this is a highly
justified approach. Most (all?) of our experimental results are
ensemble averages.
As an example, we compute the center and the variance of
sqrt( (Lw - 1/Cw)^2 + R^2)
which is the impedance of an electric circuit, with an inductance
being a normal variable (center 100, variance 1), frequency being
uniformly distributed from 10 through 20 kHz, resistance is also
uniformly distributed from 10 through 50 kOm and the capacitance has a
square root of the normal distribution. I admit, the latter is just to
make the example more interesting.
The framework first:
import Monad
-- An integer uniformly distributed within [0..32767]
type RandomState = Int
random_state_max::RandomState
random_state_max = 32767
type RandomGen = RandomState -> RandomState
newtype RandomVar a = RandomVar ( RandomGen -> RandomState -> (a,RandomState) )
instance Monad RandomVar where
return x = RandomVar $ \gen state -> (x,state) -- deterministic variable
(RandomVar x) >>= f = RandomVar $
\gen state -> let (v,state') = x gen state
RandomVar x' = f v
in x' gen state'
-- Other morphisms and monad constructors
-- A random number uniformly distributed within [0,1]
uniform_stan:: RandomVar Float
uniform_stan = RandomVar $
\gen state ->
let state' = gen state
in ((fromInt state') / (fromInt random_state_max), state')
-- A random number uniformly distributed on [a, b]
uniform a b = do { x <- uniform_stan; return $ a+x*(b-a) }
add:: (Num a) => (RandomVar a) -> (RandomVar a) -> (RandomVar a)
add = liftM2 (+)
sub:: (Num a) => (RandomVar a) -> (RandomVar a) -> (RandomVar a)
sub = liftM2 (-)
mul:: (Num a) => (RandomVar a) -> (RandomVar a) -> (RandomVar a)
mul = liftM2 (*)
divM:: (Num a,Fractional a) => (RandomVar a) -> (RandomVar a) -> (RandomVar a)
divM = liftM2 (/)
sqrM:: (Num a) => (RandomVar a) -> (RandomVar a)
sqrM = liftM $ \x -> x*x
sqrtM:: (Num a, Floating a) => (RandomVar a) -> (RandomVar a)
sqrtM = liftM sqrt
-- standardized normal distribution N(0,1)
-- This is a well-known binomial approximation
normal_stan = sub (foldM (\acc _ -> add (return acc) uniform_stan) 0 [1..12])
(return 6)
-- N(z,var)
normal z var = do { x <- normal_stan; return $ x*var + z }
-- square_root distribution
sqrt_dist dist = do { x <- dist; y <- dist; return $ max x y }
-- Run the RandomVar monad n times and estimate mathematical expectation and
-- the variance
run_random:: Int -> (RandomVar Float) -> (Float,Float)
run_random n (RandomVar mx) =
let state = 5 -- initial state
gen::RandomGen -- the pathetic random number generator
gen state = ((state * 1103515245) + 12345) `mod` (random_state_max+1)
samples = take n $ iterate (\(_,state) -> mx gen state) $ mx gen state
(m0,m1) = foldl (\(acc_0,acc_1) (v,_) -> (acc_0+v,acc_1+v*v))
(0,0) samples
in (m0/(fromInt n),sqrt( (m1 - m0*m0/(fromInt n))/(fromInt (n-1)) ))
Now let's see what we've got
Main> run_random 5000 uniform_stan
(0.49526,0.286909)
well, close enough to the expected result.
Main> run_random 5000 normal_stan
(-0.00557715,0.997899)
Again, normal_stan indeed looks like a normally-distributed variable
with center 0 and variance 1.
We can now do many wonderful things, such as
Main> run_random 1000 $ add (normal 10 1) (normal 20 2)
(29.9166,2.17711)
Back to the impedance however
-- compute the impedance of an electric circuit
-- z = sqrt( (Lw - 1/Cw)^2 + R^2)
impedance inductance capacitance freq resistance =
sqrtM $
(sqrM
(freq `mul` inductance) `sub` ((return 1) `divM` (capacitance `mul` freq)))
`add`
(sqrM resistance)
-- inductance
ct_l = normal 100 1
-- Frequency 10kHz - 20kHz
ct_freq = uniform 10000.0 20000.0
-- capacitance
ct_cap = sqrt_dist $ normal 1000 10
-- resistance
ct_res = uniform 10000.0 50000.0
z = impedance ct_l ct_cap ct_freq ct_res
Main> run_random 1000 z
(1.51302e+06,286369.0)
Here you have it: the average value and the variance. We can repeat
the sample with bigger n and watch the convergence -- just as in the
regular Monte Carlo method.
George Russell
it just produces upper and lower bounds. Unfortunately for non-trivial calculations
the upper and lower bounds are nearly always dreadfully pessimistic, which is why
interval-analysis algorithms frequently work by computing the answer non-rigorously, and
then finding some cunning technique with not much computation to bound it.
If you want to compute probabilities, this sounds horrendously complicated, because you
have to take into account all the potential interactions of the variables. So if you
compute
X = Y + Z
to compute the probability distribution for X, you need not just the probability distributions
for Y and Z, but their joint probability distribution.
The original post suggested representing the distribution as a polynomial, but how can
this be? The distribution must go to zero as you go to infinity in either direction,
but all non-trivial polynomials go to + or - infinity. But you might be able to do
something of this sort using something like wavelets. However I'm sceptical of the
merits of this approach. If you just want to compute the sensitivity of your output
to error in the original input, I would have thought just computing the partial derivatives
would be a lot easier and a lot cheaper.
Jerzy Karczmarczuk
1. Somebody wants the references concerning a language which manipulates
"fuzzy", probabilistic numbers.
2. Plenty of reaction. Jón Fairbarn suggests using functions, and says that
laziness may be useful.
3. I disagree (mildly) on the last point, since my own struggle with numerical
applications of lazy languages convinced me that if laziness in this
business can be avoided, then it should. (It doesn't change the fact that
I love and use lazy structures quite often).
4. Oleg Kiselyov is more positive on this:
> Actually, that _can_ be done. It is slow, unwieldy -- but eminently
> possible. Monads and lazy evaluation come in very handy. This article
> shows the example.
>
> The key insight is that we don't need to "symbolically" manipulate
> distributions. We represent a random variable by a function, which,
> when evaluated, generates one number. If evaluated many times, the
> function generates many numbers -- a sample of a given distribution.
> We can combine such functions into expression. At the very end, we
> evaluate that expression many times -- which gives us the final
> sample. From that sample, we can compute the estimate of the
> mathematical expectation, and the estimate of the variance (and of
> higher moments if needed). This is basically a Monte Carlo
> approach. True, Willem Kahan doesn't like it -- when it is applied to
> numerical instabilities that is. In physics, this is a highly
> justified approach. Most (all?) of our experimental results are
> ensemble averages.
[...]
An example folows. Nice and inspiring. And of course *may* be used.
But...
a. It hasn't much to do with laziness, not directly anyway. All your
monads seem strict, you use also foldl.
b. Saying that in physics this is "highly justified" may be wrongly
understood, and I protest. Monte-Carlo is used in physics where
it *IS* efficient, where you cannot do your integrals otherwise,
because e.g., the space dimension is too big to give any efficiency to
regular integration schemes. This is not the case when convoluting
expressions. Integrals are 1-dim, iterated.
c. Sampling must be very thoroughly controlled, otherwise, if you
compute A/B with B range containing zero, you may divide by zero
in the middle of the process, and you will be prosecuted.
=====
My positive comment: Yes, use functions and laziness, but in a semi-
numerical framework. Develop polynomially your distributions. No, it
is not stupid, as somebody remarked, observing that they should vanish
far away. Often the initial distributions *are* compact, and you may
use Chebyshev expansions and so on.
For a physicists natural distribution is Gaussian, natürlich.
+ If you feel that your distributions are more complex than that,
find out whether they are uni- or multi-modal; combine gaussians.
If you have only one, which is typical:
+ Evaluate the skewness and curtosis, and expand the distributions
*exp-polynomially around a gaussian* (like in Laplace steepest-descent
algorithm).
+ Use those approximations to convolve the distributions. Make the
asymptotic evaluation of your convolutions, so your functions will
have always such expanded forms.
+ Automatic differentiation, especially lazy automatic differentiation
may be quite handy. (Of course this is my personal, shameless plug-in...)
Enough of this Halloween subject...
Jerzy Karczmarczuk
Caen, France