Exemplifying the cultural differences between Machine Learning and Statistics

Anton Antonov

Accendo Data LLC

WTC-2016
October 20, Champaign, IL

Join the Conversation #WolframTechConf

Abstract

This talk is going to examine differences between doing Machine Leaning and Statistics. We assume that these two cultures tackle a large set of the same or similar problems hence a comparison is justified.

The differences between the two cultures can be observed by comparing
  1) what is considered important while participating in them, and
  2) the methodologies employed for solving "real life", practical problems.

Obviously, the cultural differences stem from practitioners backgrounds and schooling and we are going to point major differences of educational and experience backgrounds. (Also, jobs.)

General observations will be stated and an overview of concrete examples (in Mathematica) is going to be presented. The concrete examples are from the (sub-)fields of data analysis, prediction, recommendations, time series, gambling mathematics.

Two examples are going to be discussed in more detail:
  1) mashing-up algorithms,
  2) using collaborative filtering instead of regression,
  3) using natural language processing techniques in time series prediction.

As a conclusion we are going to consider how to look for or come up with the new examples showing the cultural differences. (Some of) the major approaches and tricks of each culture are going to be identified and listed.

Background bias

Background

I have strong background in simulations of nature processes. (Industrial problem solving research PhD in large scale air-pollution simulations.)

This means relatively strong background in numerical analysis, functional analysis, and linear algebra. (Also programming and supercomputing.)

I have also been very interested in optimization and operations research. (Delivered combinatorial optimization solutions in different companies in Europe and USA.)

Bias

I have somewhat of disregard for Statistics and works of statisticians.

For the many cases which I was forced to study I did not see the need for the statistical language. It seemed that optimization and linear algebra would suffice.

And I found (find) the inductive jumps made by statisticians annoying.

I do not see this as a good attitude. Merely a valid one.

Role playing

Having a (somewhat) informed hate about Statistics had the effect that in many cases I have/had to play the role of “the statistician“ in companies comprised mostly of software engineers and/or (young) data scientists.

Timid professionals

I tend to think that successful statisticians are more likely to be “timid professionals”, [13], and successful machine learning practitioners tend to be more maverick.

Not Statistics vs. Data science

This talk is not about Statistics vs Data science. (But see some of the references.)

The scope is much narrower : the overlap of Statistics and Machine learning problems.

More precisely we focus here on prediction, classification, recommendations.

Remark: In some sense Data science and Machine learning are part of Statistics -- we deal with data.

Déjà vu all over again

“It’s like déjà vu all over again.”
-- Yogi Berra

Spectral and finite element methods

In many ways this Statistics main MO is like the application of Finite Element Methods (FEM).

Often Statistics we follow this path when dealing with a problem:

1. Get data, make initial observations.

2. Make a model -- say a weighted sum of functions.

3. Fit the model -- extract model parameters by some sort of optimization procedure.

4. Check validity of obtained results

This is also kind of what we do in FEM. The linear model is a weighted sum of vectors in some Hilbert space (i.e. functions). We solve a minimization problem, fitting problem through finding the solution of a large system of linear equations.

In statistics the model fitting happens through Least squares, justified by the Gauss-Markov theorem.

Music

Similar dynamic between new and old practitioners can be observed in music.

Equal temperament, Bach’s “The Well-Tempered Clavier”, correspond to Kolmogorov axiomatic probability theory based on sigma-algebras.

The Pareto law manifestation

Pareto law (Pareto principle, 80/20 rule) is manifest in many human activities and distributions. (Also in nature processes. See this blog post.)

This means that small fraction of the statisticians make the major part of the contributions in Statistics.

Same for Machine learning practitioners.

The best practitioners of both sides have similar effectiveness. We will say that these occupy the top Pareto fraction. The rest are in the corresponding (much larger) bottom Pareto fractions.

Sometimes I might call the bottom Pareto fraction long tail.

We can say that the best Statisticians and ML practitioners are similarly effective using different ideologies, and the long tail Statisticians and ML practitioners have similar attitudes using different tools.

