# Problem fitting density function data
**URL:** https://discourse.julialang.org/t/problem-fitting-density-function-data/97085
**Category:** Statistics
**Created:** [April 4, 2023, 8:53pm UTC](https://discourse.julialang.org/t/problem-fitting-density-function-data/97085 "2023-04-04T20:53:38Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@daro.slai](https://discourse.julialang.org/u/daro.slai)
#### Post date: [April 4, 2023, 8:53pm UTC](https://discourse.julialang.org/t/problem-fitting-density-function-data/97085/1 "2023-04-04T20:53:38Z")
I have a dataset from [Elaad](https://platform.elaad.io/analyses/index.php?url=ElaadNL_opendata.php).
Which looks like this:
the idea is to fit a given `Distribution::D` to the data. I wanted to try a couple of ideas but the fit needs to handle both the bins and the probability data.
I´ve tried using ´StatsBase.fit()´and ´Distributions.fit()´ but either they don´t handle ´weights´ or don´t have the ´suffstats´ of any of the asymmetric distributions I´ve tried.
My min work example for now is:
```julia
using StatsBase, StatsPlots, Distributions, Plots, CSV, DataFrames, LaTeXStrings
arr_times_df = CSV.read("distribution-of-arrival-weekdays.csv", DataFrame);
# Normalize the arrival times to 1
for i in 2:size(arr_times_df, 2)
arr_times_df[:, i] = arr_times_df[:, i]/sum(arr_times_df[:, i])
end
# add extra column with arrival times in secs Float64
arr_times_df[!, :arrival_times] = (0:0.25:23.75)*3600.0;
# Convert data to arrays
times = convert(Vector, arr_times_df.arrival_times)
prob = convert(Vector, arr_times_df.private)
# Fit a probability distribution to the data
dist = fit(TriangularDist, times, weights(prob))
```
which of course doesn´t work.
---
### Author:  [@tchebycheff](https://discourse.julialang.org/u/tchebycheff)
#### Post date: [April 5, 2023, 2:58pm UTC](https://discourse.julialang.org/t/problem-fitting-density-function-data/97085/2 "2023-04-05T14:58:04Z")
You already have the density from the data, although in a histogram format. If you are looking for a smooth curve, perhaps you can try [GaussianMixtures](https://github.com/davidavdav/GaussianMixtures.jl) package.
Looking at the shape, you could fit a mixture of 2 or 3 Gaussians. Construct a GMM instance and fit to the observed likelihood?
---
### Author:  [@daro.slai](https://discourse.julialang.org/u/daro.slai)
#### Post date: [April 13, 2023, 7:58pm UTC](https://discourse.julialang.org/t/problem-fitting-density-function-data/97085/3 "2023-04-13T19:58:30Z")
Well it actually took me a lot of time to figure out `GaussianMixtures.jl` and the rest of the `StatsBase.jl` environment, but this is what I came out with.
```julia
using StatsBase, StatsPlots, LinearAlgebra, IterTools, Distributions, Plots, CSV, DataFrames,
using GaussianMixtures
function make_mixture(h)
MixtureModel([Uniform(e...) for e in partition(h.edges[1], 2, 1)],
normalize(h.weights, 1))
end
edges=collect(0.0:0.25:24)
weights=arr_times_df.private;
h=Histogram(edges, weights)
m = make_mixture(h) # create the mixture empirical models of bins
gmm = MixtureModel(GMM(2, rand(m, 5000))); # create the GMM model with 2 peaks and 5000 samples out of the mixture model
print(gmm)
K-means converged with 12 iterations (objv = 15000.248816129308)
┌ Info: Initializing GMM, 2 Gaussians diag covariance 1 dimensions using 5000 data points
└ @ GaussianMixtures C:\Users\daros\.julia\packages\GaussianMixtures\zDaBV\src\train.jl:79
┌ Info: K-means with 2000 data points using 12 iterations
│ 500.0 data points per parameter
└ @ GaussianMixtures C:\Users\daros\.julia\packages\GaussianMixtures\zDaBV\src\train.jl:141
MixtureModel{Normal{Float64}}(K = 2)
components[1] (prior = 0.8222): Normal{Float64}(μ=19.447283283332077, σ=2.401753639052661)
components[2] (prior = 0.1778): Normal{Float64}(μ=10.907372546632375, σ=5.855637979445984)
```
Sadly I think the fit is not very good
```julia
plot(h)
plot!(gmm, func=pdf,linestyles=[:dash], linecolor=:lightblue, label="Components")
density!(rand(gmm, 100_000), linecolor=:lightblue, linewidth=2, xlabel = "Value", label = "Mixture model")
```
if i normalize `h` with`normalize(h, mode=:density)` then the fit looks better:
but then `sum(h.weights)=4.0` so…
---
### Author:  [@tbeason](https://discourse.julialang.org/u/tbeason)
#### Post date: [April 13, 2023, 8:11pm UTC](https://discourse.julialang.org/t/problem-fitting-density-function-data/97085/4 "2023-04-13T20:11:02Z")
Since you are working with a time, and time cannot be negative, you should incorporate this information into the estimation somehow. A convenient way, if time is always nonzero, is to work with log(time).
---
### Author:  [@daro.slai](https://discourse.julialang.org/u/daro.slai)
#### Post date: [July 20, 2023, 1:43pm UTC](https://discourse.julialang.org/t/problem-fitting-density-function-data/97085/5 "2023-07-20T13:43:33Z")
Well my application is not super strict so I can just sample from the model, check the value and if its outside of the limits [0:24] I re-sample.
Now I’m trying to make forecasts out of this data. I’m more used to forecasting timeseries/statespace models, so any recommendation is more than welcomed.