# MLJ probabilistic models basics
**URL:** https://discourse.julialang.org/t/mlj-probabilistic-models-basics/36009
**Category:** General Usage
**Tags:** question, mlj
**Created:** [March 15, 2020, 1:15pm UTC](https://discourse.julialang.org/t/mlj-probabilistic-models-basics/36009 "2020-03-15T13:15:33Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@jw3126](https://discourse.julialang.org/u/jw3126)
#### Post date: [March 15, 2020, 1:15pm UTC](https://discourse.julialang.org/t/mlj-probabilistic-models-basics/36009/1 "2020-03-15T13:15:33Z")
I am new to MLJ and trying to do a linear regression:
```julia
using MLJ
nexamples = 10
nfeatures = 3
X = randn(nexamples, nfeatures)
A = randn(nfeatures, 1)
b = randn()
noise = 0.1*randn(nexamples)
y = vec(X * A .+ b) + noise
X = MLJ.table(X)
model = @load LinearRegressor pkg = GLM
mach = machine(model, X, y)
fit!(mach)
y_hat = predict(mach)
```
This gives me a vector of probability distributions, which is really cool.
```julia
5-element Array{Distributions.Normal{Float64},1}:
Distributions.Normal{Float64}(μ=-5.217306370475575, σ=0.04688655865968578)
Distributions.Normal{Float64}(μ=-1.6777270495732746, σ=0.04688655865968578)
Distributions.Normal{Float64}(μ=-1.4193911630828542, σ=0.04688655865968578)
Distributions.Normal{Float64}(μ=-0.05493986909253701, σ=0.04688655865968578)
Distributions.Normal{Float64}(μ=-1.7665791465850154, σ=0.04688655865968578)
```
However I would like to `evaluate` the model and also compare it against deterministic models.
How to do that? One idea would be to make the model deterministic, e.g. by taking the mean. How to code that?
---
### Author:  [@tlienart](https://discourse.julialang.org/u/tlienart)
#### Post date: [March 15, 2020, 5:22pm UTC](https://discourse.julialang.org/t/mlj-probabilistic-models-basics/36009/2 "2020-03-15T17:22:32Z")
you can use `predict_mean` (also in pipelines)
---
### Author:  [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)
#### Post date: [March 15, 2020, 11:24pm UTC](https://discourse.julialang.org/t/mlj-probabilistic-models-basics/36009/3 "2020-03-15T23:24:39Z")
@jw3126 @tlienart
This looks really cool.
In the example above only the conditional expectation `E[Y|X]` changes w/ X.
This is assuming
1. homoskedasticity `Var[Y|X]` is constant.
2. the conditional distribution `Y|X` is normal (due to GLM)
In reality that’s usually not likely.
I tried this w/ the Boston housing data & got a constant `skedastic` function
```julia
X, y = @load_boston;
model = @load LinearRegressor pkg = GLM
mach = machine(model, X, y)
fit!(mach)
y_hat = predict(mach)
Distributions.Normal{Float64}(μ=30.212372064783466, σ=4.787101274343532)
Distributions.Normal{Float64}(μ=25.267233800068816, σ=4.787101274343532)
Distributions.Normal{Float64}(μ=30.849358585028362, σ=4.787101274343532)
```
Assuming you have enough data (you’re gonna need a lot), is it possible to use MLJ to get heteroskedastic predictions?
This is particularly useful in finance/insurance where users care alot more about `σ` then `μ`.
Is it possible to get predictions which are not always normal? Either best parametric fit (from some set of models, or nonparametric)
(I realize we’re prob gonna have to depart from linear models…)
For example suppose X is 1-dimensional:
```julia
For X ∈ (1.0,1.05), Y|X ∼ Distributions.Normal{Float64}(μ=3.3, σ=4.8)
For X ∈ (2.2,2.35), Y|X ∼ Distributions.Normal{Float64}(μ=2.3, σ=1.6)
For X ∈ (2.8,3.10), Y|X ∼ Distributions.LogNormal{Float64}(μ=5.0, σ=3.6)
For X ∈ (3.5,3.72), Y|X ∼ NonParam{Float64}(μ̂ =4.0, σ̂ =3.3)
```
In grad school it is very routine to estimate `Var[b|X]` assuming heteroskedasticity using robust EHW standard errors, or clustered se, (or FGLS w/ ML methods), then we can infer `E[Y|X]`. But only for linear models.
Is there a literature on this I’m not aware of? @Tamas_Papp @fipelle
I’d ultimately like to produce something like @mthelm85’s elegant figure:
**update 1** : @fipelle suggests Adrian et al 2019. ([code in Matlab](https://www.openicpsr.org/openicpsr/project/113169/version/V1/view?flag=follow&pageSize=100&sortOrder=(?title)&sortAsc=true))
Figure 1: One-year-ahead predictive distribution of real GDP growth, based on quantile regressions with current real GDP growth and NFCI as conditioning variables
A shortcoming of this method is that it only produces unimodal predictions…
Technically this is doable in MLJ which has options for quantile regression.
In addition, the authors fit a skewed t -distribution (Azzalini & Capitanio 2003) to smooth the quantile function and recover a probability density function.
Azzalini has an [R package](https://cran.r-project.org/web/packages/sn/index.html) for the skew-normal/t-distribution.
The skewed-t is currently not in [Distributions.jl](https://github.com/JuliaStats/Distributions.jl) (@simonbyrne @johnmyleswhite @andreasnoack)
**Update 2** : it looks like @oxinabox’s [DensityEstimationML.jl](https://github.com/oxinabox/DensityEstimationML.jl) can be helpful as well but I don’t see any concise examples right now.
**Update 3** : an emerging literature extends boosting to probabilistic forecasting.
[XGBoostLSS](https://github.com/StatMixedML/XGBoostLSS), [CatBoostLSS](https://github.com/StatMixedML/CatBoostLSS), [Gamlss](https://cran.r-project.org/web/packages/gamlss/index.html), [GamboostLSS](https://cran.r-project.org/web/packages/gamboostLSS/index.html), [bamlss](https://cran.r-project.org/web/packages/bamlss/index.html), [disttree](https://rdrr.io/rforge/disttree/), [ngboost](https://github.com/stanfordmlgroup/ngboost/).
Beautiful: [Slides](https://drive.google.com/file/d/183BWFAdFms81MKy6hSku8qI97OwS_JH_/view) & [Docs](https://github.com/stanfordmlgroup/ngboost) & [Guide](https://stanfordmlgroup.github.io/ngboost/intro.html)
[XGBoostLSS Paper](https://arxiv.org/pdf/1907.03178.pdf):
“The ultimate goal of regression analysis is to obtain information about the **[entire] conditional distribution** of a response given a set of explanatory variables.” (Hothorn et al., 2014, emphasis added)
Julia can really shine here!
---
### Author:  [@fipelle](https://discourse.julialang.org/u/fipelle)
#### Post date: [March 16, 2020, 1:33am UTC](https://discourse.julialang.org/t/mlj-probabilistic-models-basics/36009/4 "2020-03-16T01:33:40Z")
Recently, I have seen a lot of empirical work (in macroeconomics and macro-finance) on quantile regressions (for example: [Adrian, Boyarchenko and Giannone, 2019](https://www.aeaweb.org/articles?id=10.1257/aer.20161923)). However, I am not exactly sure from your post if you are more interested in estimating the conditional mean or variance. If you can clarify further on that I can try to be more helpful.
---
### Author:  [@Albert\_Zevelev](https://discourse.julialang.org/u/Albert_Zevelev)
#### Post date: [March 16, 2020, 1:52am UTC](https://discourse.julialang.org/t/mlj-probabilistic-models-basics/36009/5 "2020-03-16T01:52:46Z")
I’d like to estimate the entire conditional distribution when possible.
If not, @ least a heteroskedastic normal.
Thanks for link.
Looks like they use a skewed t-distribution developed by Azzalini and Capitanio (2003), coded in Matlab.
[https://www.openicpsr.org/openicpsr/project/113169/version/V1/view?path=/openicpsr/113169/fcr:versions/V1/azzalini&type=folder](https://www.openicpsr.org/openicpsr/project/113169/version/V1/view?path=/openicpsr/113169/fcr:versions/V1/azzalini&type=folder)
I’ve never seen this.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [March 16, 2020, 8:22am UTC](https://discourse.julialang.org/t/mlj-probabilistic-models-basics/36009/6 "2020-03-16T08:22:26Z")
> [@Albert\_Zevelev](#):
>
> I’d like to estimate the entire conditional distribution when possible.
I would go with a parametrized Bayesian model (then model checking, and extending as necessary, repeated until satisfied), but that has nothing to do with MLJ so that would be off-topic here. Also, it is something that you can easily do with existing MCMC tools in Julia.
---
### Author:  [@jw3126](https://discourse.julialang.org/u/jw3126)
#### Post date: [March 16, 2020, 10:03am UTC](https://discourse.julialang.org/t/mlj-probabilistic-models-basics/36009/7 "2020-03-16T10:03:09Z")
Thanks! I tried the following, but it just makes julia hang without output:
```julia
using Pkg
pkg"""
up
add MLJ MLJModels GLM MLJBase MLJModelInterface
"""
using MLJ
nfeatures = 3
nexamples = 10
X = randn(nexamples, nfeatures)
A = randn(nfeatures, 1)
noise = zeros(nexamples)
noise = randn(nexamples)
y = vec(X * A) + noise
X = MLJ.table(X)
Xs = source(X)
ys = source(y, kind=:target)
pmodel = @load LinearRegressor pkg=GLM
pmach = machine(pmodel, Xs, ys)
y_hat = predict_mean(pmach)
model = @from_network Det(pmodel=pmodel) <= y_hat
evaluate(model, X, y, measure=rms)
```
---
### Author:  [@tlienart](https://discourse.julialang.org/u/tlienart)
#### Post date: [March 16, 2020, 7:08pm UTC](https://discourse.julialang.org/t/mlj-probabilistic-models-basics/36009/8 "2020-03-16T19:08:47Z")
You forgot a fit in the mix, how about:
```julia
using MLJ
n, p = 10, 3
X = randn(n, p)
A = randn(p)
y = X*A + randn(n)
Xt = MLJ.table(X)
@load LinearRegressor pkg=GLM
pmach = machine(LinearRegressor(), Xt, y)
fit!(pmach)
y_hat = predict_mean(pmach)
```
---
### Author:  [@jw3126](https://discourse.julialang.org/u/jw3126)
#### Post date: [March 16, 2020, 7:56pm UTC](https://discourse.julialang.org/t/mlj-probabilistic-models-basics/36009/9 "2020-03-16T19:56:06Z")
Thanks, that does work. However it is not quite what I am looking for. Really I would like to build a first class deterministic model out of the probabilistic one. For instance I would like to `evaluate` it and compare against other deterministic models using `rms` measure.
---
### Author:  [@oxinabox](https://discourse.julialang.org/u/oxinabox)
#### Post date: [March 27, 2020, 9:26am UTC](https://discourse.julialang.org/t/mlj-probabilistic-models-basics/36009/10 "2020-03-27T09:26:29Z")
> [@Albert\_Zevelev](#):
>
> Update 2: it looks like @oxinabox’s [DensityEstimationML.jl](https://github.com/oxinabox/DensityEstimationML.jl) can be helpful as well but I don’t see any concise examples right now.
That was a research project that was never completed,
I think the math in there all holds and may actually be at least partially novel.
(though on large bit of it I found out after writing 80% of the paper on it, definately is not novel)
I have now archieved the repo.
Idk how much of the stuff that ended up in it actually worked or was usable.