gdps={CountryData[#,"Name"],CountryData[#,"GDP"]}&/@CountryData["Countries"];
gdps=DeleteCases[gdps,{_,_Missing}]/.Quantity[x_,_]:>x;
ParetoLawPlot[gdps[[All,2]],ImageSize->400]
AutoCollapse[]

This sword is not called hate, it is called love

“ This sword is not called hate, it is called love.” 
				-- Paul Claudel, “Jeanne D’Arc Au Bûcher”

There are two ways to deal with the long tail, one is to make it thicker, the other is to cut it.

This talk is more about making the long tail thicker.

I cannot be fair and balanced, but I can be unfair and balanced. That is not something I want to do with this talk. I want it to be on constructive side, ideally this kind of discussions would makes better in practicing Statistics and ML.

The Machine learning protagonists

There will be three ML protagonists in this talk.

Simple Nuclear Physicist (SNP)

Our protagonist is a “Simple Nuclear Physicist” (SNP) -- someone who is accustomed to obtaining a lot of data that has to be analyzed and mined sometimes very deeply, rigorously, and from a lot of angles, for different hypotheses. SNP is fairly adept in programming and critical thinking, but he does not have or care about deep knowledge of statistics methods or machine learning algorithms. SNP is willing and capable to use software libraries that provide algorithms for statistical methods and machine learning.

Common Operations Research Analyst (CORA)

This guy is shared by different departments in some large company like Big Amusement Parks Corporation or Large Vacation Cruises Company. CORA is an internal consultant providing different kinds of decision support for different types of stakeholders. In the Large Vacation Cruises Company CORA is responsible for ships routes design and rerouting, price elasticity models for vacation packages, products and services placement within a ship.

Standard Machine Leaning Engineer (SMLE)

This guy is a typical ML practitioner. Uses ML tools and mindset, might be called data scientist by some.

The major methodological differences to watch for

These are the principles and lesson to watch for in this presentation.

ML lessons

As discussed in Leo Breiman’s paper [1] the major lessons from ML and the algorithmic culture are:

1. Multiplicity of good models (Rashomon).

2. No simple models for good precision results (no Occam).

3. Embracing high dimensionality -- not just a curse. (Bellman).

ML MO

4. Using simulations instead of models.

5. Mash-up of algorithms.

7. Collaborative filtering application.

The examples to examine

We are going to see the ML

The overview list

1. Quantile regression instead of linear regression

2. Using graph based recommendations instead of statistics

3. External data utilization

3a. Time series dependence: formal and informal

4. Computational linguistics methods for time series prediction and analysis

The more detailed exposition list

1. Receiver Operating Characteristic (ROC) for classifier ensembles, damaging, and interpolation

2. Progressive jackpot simulation

3. Extraction of the variable distribution from a noised mixture

Used packages

These commands load the used Mathematica packages used in this talk:

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/MathematicaForPredictionUtilities.m"]
Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/ROCFunctions.m"]
Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/ClassifierEnsembles.m"]
Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/QuantileRegression.m"]
Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/MosaicPlot.m"]

The ROC example

SNP is capable of coming up with Receiver Operating Characteristic (ROC) if he is not aware of it already. ROC is very similar to the so called phase space diagrams physicists do. For a full exposition see  “ROC for classifier ensembles, bootstrapping, damaging, and interpolation”, [7].

Set up

ROC plots evaluation (in brief)

Assume we are given a binary classifier with the class labels P and N (for “positive” and “negative” respectively).

Consider the following measures True Positive Rate (TPR):

TPR:=((correctly classified positives)/(total positives)).

and False Positive Rate (FPR):

FPR:=((incorrectly classified negatives)/(total negatives)).

Assume that we can change the classifier results with a parameter θ and produce a plot like this one:

For each parameter value (Subscript[θ, i]) the point {TPR((Subscript[θ, i])),FPR((Subscript[θ, i]))} is plotted; points corresponding to consecutive (Subscript[θ, i])’s are connected with a line. We call the obtained curve the ROC curve for the classifier in consideration. The ROC curve resides in the ROC space as defined by the functions FPR and TPR corresponding respectively to the x-axis and the y-axis.

The ideal classifier would have its ROC curve comprised of a line connecting {0,0} to {0,1} and a line connecting {0,1} to {1,1}.

Given a classifier the ROC point closest to {0,1}, generally, would be considered to be the best point.

Data used: the Titanic dataset

These commands load the Titanic data (that is shipped with Mathematica).

data=ExampleData[{"MachineLearning","Titanic"},"TrainingData"];
columnNames=(Flatten@*List)@@ExampleData[{"MachineLearning","Titanic"},"VariableDescriptions"];
data=((Flatten@*List)@@@data)[[All,{1,2,3,-1}]];
trainingData=DeleteCases[data,{___,_Missing,___}];
Dimensions[trainingData]

RecordsSummary[trainingData,columnNames]

data=ExampleData[{"MachineLearning","Titanic"},"TestData"];
data=((Flatten@*List)@@@data)[[All,{1,2,3,-1}]];
testData=DeleteCases[data,{___,_Missing,___}];
Dimensions[testData]

RecordsSummary[testData,columnNames]

Replace categorical with numerical values

nTrainingData=trainingData/.{"survived"->1,"died"->0,"1st"->0,"2nd"->1,"3rd"->2,"male"->0,"female"->1};

nTestData=testData/.{"survived"->1,"died"->0,"1st"->1,"2nd"->2,"3rd"->3,"male"->0,"female"->1};

Ensembles

SeedRandom[8989]
aCLs=EnsembleClassifier[Automatic,trainingData[[All,1;;-2]]->trainingData[[All,-1]]];

With Automatic the following built-in classifiers are used:

Keys[aCLs]

If we want to compute ROC of a range of probability thresholds :

rocRange=Range[0,1,0.025];
aROCs=Table[(
cres=EnsembleClassifyByThreshold[aCLs,testData[[All,1;;-2]],"survived"->i,"ProbabilitiesMean"];ToROCAssociation[{"survived","died"},testData[[All,-1]],cres]),{i,rocRange}];
rocEnGr=ROCPlot[rocRange,aROCs,"PlotJoined"->Automatic,PlotLabel->"Classifier ensemble",GridLines->Automatic]

Bootstrapping ensemble with a single classification method

First we select a classification method and make a classifier with it.

clMethod="SupportVectorMachine";
sCL=Classify[trainingData[[All,1;;-2]]->trainingData[[All,-1]],Method->clMethod];

The following code makes a classifier ensemble of 12 classifier functions using resampled, slightly smaller (10%) versions of the original training data (with RandomChoice).

SeedRandom[1262];
aBootStrapCLs=Association@Table[(
inds=RandomChoice[Range[Length[trainingData]],Floor[0.9*Length[trainingData]]];
ToString[i]->Classify[trainingData[[inds,1;;-2]]->trainingData[[inds,-1]],Method->clMethod]),{i,12}];

Let us compare the ROC curves of the single classifier with the bootstrapping derived ensemble.

rocRange=Range[0.1,0.9,0.025];
AbsoluteTiming[
aSingleROCs=Table[(
cres=ClassifyByThreshold[sCL,testData[[All,1;;-2]],"survived"->i];ToROCAssociation[{"survived","died"},testData[[All,-1]],cres]),{i,rocRange}];
aBootStrapROCs=Table[(
cres=EnsembleClassifyByThreshold[aBootStrapCLs,testData[[All,1;;-2]],"survived"->i];ToROCAssociation[{"survived","died"},testData[[All,-1]],cres]),{i,rocRange}];
]

Legended[
Show[{
ROCPlot[rocRange,aSingleROCs,"ROCColor"->Blue,"PlotJoined"->Automatic,GridLines->Automatic],
ROCPlot[rocRange,aBootStrapROCs,"ROCColor"->Red,"PlotJoined"->Automatic]}],
SwatchLegend@@Transpose@{{Blue,Row[{"Single ",clMethod," classifier"}]},{Red,Row[{"Boostrapping ensemble of\n",Length[aBootStrapCLs]," ", clMethod," classifiers"}]}}]

We can see that we get much better results with the bootstrapped ensemble.

Damaging: decreasing proportions of females

First, let us find all indices corresponding to records about females.

femaleInds=Flatten@Position[trainingData[[All,3]],"female"];

The following code standardizes the training data corresponding to females, finds the mean record, computes distances from the mean record, and finally orders the female records indices according to their distances from the mean record.

t=Transpose@Map[Rescale@*Standardize,N@Transpose@nTrainingData[[femaleInds,1;;2]]];
m=Mean[t];
ds=Map[EuclideanDistance[#,m]&,t];
femaleInds=femaleInds[[Reverse@Ordering[ds]]];

The following plot shows the distances calculated above.

ListPlot[Sort@ds,PlotRange->All,PlotTheme->"Detailed"]

The following code removes from the training data the records corresponding to females according to the order computed above. The female records farthest from the mean female record are removed first.

AbsoluteTiming[
femaleFrRes=Association@
Table[cl->
Table[(
inds=Complement[Range[Length[trainingData]],Take[femaleInds,Ceiling[fr*Length[femaleInds]]]];
cf=Classify[trainingData[[inds,1;;-2]]->trainingData[[inds,-1]],Method->cl];cfPredictedLabels=cf/@testData[[All,1;;-2]];
{fr,ToROCAssociation[{"survived","died"},testData[[All,-1]],cfPredictedLabels]}),
{fr,0,0.8,0.05}],
{cl,{"NearestNeighbors","NeuralNetwork","LogisticRegression","RandomForest","SupportVectorMachine","NaiveBayes"}}];
]

The following graphics grid shows how the classification results are affected by the removing fractions of the female records from the training data. The results for none or small fractions of records removed are more blue.

GraphicsGrid[ArrayReshape[
Table[
femaleAROCs=femaleFrRes[cl][[All,2]];
frRange=femaleFrRes[cl][[All,1]];ROCPlot[frRange,femaleAROCs,PlotRange->{{0.0,0.25},{0.2,0.8}},PlotLabel->cl,"ROCPointColorFunction"->(Blend[{Blue,Red},#3/Length[frRange]]&),ImageSize->300],
{cl,Keys[femaleFrRes]}],
{2,3}],Dividers->All]

We can see that removing the female records outliers has dramatic effect on the results by the classifiers “NearestNeighbors” and “NeuralNetwork”. Not so much on “LogisticRegression” and “NaiveBayes”.

Simple math

See the article “An introduction to ROC analysis” by Fawcett. The math is very simple (pre-high school level, ratios). I.e. not statistics hypotheses, i.i.d. etc. is required to do ROC. Just common sense.

The ROC example: Linear regression comparison

Do linear regression

lfmFunc=LinearModelFit[{nTrainingData[[All,1;;-2]],nTrainingData[[All,-1]]}]

It is a good idea to evaluate the regression results. One way is to verify the assumption of i.i.d of errors in the equation / model used for Linear regression:

(Subscript[y, i])==α+(Subscript[β, i])(Subscript[x, i]) + (Subscript[ϵ, i]), i∈[1,n],

where (Subscript[ϵ, i]) are assumed to be i.i.d. (See also the Gauss-Markov theorem.)

In matrix-vector notation:

y==X β + ϵ, y∈(ℝ^(n×1)), X∈(ℝ^(n×p)),β∈(ℝ^(p×1)), ϵ∈(ℝ^(n×1)).

Evaluation of the linear regression

(*Predicted vs.residual*)
ListPlot[Transpose[{lfmFunc["PredictedResponse"],lfmFunc["FitResiduals"]}],Frame->True,FrameLabel->{{"Residual",""},{"Predicted","Predicted vs. Residual"}}]

(*Quantile plot*)
QuantilePlot[lfmFunc["FitResiduals"]]

(*Histogram of residuals*)
Show[Histogram[lfmFunc["FitResiduals"],"FreedmanDiaconis","PDF",Frame->True,FrameLabel->{{"Probability density",""},{"Residual","Histogram"}}],Plot[PDF[NormalDistribution[0,StandardDeviation[lfmFunc["FitResiduals"]]]][x],{x,-150,150}]]

We can use Cook’s distances to evaluate the influence of outliers of the data to the

ListPlot[lfmFunc["CookDistances"],PlotRange->All,PlotTheme->"Detailed"]

One prescription is to watch for and explain anything above this threshold:

4./Length[trainingData]

Progressive jackpot modeling

In this case CORA is working on a Boat Casino (somewhere few hours away from here.)

The exposition of this too long to include in this presentation notebook. See this Markdown file “Progressive jackpot modeling.md”.

Notebook link, “Progressive jackpot modeling.nb”.

Problem set-up

Neat math modeling vs simulations

“We are going to abandon this “neat math” approach in favor of a different one based on stochastic processes simulation. The stochastic simulation approach is more useful if the model has to be expanded in order to incorporate different features and sub-processes.”

Distribution extraction from noised mixture

This is a task SMLE is had to work on for some governmental department.

Problem formulation

This document started as an answer for the Mathematica StackExchange discussion:

http://mathematica.stackexchange.com/questions/108182/extracting-signal-from-gaussian-noise

We want to find the distribution of the variable X that used to generate the data Z := X +Y, where is Y has the distribution N(0,σ). The standard deviation σ is known. We are interested in finding the CDF and/or PDF of X.

Here is a surrogate distribution for the type of X that is expected. This is defined:

SeedRandom[1256]
fSurvivalGompertzDistRand[α_,β_]:=ProbabilityDistribution[(1/((E^(α/β) Gamma[0,α/β])/β) E^(((1-E^(t β)) α)/β)),{t,0,∞}]

This generates the part for X:

dataX=RandomVariate[fSurvivalGompertzDistRand[0.016,0.65],{20000}];
Histogram[dataX,50,PlotTheme->"Detailed"]

Then we add some Gaussian noise:

Manipulate[(
dataNoise=dataX+RandomVariate[NormalDistribution[0,σ],{20000}];Histogram[dataNoise,60,PlotLabel->Row[{"σ=",σ}],PlotTheme->"Detailed",PlotRange->{0,1500}]),
{{σ,1},0.1,2.5,0.1}]

snapshot

DynamicModule[{σ=0.2`},dataNoise=dataX+RandomVariate[NormalDistribution[0,σ],{20000}];Histogram[dataNoise,60,PlotLabel->Row[{"σ=",σ}],PlotTheme->"Detailed",PlotRange->{0,1500}]]

General procedure

The procedure is as follows:

1. Find a function F that applied to a collection of real numbers produces one value (say, 0) for normally distributed data and other different values for non-normally distributed data. (F can be Kurtosis.)

2. Pick a sample (Subscript[{(Subscript[s, i])}, i=1]^n) of the noised data with a relatively small number of n points (say n ∈ [20,40]).

3. Formulate an optimization problem with n variables (Subscript[{(Subscript[v, i])}, i=1]^n) that maximizes |F(s-v)| subject to constrains that would hold if (Subscript[{(Subscript[v, i])}, i=1]^n) come from a normal distribution (with known parameters).

4. Solve the optimization problem several times with different samples and accumulate the sets (Subscript[{(Subscript[s, i])-(Subscript[v, i])}, i=1]^n)  and (Subscript[{(Subscript[v, i])}, i=1]^n). Monitor the (χ^2) test over (Subscript[{(Subscript[v, i])}, i=1]^n) .

5. Reconstruct the original distribution (“signal”) CDF and PDF using quantiles of union of (Subscript[{(Subscript[s, i])-(Subscript[v, i])}, i=1]^n) from all optimization runs.

Full exposition of the solution

The full solution with graphs an code can be read here : “Distribution extraction from a Gaussian noised mixture”, [8].

Notebook link “Distribution-extraction-from-a-Gaussian-noised-mixture”.

Overview example 1: Quantile regression

See this MathematicaForPrediciton at WordPress category list for discussions of implementing and using Quantile Regression with Mathematica; or [11]. Most relevant for this discussion is probably “Estimation of conditional density distributions”, [14].

location={"Atlanta","GA"};
(*location="Melbourne";*)
tempData=WeatherData[location, "Temperature",{{2010, 1, 1}, {2016, 3, 12}, "Day"}];
tempData//Length

tempPData=Partition[Normal@tempData["Values"],2,1];
tempPData=tempPData/.Quantity[q_,_]:>q;
tempPData=Select[tempPData,VectorQ[#,NumberQ]&];
tempPData//Dimensions

grDLP=ListLinePlot[tempData,PlotRange->All,AspectRatio->1/3,PlotTheme->"Detailed"]

tempDataArray=Transpose[{tempData["TimeList"][[1]],Normal[tempData["Values"]]/.Quantity[x_,_]:>x}];

Find regression quantiles:

qs=Join[{0.02},FindDivisions[{0,1},6][[2;;-2]],{0.98}]//N
qs//Length

AbsoluteTiming[
qFuncsDLP=QuantileRegression[tempDataArray,26,qs,Method->{LinearProgramming,Method->"CLP"}];
]

Plot regression quantiles and time series data:

AbsoluteTiming[grQFuncsDLP=Plot[Evaluate[Through[qFuncsDLP[x]]],{x,Min[tempDataArray[[All,1]]],Max[tempDataArray[[All,1]]]},PerformanceGoal->"Speed",PlotPoints->130];
]

Show[{grDLP,grQFuncsDLP}]

Overview example 2: NLP for time series collections

Instead of doing time series analysis we can try to convert the time series problems into Natural Language Processing (NLP) problems.

More precisely, using suitable discretization we convert time series problems into sequence mining problems, in some cases with large enough collections we can go further at apply NLP algorithms.

Regression quantiles generated sequence of states

Quantile regression fitting

We have time series (Atlanta, GA, temperature) with fitted regression quantiles.

States

This procedure is somewhat reminiscent to making dimensionless models in Physics and Chemistry.

The regression quantiles generate bands. (Seven in this case.) We replace each time series point with the band number it belongs to.

Shallow[qstates]

Now we can do sequence mining.

Here is the mosaic plot picturing the conditional probabilities of the state occurrence:

MosaicPlot[Partition[qstates,2,1]]

Overview example 3: Outside data, Elance vs. oDesk

Set-up

Elance connects employers from around the world with free lancers from around world.

Each Elance customer who wants to employ a free-lancer (or five) posts a job description and receives bids for it from freelancers.

Elance takes a percentage of financial transactions.

Several questions are of interest:

1. How employers select from the bids? I.e. given a job and, say, 20 bids which 2-3 bids would be selected to final evaluation and awarding by the job poster?

2. How the previous performance of a freelancer -- money amount earned, star rating, number of successful completions, test scores -- would influence freelancer’s selection for a job?

3. How to devise a recommendation system for successful bids and evaluate it?

4. How to recommend freelancers for job?

5. How to deal with plagiarism in bids or jobs descriptions?

The solution

The original idea was to use regression and classification. I had to normalize the data in a couple of different ways in order to get reasonable performance. Which was still not considered satisfactory. I tried to explain that given the data this is the best we can do, which was not accepted by stakeholders much since it seemed contrived and the stakeholders were very confident that the necessary signals were in the data Elance has.

What worked was to show that that last assumption is not true. Since I knew from the regression / classification efforts that in many cases freelancers with zero performance (no money earned, no star rating, etc.) are selected for jobs, I looked into the competitor’s site, oDesk. Found that many of Elance-zero-performance freelancers have stellar reputation in oDesk. The employers simply preferred Elance for hiring, but were doing the search and selection of candidates at oDesk.

After explaining this to stakeholders I somewhat jokingly mentioned that oDesk and Elance should share their data, or join forces in other way. Six-seven months after I finished my contract with Elance, Elance and oDesk announced that they are merging into one company.

How it was done?

I had a relatively simple Mathematica script that was gathering data from Elance and oDesk for a given occupation tier (e.g. Accounting, Software Engineering, etc.) and calculates a matching score for a pair of profiles (one from each site). I did this for around 600 profiles across several different tiers -- I just had to prove a point.

Turned out that a more complete and systematic version of this kind competitor site mining is done a lot by dating sites. For them it is important to know which / how many people have “found someone” or simply switched to another site.

Why is this mentioned here?

Well ultimately it was not Statistical methods that provided the answer -- just data science and machine learning was used.

Overview example 4: Using graphs based recommendations instead of regression

I am not comfortable showing the actual data set-up of this project the results of the solution that worked well. I am going to describe the set-up and the algorithm.

For an example of the bi-partite graphs this recommender works on see the neat example in this blog post: “RSparseMatrix for sparse matrices with named rows and columns”.

Import["https://mathematicaforprediction.files.wordpress.com/2015/10/movies-actors-graph.png"]

Morphological analysis of behavior

Statistics

Frequentist/objective statistics tells us what to do if we have a model. Bayesian/subjective statistics tells us what to do if we know (or have premonition) of the probabilities priors.

1. Gather, transform, clean, normalize data.

2. Look how to apply a model with i.i.d assumption. 
    (Come up with a model.)

3. Extract parameters by some fitting or optimization procedure.
    (Model fitting.)

4. Verify model assumptions from obtained results.
     I.i.d., white noise, quantile plots.

5. Test and verify prediction.
     Goodness of fit tests.

Machine learning

See “Morphological analysis of algorithms vs problems” and “Morphological analysis of algorithms vs data” .

Building up an example from scratch

Statistics

1. Get time series retail data from different geographic locations.

2. Aggregate time series spatially so the law of large numbers takes effect.

3. Average the time series into one for each cluster (with some normalization applied).

4. Apply time series methods (ARMA, ARIMA, etc.)
    Certain monthly or weekly seasonality would be observed.

5. Verify prediction results.

Machine learning

1. Get time series retail data from different geographic locations.

2. Cut the time series into weekly wave fragments.

3. Cluster the weekly wave fragments or apply some other unsupervised learning process.
    E.g. Spearman correlation combined with NLP application.

4. Make predictions based on cluster membership of Nearest Neighbors (NNs).

5. Verify prediction results.
    E.g. A/B testing.

Summary

These two MNIST recognition bases can be used to summarize differences between Statistics and Machine Leaning.

SVD/PCA -- Statistics

Quick algorithm, theoretically well justified. We have a mean image that is modified with difference images into an unknown image.

NNMF -- Machine Learning

Using an optimization algorithm with lots of local minima. Easy to interpret results. An unknown digit is approximated as most likely combination of several known ways to write “5”.

References

[1] George E.P. Box, Science and Statistics, Journal of the American Statistical Association, Vol. 71, No. 356. (Dec., 1976), pp. 791-799. URL: http://mkweb.bcgsc.ca/pointsofsignificance/img/Boxonmaths.pdf

[2] Leo Breiman, Statistical Modeling: The Two Cultures with Comments, Statistical Science, Statistical Science 2001, Vol. 16, No. 3, 199–231.  URL: https://projecteuclid.org/euclid.ss/1009213726

[3] David Donoho, 50 years of Data Science, 2015, URL: http://courses.csail.mit.edu/18.337/2015/docs/50YearsDataScience.pdf

[4] John W. Tukey, The Future of Data Analysis, The Annals of Mathematical Statistics Vol. 33, No. 1 (Mar., 1962), pp. 1-67, URL: https://www.jstor.org/stable/2237638

[5] Brendan O'Connor, "Statistics vs. Machine Learning, fight!", (2008), URL: http://brenocon.com/blog/2008/12/statistics-vs-machine-learning-fight/

[6] "The Two Cultures: statistics vs. machine learning?", Cross Validated, URL: http://stats.stackexchange.com/questions/6/the-two-cultures-statistics-vs-machine-learning

[7] Anton Antonov, “ROC for classifier ensembles, bootstrapping, damaging, and interpolation”, (2016), MathematicaForPrediction at GitHub, folder Documents. WordPres blog post.

[8] Anton Antonov, “Distribution extraction from a Gaussian noised mixture”, (2016), MathematicaVsR at GitHub, project DistributionExtractionAFromGaussianNoisedMixture .

[9] Tom Fawcett, An introduction to ROC analysis, (2006), Pattern Recognition Letters, 27, 861–874. (Link to PDF.)

[10] Anton Antonov, “RSparseMatrix for sparse matrices with named rows and columns”, (2015), MathematicaForPrediction at WordPress.

[11] Anton Antonov, “Quantile regression through linear programming”, (2013), MathematicaForPrediction at WordPress.

[12] Anton Antonov, “Morphological analysis of data mining and prediction algorithms”, (2013), MathematicaForPrediction at WordPress.

[13] Jeff Schmidt, ‪Disciplined Minds: A Critical Look at Salaried Professionals and the Soul-battering System that Shapes Their Lives‬, ‪Rowman & Littlefield, 2001‬.

[14] Anton Antonov, “Estimation of conditional density distributions”, (2014), MathematicaForPrediction at WordPress.