New Function Recommendations

[NunoSempere]

[From Ozzie]
Recommend new functions for Squiggle here.

Checklist

All recommendations are appreciated! However, it's extra helpful if you could include any/all of the following:

1. What is the code for? What use cases might it support?
2. How important/critical is it?
3. Do other related tools have similar functions? If so, please link to them. If the API you suggest is different from theirs, explain why. They main related tools would be Python, Julia, Mathematica, R, Rescript, Guesstimate, PyMC3.
4. What should the namespace be? (Dict.X, Duration.X, Dist.X).
5. What is the type definition? Like,SampleSet.fromDist: (list<number>) => sampleSet
6. Should we begin this as an experimental feature? (Put it in the namespace Danger.x for a while, until we're more comfortable with it)
7. How can we best implement this? (Pseudo-code or working code is useful)

Checklist template

(Copy and paste this into your comment)

## [Function Name]
```
[fnName with namespace and definitions]
```
### Description
### Importance
### Similar Implementations / Inspiration
### Should this begin experimentally?
### How to implement

Considerations

New functions are fairly costly. They complicate the code, and because it's bad for us to remove functionality (this would break people's code), we need to support them indefinitely. One exception is that we can remove them in major releases, but we shouldn't do this often.

What's particularly bad is changing function definitions to mean different things. For example, if we support random() to return a random number from 0 to 1, then change it to return a number from -Infinity to Infinity, later. This would break people's code in ways that are difficult to notice, which is basically catastrophically bad.

[NunoSempere]

choose(n,k) = factorial(n)/(factorial(n-k) * factorial(k))

[OAGr]

You can do this, but it's called combinations.

[NunoSempere]

The implementation of this doesn't work for large numbers, perhaps because it requires large factorials. Ended up using a normal approximation instead, making a note for future use.

[NunoSempere]

laplace(s,n) = (s+1)/(n+2)

[OAGr]

One issue here is to distinguish this from the laplace distribution.

[quinn-dougherty]

laplaceSucc

[NunoSempere]

binomial(k,n,p) = choose(n,k) * p^k * (1-p)^(n-k)

[OAGr]

Also needs to be distinguished from the distribution.

[OAGr]

Listing standard functions is good. Also, we should organize them into packages, like, Math.binomial.

[OAGr]

Anything here would be fairly easy to add, though each item would take a bit of work: https://stdlib.io/docs/api/latest/@stdlib/math/base/special

