Automatic Formula Discovery in the Wolfram Language

Giorgia Fortuna, Wolfram Research

Join the Conversation #WolframTechConf

Introduction

Discover the hidden structure in some data

What is the probability that a word beginning with “a” has length 5?

{"a","aah","aardvark","aardvarks","abaci","aback","abacus",«4486»,"azalea","azaleas","azimuth","azimuthal","azimuths","azure","azures"}

Grid[{{WordCloud[StringLength /@ DictionaryLookup["a" ~~ ___], ImageSize->{200, 200}], Histogram[StringLength /@ DictionaryLookup["a" ~~ ___], ImageSize->{{200, 200}}]}}, Frame->All]

What is the expected weight of a person 75 inches tall?

{(Quantity[58, "in", inches, "Inches"])->(Quantity[115, "lb", pounds, "Pounds"]),(Quantity[59, "in", inches, "Inches"])->(Quantity[117, "lb", pounds, "Pounds"]),(Quantity[60, "in", inches, "Inches"])->(Quantity[120, "lb", pounds, "Pounds"]),(Quantity[61, "in", inches, "Inches"])->(Quantity[123, "lb", pounds, "Pounds"]),(Quantity[62, "in", inches, "Inches"])->(Quantity[126, "lb", pounds, "Pounds"]),(Quantity[63, "in", inches, "Inches"])->(Quantity[129, "lb", pounds, "Pounds"]),(Quantity[64, "in", inches, "Inches"])->(Quantity[132, "lb", pounds, "Pounds"]),(Quantity[65, "in", inches, "Inches"])->(Quantity[135, "lb", pounds, "Pounds"]),(Quantity[66, "in", inches, "Inches"])->(Quantity[139, "lb", pounds, "Pounds"]),(Quantity[67, "in", inches, "Inches"])->(Quantity[142, "lb", pounds, "Pounds"]),(Quantity[68, "in", inches, "Inches"])->(Quantity[146, "lb", pounds, "Pounds"]),(Quantity[69, "in", inches, "Inches"])->(Quantity[150, "lb", pounds, "Pounds"]),(Quantity[70, "in", inches, "Inches"])->(Quantity[154, "lb", pounds, "Pounds"]),(Quantity[71, "in", inches, "Inches"])->(Quantity[159, "lb", pounds, "Pounds"]),(Quantity[72, "in", inches, "Inches"])->(Quantity[164, "lb", pounds, "Pounds"])}

({{, }})

Have an understandable, simple and explicit model

FindDistribution & FindFormula

FindDistribution and FindFormula are meant to find good models that fit (but not overfit) the data, which are simple, explicit, and easy to understand.

({{({{, }}), }})

FindDistribution can be used to solve the “letter a” example.

FindFormula can be used to solve the “height-weight” example.

The main difference between these two functions is that the models to analyze in the first one are finite (possible distributions) whereas there are infinitely many models (functions) for the second one.

({{, }})

FindDistribution in Practice

Text Frequency

Given a text, what is the distribution of occurrences of words in it?

the frequency of any word is inversely proportional to its rank in the frequency table.Thus the most frequent word will occur approximately twice as often as the second most frequent word

PacletDirectoryAdd["~/git/ARG/SymbolicMachineLearning"];
Get["SymbolicMachineLearning`"]

text = WikipediaData["London"];
words = TextCases[DeleteStopwords[text], "Word"];

wordCount=Tally[words][[All,2]];

wordCount

Reverse[Sort@wordCount]

dist = FindDistribution[wordCount]

Show[Histogram[wordCount,{1,40,0.9},"ProbabilityDensity"],DiscretePlot[PDF[dist,x],{x,1,40},PlotStyle->PointSize[Medium], PlotRange->All]]

Earthquake

What is the probability of not having any earthquake for the next 500 days?

({{, , , }})

ExampleData[{ "Statistics","EarthquakeWaitingTimes"}, "LongDescription"]
data = ExampleData[{ "Statistics","EarthquakeWaitingTimes"}]

distribution = FindDistribution[data]

Back to the question

1-CDF[distribution, 500]

Singers Hight

({{, , }})

ExampleData[{ "Statistics","SingerHeights"}, "LongDescription"]
data = ExampleData[{"Statistics","SingerHeights"}];