(really, most specific functions in stdlib that don't involve other types should be okay)

[OAGr]

Here's the library of functions we currently import, though I want to have them be manual, later on.
https://mathjs.org/docs/reference/functions.html

[NunoSempere]

integrate(fun, min, max) = {
    // assume that min and max are integers.
    epsilon = 0.001
    l = max - min
    meanF(t) = mean(fun(t))
    intervals = map(upTo(0, (l/epsilon)), ({|n| min + n*epsilon}))
    values = map(intervals, ({ |x | meanF(x)})) 
    result = reduce(values, 0, ({|acc, x| acc + x})) * epsilon
    result
}

// Example
f(t) = normal(1, 2)
integrate(f, 1, 10)

[OAGr]

Good point on the last bit. I like giving users the option at this point.

[quinn-dougherty]

We should be serving at least a few integrals. A whole namespace, even, Integrate.simpson, Nuño's code here is a little Riemanny so we could call it Integrate.riemann. Integrate.trapezoid

[NunoSempere]

L mentions using:
integrate(fun, min, max) = {
// assume that min and max are integers.
epsilon = 1
l = max - min
meanF(t) = fun(t)
intervals = map(List.upTo(0, (l/epsilon)), ({|n| min + n*epsilon}))
values = map(intervals, ({ |x | meanF(x)}))
result = reduce(values, 0, ({|acc, x| acc + x})) * epsilon
result
}

[renxida]

Oh hey the ability to basically numerically solve stochastic differential equations would be AMAZING to have; but also, how costly is this going to be?

[OAGr]

Yea; differential equation solving is it's whole can of worms, I doubt we'd get to that for a while. Really neat though.

[NunoSempere]

powIndependentSamples(dist, exponent) = {
  n = floor(exponent)
  remainder = exponent - n
  helperList = upTo(1,n)
  dist_to_the_nth_power = reduce(helperList, 1, {|acc, _| acc * toSampleSet(dist) }) 
  dist_to_the_remainder = toSampleSet(dist)^remainder
  result = dist_to_the_nth_power * dist_to_the_remainder
  result
}

// Example

growth = 1 + (-0.01 to 0.05)
powIndependentSamples(growth, 100)

[OAGr]

One thing we could do here is have alternative functions for these calls. As a bad example:
independentProduct, independentExponent.

Separately, I'm thinking of a different way of writing this as something like,

a = normal(5,2)
sampleNTimesFn(t) = {|| => sampleN(a,t) |> product}
f(t) = SampleSet.make(sampleNTimesFn(t))

Here, whenever SampleSet.make is called with a zero-argument function that returns a number, it calls it a bunch of times in order to fill the distribution.

[NunoSempere]

loguniform(a, b) = exp(uniform(log(a), log(b)))

[NunoSempere]

Surprised that this isn't in the standard library already

[OAGr]

One annoying thing is that we'd need more work to handle this symbolically; we'd want the formulas for mean/stdev/percentiles and similar. We could do it without that, but it would be unusual. The method you suggest would be fairly low fidelity; requires sampling.

[NunoSempere]

The method you suggest would be fairly low fidelity; requires sampling.

Yeah, I agree; it would be better to have it built-in.

[NunoSempere]

geometric mean of probabilities, and other methods implemented in https://www.npmjs.com/package/@forecasting/aggregation

[berekuk]

List.upTo name is annoying. In Ruby it's x.upTo(y), which makes sense, but in Squiggle it's upTo(from, to), and "up to... from!" trips me over every time.

I'd prefer it to be called range, like in Python. Though Python's ranges are a bit different, they don't include the right boundary (range(3,5) is [3,4]). And they're also lazy iterators, but that doesn't seem too important.

[OAGr]

I don't think we have simple random() functions yet.

[NunoSempere]

// Example distribution creation from an inverse cdf
// also documents loguniform.
createInverseCdf(min, max) = {|y| min*(max/min)^(y)}

inverseCdfToSampleSet(inverseCdf) = {
    seedSamples = sampleN(uniform(0,1), 10000)
    samples = map(seedSamples, inverseCdf)
    SampleSet.fromList(samples)
}

inverseCdf = createInverseCdf(1, 50)
inverseCdfToSampleSet(inverseCdf)

[OAGr]

Btw, for the future, I recommend adding the namespaces as well. All functions should have default namespaces, though some are available to use without these. (It's confusing, but I think less confusing/annoying than alternatives)

PointSet.fromInverseCdfFn(number=>number) (We can then call this from Dist.fromInverseCdfFn later)

If we have the inverse CDF, I'd prefer to just take all xy coordinates from x=0 to 1. This will be more accurate than sampleset stuff.

[NunoSempere]

random() = sample(uniform(0,1))
// random()

[OAGr]

Bonus points if you could link to implementations in Python / Julia, so we could make sure to be standard with existing work. I think some groups do random() from 0 to 1, but I'd want to check to make sure.

When deciding on a new key function for the main library, it seems important to be fairly thorough about making sure we're confident in the API. It's pretty bad to change later. Like, if we changed random() from 0 to 1, to 0 to 10, at any point, that would really mess with people's work.

[NunoSempere]

Ozzie mentions: having a Danger namespace, so that we can call Danger.newFunction, but know that it's pretty experimental.

[OAGr]

System.env()

Would return a record, like,

{
  sampleCount: 1000,
  xyPointLength: 1000
}

Optionally we could send in the other view data too, but that seems pretty unnecessary.

We could also have separate functions like System.sampleCount and System.xyPointLength. This might be better because we might want to later rename xyPointLength, and this way, we could just redirect this to the new name for a while if needed.

[renxida]

I'm extremely new to Squiggle, but I would love a suite of max-entropy distribution inference functions.

Something along the lines of

1. I know the min, max, median, and number of samples after sampling an unknown distribution. What is something that this distribution may look like?
2. I know the a, b, cth percentiles of a distribution. What is something this distribution could look like (e.g. some data sources only give you 25, 50, 75th percentiles)

[OAGr]

1. Thanks! I like the idea of this, but not sure how to best do it. I imagine it would take some time to figure out. If someone suggests a specific implementation, that would help a lot here.
2. I think the metalog distribution could help here, but we don't yet have a JS implementation.

[NunoSempere]

To do: Figure out a neat/short way to express the following type of relationship: Given a distribution (e.g., funds), combine it with a function which goes from each point of a distribution to another distribution.

E.g., you have some uncertainty about how many funds EA will have, and some uncertainty of how much impact to achieve per unit of money, but you also have diminishing returns in terms of money.

ea_capital = 15B to 100B
value_of_capital(c) = {
    utilons = to(log(c), log(c) + c/2)
    utilons
}

I can't find the function to manipulate xyShapes, but if I find it I'll give an example.

[leopoldaschenbrenner]

agreed this would be nice. So far just doing this via SampleSet.map

[OAGr]

[Edit: Actually, this specific example would be tough, because it returns a distribution.]
We'll add PointSet.mapY and PointSet.mapX sometime, they should handle this. (Would work well for point set distributions and okay for symbolic distributions, not so much for SampleSet distributions)

[OAGr]

Can you describe what you want to do with this? I think that typically Monte Carlo sampling is the right answer here. It seems tricky to me to have a pointSet "distribution" where every point is itself a distribution. I guess I could see it, but it seems complicated, I'm not sure if it's really needed.

Typically the thing to do is to collapse a "distribution of distributions" into a single distribution.

[NunoSempere]

Yes. Suppose you have Experimental.PointSetToXYShape. Then, conceptually, I want to do something like:

value_of_capital(c) = {
    utilons = to(log(c), log(c) + sqrt(c))
    utilons
}

ea_capital = 15B to 100B
ea_impact_dist_of_dists = map(Experimental.PointSetToXYShape(ea_capital), {|p| 
  x = p.x
  y = p.y
  ea_impact_if_x = value_of_capital(x)
  result = {dist: ea_impact_if_x, y: y}
})
dists_array = map(ea_impact_dist_of_dists, {|d| d.dist})
dists_weights = map(ea_impact_dist_of_dists, {|d| d.y})
ea_impact_dist = mx(dists_array, dists_weights)
ea_impact_dist

Caveat: I'm guessing that Squiggle isn't able to deal with the giant mx at the end, so I would have to think about a more efficient implementation, but that is the general idea.

In contrast, what I can currently do is:

value_of_capital = 30 to 200k
ea_capital = 15B to 100B
ea_impact = 15B to 100B / value_of_capital

But that isn't the same thing, and in particular it doesn't take into account marginal returns, which bothers me a bit.

[NunoSempere]

Thinking about this a bit more, I think I would end up using sample sets anyways, so it doesn't seem so urgent.

[NunoSempere]

truncate(d, min, max) = truncateLeft(truncateRight(d, max), min)

I remember this existing, but can't find it

[OAGr]

Yea, we had something like this in a much earlier version of squiggle, but it wasn't added to newer versions.

[OAGr]

It looks like it's still there, just not documented. I'm adding now.

[NunoSempere]

Generalized correlations:

q = uniform(0, 1) |> SampleSet.fromDist

f(t) = SampleSet.map(q, {|x| quantile(beta(10, 5 + t), x)})

f

I.e., if we have a distribution, and each point of that distribution is associated with a different distribution.

[NunoSempere]

More generally, there are a number of solutions to this:

• Chaining distributions (might involve some elegant math trick, like convolution?)
• Chaining samplesets (might be more efficiently doable if we know beforehand that this is explicitly what we are doing).

I think that there is an elegant solution here, but it would take some thought.

[NunoSempere]

Following up on this:

dist = 1 to 10 
f(point) = to(0, 1+point)
recurseDist(dist, f) = {
  distAsSampleSet = SampleSet.fromDist(dist)
  result  = SampleSet.map(distAsSampleSet, {|point|  sampleN(f(point), 1)})
  result 
}

[Xodarap]

zip(left, right) = map(List.upTo(0, List.length(left)), {|idx| [left[idx], right[idx]] })

I think I'm misunderstanding something because this doesn't work, it says length is not a function. But this is a pretty standard functional programming function; I needed in the very first squiggle program I tried to write, for what that's worth.

[berekuk]

Yes, List.length is currently missing. We'll add it in the next release, #1217.

[berekuk]

For now you can simulate it with length(list) = List.reduce(list, 0, {|acc,v|acc+1}).

Btw, your upper upTo boundary should be length - 1; List.upTo includes both endpoints.

[joel-becker]

student t distributions!

my current use case is estimating the posterior distribution of a mean.