datagathered = GatherBy[data, First[Characters[#][[2]]]&][[All, All, 1]];

Histogram[Flatten@datagathered]

distributions = FindDistribution[#, TargetFunctions->"Continuous"]&/@datagathered

Grid[{MapThread[Labeled[Histogram[#1], #2 ]&, {datagathered,  {"Soprano", "Alto", "Tenor", "Bass"}}]}, Frame->All]

FindFormula in Practice

Female Heights and Weights

Quit[]

ExampleData[{ "Statistics","FemaleHeightsAndWeights"}, "ColumnDescriptions"]//Column
data = ExampleData[{ "Statistics","FemaleHeightsAndWeights"}]

ListPlot[data]

fit = FindFormula[data]

Show[Plot[fit[x], {x, 50, 72}], ListPlot[data, PlotStyle->Gray]]

What is the expected weight of a woman 80 inches tall?

fit[80]

What is the distribution of weights for soprano singers?

sopranos = {64,62,66,65,60,61,65,66,65,63,67,65,62,65,68,65,63,65,62,65,66,62,65,63,65,66,65,62,65,66,65,61,65,66,65,62,63,67,60,67,66,62,65,62,61,62,66,60,65,65,61,64,68,64,63,62,64,62,64,65,60,65,70,63,67,66};

weightss = fit[sopranos]

FindDistribution[weightss]

Look at more results

SeedRandom[989]
more = FindFormula[data, x, 4, Method->"ParallelTempering"]

Multicolumn[Table[Show[Plot[i, {x, 50, 82}], ListPlot[data, PlotStyle->Red],ImageSize->350 ], {i, more}],2, Frame->All]

Prime numbers

What is the number of primes less or equal of a given number?

The frequency of primes can be asymptotically be aproximates with n/Log[n]

primep = Table[PrimePi[n],{n,80}]

Grid[{{Labeled[ListPlot[primep, ImageSize->400], "plot of primes less or equal than ..."],    Labeled[Show[ListPlot[primep, PlotStyle->Gray], Plot[{x/Log[x]}, {x, 3, 80}], ImageSize->400], "comparison with the function x/Log[x]"]}}, Spacings->4]

fitprimes  =  FindFormula[primep]

This function is NOT  x/Log[x], but it fits the points

Show[ListPlot[primep, PlotStyle->Gray], Plot[fitprimes[x], {x, 1, 100}]]

What is the the number of primes less than 90?

fitprimes[90]

PrimePi[90]

The fit should generalize to unseen data

remprimep = Table[PrimePi[n],{n,80, 200}];
remprimep = Transpose[{Range[80, Length[remprimep]+79], remprimep}];
lm = ListPlot[{primep, remprimep}, PlotLegends->{"seen data", "unseen data"}]

Show[lm, Plot[fitprimes[x], {x, 1, 200}]]

Monthly Cow Milk over 13 years

({{, , }})

milk = {{"1962-01",578.3`},{"1962-02",609.8`},{"1962-03",628.4`},{"1962-04",665.6`},{"1962-05",713.8`},{"1962-06",707.2`},{"1962-07",628.4`},{"1962-08",588.1`},{"1962-09",576.3`},{"1962-10",566.5`},{"1962-11",561.1`},{"1962-12",571.4`},{"1963-01",589.1`},{"1963-02",615.3`},{"1963-03",641.2`},{"1963-04",682.8`},{"1963-05",728.5`},{"1963-06",726.4`},{"1963-07",648.`},{"1963-08",605.8`},{"1963-09",591.5`},{"1963-10",576.3`},{"1963-11",573.2`},{"1963-12",587.1`},{"1964-01",616.6`},{"1964-02",648.6`},{"1964-03",675.5`},{"1964-04",715.3`},{"1964-05",756.`},{"1964-06",746.7`},{"1964-07",665.7`},{"1964-08",627.4`},{"1964-09",612.8`},{"1964-10",599.9`},{"1964-11",602.7`},{"1964-12",622.5`},{"1965-01",646.1`},{"1965-02",676.1`},{"1965-03",696.1`},{"1965-04",732.5`},{"1965-05",767.8`},{"1965-06",767.`},{"1965-07",689.3`},{"1965-08",641.2`},{"1965-09",624.`},{"1965-10",609.7`},{"1965-11",610.8`},{"1965-12",623.5`},{"1966-01",664.7`},{"1966-02",690.3`},{"1966-03",722.6`},{"1966-04",766.`},{"1966-05",796.3`},{"1966-06",809.6`},{"1966-07",721.7`},{"1966-08",684.4`},{"1966-09",670.6`},{"1966-10",654.9`},{"1966-11",654.4`},{"1966-12",675.5`},{"1967-01",700.1`},{"1967-02",725.1`},{"1967-03",748.2`},{"1967-04",795.4`},{"1967-05",821.8`},{"1967-06",828.9`},{"1967-07",753.1`},{"1967-08",708.9`},{"1967-09",690.9`},{"1967-10",674.5`},{"1967-11",669.6`},{"1967-12",685.3`},{"1968-01",704.`},{"1968-02",730.5`},{"1968-03",760.9`},{"1968-04",807.6`},{"1968-05",842.4`},{"1968-06",838.`},{"1968-07",768.8`},{"1968-08",726.6`},{"1968-09",711.2`},{"1968-10",693.2`},{"1968-11",686.9`},{"1968-12",698.1`},{"1969-01",720.7`},{"1969-02",750.1`},{"1969-03",770.8`},{"1969-04",816.7`},{"1969-05",855.2`},{"1969-06",857.3`},{"1969-07",786.5`},{"1969-08",750.1`},{"1969-09",735.6`},{"1969-10",709.9`},{"1969-11",700.1`},{"1969-12",720.7`},{"1970-01",736.4`},{"1970-02",768.5`},{"1970-03",792.4`},{"1970-04",836.`},{"1970-05",869.9`},{"1970-06",871.5`},{"1970-07",804.1`},{"1970-08",768.8`},{"1970-09",750.8`},{"1970-10",733.4`},{"1970-11",721.4`},{"1970-12",737.4`},{"1971-01",789.4`},{"1971-02",821.8`},{"1971-03",844.4`},{"1971-04",890.8`},{"1971-05",924.9`},{"1971-06",926.3`},{"1971-07",853.2`},{"1971-08",818.9`},{"1971-09",801.5`},{"1971-10",785.5`},{"1971-11",774.1`},{"1971-12",785.5`},{"1972-01",811.`},{"1972-02",838.6`},{"1972-03",873.9`},{"1972-04",913.1`},{"1972-05",943.6`},{"1972-06",948.6`},{"1972-07",877.8`},{"1972-08",839.5`},{"1972-09",820.8`},{"1972-10",795.3`},{"1972-11",777.2`},{"1972-12",790.4`},{"1973-01",806.1`},{"1973-02",840.3`},{"1973-03",867.`},{"1973-04",911.1`},{"1973-05",939.6`},{"1973-06",937.5`},{"1973-07",865.`},{"1973-08",821.8`},{"1973-09",795.4`},{"1973-10",776.6`},{"1973-11",771.1`},{"1973-12",787.4`},{"1974-01",813.`},{"1974-02",845.7`},{"1974-03",872.9`},{"1974-04",915.2`},{"1974-05",951.4`},{"1974-06",960.8`},{"1974-07",891.5`},{"1974-08",851.3`},{"1974-09",826.9`},{"1974-10",797.3`},{"1974-11",784.3`},{"1974-12",798.2`}};

r = Sort@Range[0, 10, 10/(Length[milk]-1)];
data =N@Transpose[{r, milk[[All, 2]]}];

ListPlot[data]

Quit[]

fit1 = FindFormula[data, x]

Show[Plot[fit1, {x, 0, 10}], ListPlot[data, PlotStyle->Gray]]

fit2 = FindFormula[data,x,  TargetFunctions->{Plus, Times, Sin, Cos}]

Show[Plot[{fit1, fit2}, {x, 0, 10}], ListPlot[data, PlotStyle->Gray]]

The milk production is increasing

The milk production is periodic. What is the period?

period =2 Pi/8

Solve[N@period/10 == x/156, x]

Differential Equations

diffEquation = {y''[x]+y[x]==Cos[x],y[0]==y'[0]==0};
diffEquation//Column

sol = NDSolve[{y''[x]+y[x]==Cos[x],y[0]==y'[0]==0},y,{x,0,20}]

r = RandomReal[{0, 20}, 1000];
data = Transpose[{r, sol[[1,1,2]][r]}];
ListPlot[data]

fit = FindFormula[data, x]

Show[Plot[fit, {x, 0, 20}], ListPlot[data, PlotStyle->Gray], ImageSize->500]

What is behind FindFormula

The Problem

In both FindDistribution and FindFormula the task is to find a good model that explains the data, but FindFormula has a hard time doing it!

There are infinitely many models to consider.

Some of them could fit the data so well that they could actually overfit.

Model Selection

FindFormula combines non linear regression with Markov chain Monte Carlo sampling methods.

Non linear regression

Select a finite number of models and estimate parameters.

Many models compete with each other and the best one is selected

({{(Labeled[, "a x"]), (Labeled[, "LogisticSigmoid[ax]"]), (Labeled[, " Sin[α x]"])}, {(Labeled[, "a x^2"]), (Labeled[, "Exp[a x]"]), (Labeled[, "Exp[a x] Sin[β x]"])}})

data2 =Table[{x, Exp[-x]Sin[20x]+RandomVariate[NormalDistribution[0, 0.002]]}, {x, RandomReal[{2, 5}, 2000]}];

fit = FindFormula[data2,x]

Markov Chain Monte Carlo sampling

Quit[]

PacletDirectoryAdd["~/git/ARG/SymbolicMachineLearning"];
Get["SymbolicMachineLearning`"]

The goal is to find the model for which the probability of it given the data is as large as possible. This optimization problem can be tackle using MCMC sampling.

Simulated Annealing: it is used to solve the maximization problem given the fact that the “space” of functions is wild and huge.

ParallelTempering: it is used to sample from the known distribution that provides good fits.

Start with N functions (called replica) e.g {x, x^2, Sin[x], ...}

At each step you propose a move: a move is a choice of N other functions in a neighborhood of them e.g. {2 x, x^2+x, Sin[1.5 x], ...}

For each entry you either accept of reject the corresponding entry of the move with a given probability. This probability depends on the energy function we are trying to maximize/minimize (something related to the mean error for this problem) and on the position of the entry.

{({{(x^2), , x+(x^2)}}),({{x, , 2 x}}),({{Sin[x], , Sin[1.5` x]}}),...}

Once you have determined the new functions, you can choose to swap two of them or not.

({{{(x^2),x,Sin[x],...}, , {(x^2),Sin[x],x,...}}})

SeedRandom[2]
data =Table[{x, x+5 Sin[x]  + 7Cos[x]Sin[x]+RandomVariate[NormalDistribution[0, 0.2]]}, {x, RandomReal[{-9, 9}, 1000]}];
ListPlot[data, PlotStyle->GrayLevel[.50]]

SeedRandom[2]
fit = FindFormula[data,  x, Method->"ParallelTempering", TimeConstraint->9]

Show[ListPlot[data, PlotStyle->GrayLevel[.50]], Plot[fit, {x, -40, 40}, PlotRange->All, PlotStyle->Orange]]

SeedRandom[2]
data =Table[{x,Exp[Sin[x]]+RandomVariate[NormalDistribution[0, 0.002]]}, {x, RandomReal[{2, 40}, 2000]}];
ListPlot[data, PlotStyle->GrayLevel[.50]]

SeedRandom[2]
fit = FindFormula[data,  x, Method->"ParallelTempering", TimeConstraint->9]

Show[ListPlot[data, PlotStyle->GrayLevel[.50]], Plot[fit, {x, 2, 40}, PlotRange->All, PlotStyle->Orange]]

SeedRandom[2]
fit2 = FindFormula[data,  x, Method->"ParallelTempering", TimeConstraint->9, TargetFunctions->{Sin, Exp}]

Show[ListPlot[data, PlotStyle->GrayLevel[.50]], Plot[fit2, {x, 2, 40}, PlotRange->All, PlotStyle->Orange]]

Controlling Overfitting

data =Table[{x,x+RandomVariate[NormalDistribution[0, 1]]}, {x, RandomReal[{-2, 7}, 10]}];
ListPlot[data, PlotStyle->Directive[Gray,PointSize[Large]]]

fit1 = FindFormula[data, x]

fit2 = FindFormula[data, SpecificityGoal->"High", TimeConstraint->9]

Show[ListPlot[data, PlotStyle->Directive[Gray,PointSize[Large]]], Plot[{ fit2[x]}, {x, -2, 8}]]

What will come next?

Generalize FindFormula

Have a higher dimensional version of it.

Make it able to return differential equations solved by data

SeedRandom[2]
data =SortBy[Table[{x,Sin[2 x]+RandomVariate[NormalDistribution[0, 0.002]]}, {x, RandomReal[{2, 40}, 2000]}], First];

findDifferentialRoots[data]

1./4

Graphical Models

Nodes correspond to random variables

Edges represent statistical dependencies between the variables

Graphical models are easy to understand

the graph means P(A, B, C, D, E) = P(A) P(B) P(C|A,B)P(D|B,C)P(E|C,D)

Computing conditional probabilities becomes easier if you know the dependency graph

Make some of the algorithms we use in FindFormula top level functions

Simulated Annealing

Parallel Tempering