# Learning statistics: seeking advice on how to model a 2d-size distribution of clothing items

**URL:** <https://discourse.julialang.org/t/learning-statistics-seeking-advice-on-how-to-model-a-2d-size-distribution-of-clothing-items/71152>\
**Category:** Statistics\
**Created:** [November 8, 2021, 2:46pm UTC](https://discourse.julialang.org/t/learning-statistics-seeking-advice-on-how-to-model-a-2d-size-distribution-of-clothing-items/71152 "2021-11-08T14:46:00Z")\
**Posts on this page:** 6\
**Page:** 2

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [November 10, 2021, 4:21pm UTC](https://discourse.julialang.org/t/learning-statistics-seeking-advice-on-how-to-model-a-2d-size-distribution-of-clothing-items/71152/21 "2021-11-10T16:21:04Z")

</div>

> [@tim\_pi](#):
>
> I’m not really looking to forecast sales, rather the distribution of sizes of said sales.

Well, instead of trying to forecast the sum of a bunch of things (Which is relatively more easy) you are trying to forecast **each** of the individual things (which is harder!).

---

<div class="post-metadata">

**Author:** ![mthelm85](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mthelm85/32/224164_2.png) [@mthelm85](https://discourse.julialang.org/u/mthelm85)\
**Post date:** [November 10, 2021, 6:04pm UTC](https://discourse.julialang.org/t/learning-statistics-seeking-advice-on-how-to-model-a-2d-size-distribution-of-clothing-items/71152/22 "2021-11-10T18:04:22Z")

</div>

> [@tim\_pi](#):
>
> Yes! I think I’ll try that. Thank you for your suggestion.

Here’s a quick example using Turing and your test scores idea. Let’s say a teacher gave a particular exam to 200 students total over the course of last year. This year, the teacher wants to estimate what the year-end distribution of test scores is likely to be, given last year’s data and the test scores for 10 students that have taken the exam so far this year.

```julia
using Distributions
using MCMCChains
using StatsPlots
using Turing

d_old = Normal(70.0, 25.0)
d_new = Normal(78.0, 20.0) # this year's class is smarter

scores_old = rand(d_old, 200)
prior = fit(Normal, scores_old)

@model function test_scores(scores)
    μ ~ Normal(prior.μ, 25)
    σ ~ Uniform(0, 50)
    scores ~ Normal(μ, σ)
end

N = 10
scores_new = rand(d_new, N)

chain = sample(test_scores(scores_new), NUTS(), 1000)

μ_summary = chain[:μ]
σ_summary = chain[:σ]

plot(
    prior,
    label="Prior",
    xlabel="Score",
    ylabel="Density",
    title="N = $N \n CI: ($(Int64(round(quantile(vec(μ_summary), 0.025)))), $(Int64(round(quantile(vec(μ_summary), 0.975)))))"
)
vline!([prior.μ], label="Prior μ")
plot!(Normal(mean(μ_summary), mean(σ_summary)), label="Posterior")
vline!([mean(μ_summary)], label="Posterior μ")

```

![n_10](https://global.discourse-cdn.com/julialang/original/3X/5/4/54bb2251561e84faa3acc2250de7385602e62352.png)

If we had 100 scores from this year’s class rather than 10:

![n_100](https://global.discourse-cdn.com/julialang/original/3X/a/1/a137ae49a346253f65d404ce61c2765afede1215.png)

As you can see, the credible interval is much smaller when we have more data.

---

<div class="post-metadata">

**Author:** ![tim\_pi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim_pi/32/30601_2.png) [@tim\_pi](https://discourse.julialang.org/u/tim_pi)\
**Post date:** [November 10, 2021, 6:48pm UTC](https://discourse.julialang.org/t/learning-statistics-seeking-advice-on-how-to-model-a-2d-size-distribution-of-clothing-items/71152/23 "2021-11-10T18:48:28Z")

</div>

Hi @mthelm85 thank you so very much for putting in the effort to create an example. Much obliged. I’ll be studying this later this evening.

I had some doubts raising questions in this community (being green and all), but am very happy I did, thanks again for raising the bar on being helpful! :happy-camper: Cheers!

---

<div class="post-metadata">

**Author:** ![tim\_pi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim_pi/32/30601_2.png) [@tim\_pi](https://discourse.julialang.org/u/tim_pi)\
**Post date:** [November 11, 2021, 12:38pm UTC](https://discourse.julialang.org/t/learning-statistics-seeking-advice-on-how-to-model-a-2d-size-distribution-of-clothing-items/71152/24 "2021-11-11T12:38:27Z")

</div>

I’ve looked at your example and tried it out. I hope you don’t mind me asking some follow up questions.

1. I’m wondering why a `fit` for the prior was chosen (as opposed to just taking `70` from d\_old). And in the same regard: why only use the `prior.μ` in the model, as opposed to using the fitted σ : ie: `μ ~ Normal(prior.μ,prior.σ)`
2. What’s the reasoning behing choosing a `Uniform(0,50)` distribution for σ in the model? I think I understand why it is uniform, but wouldn’t it be a more accurate model if a different interval was chosen? (ie: 15,35)

There were a few things I didn’t understand at first glance, so I reiterated over different values of `N` to get a better grasp of what what happening. I’ve made a plot to better see what was going on (TIL). Maybe it’s useful for other people following this discussion-thread as well.

NOTE: The fitted prior was constant in every iteration, the `scores_new` were resampled.

# Plot showing different `scores_new` samples with `N=10`

![ciplot-10](https://global.discourse-cdn.com/julialang/original/3X/5/c/5c7db7009f6d406b893f49503767d38e4632afb8.gif)

# Plot showing different `scores_new` samples with `N=100`

![ciplot-100](https://global.discourse-cdn.com/julialang/original/3X/5/c/5c4bbce4ed7c6c36ee824ffb5db62c4852bf58a6.gif)

# Plot showing different `scores_new` samples with `N=1000`

![ciplot-1000](https://global.discourse-cdn.com/julialang/original/3X/e/6/e6db0187033cdcbd4da7a2af46707946b9be8fbc.gif)

One final thing I was wondering about looking at the graphs is how to compute when `N` is sufficiently high (ie: wigglyness is low), would you only take `CI` under consideration or are there other measures one can use.

Thank you again for showing me all this and taking your time to provide examples, I understand you might have better things to do than answering my questions, no pressure 😉

---

<div class="post-metadata">

**Author:** ![mthelm85](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mthelm85/32/224164_2.png) [@mthelm85](https://discourse.julialang.org/u/mthelm85)\
**Post date:** [November 11, 2021, 1:35pm UTC](https://discourse.julialang.org/t/learning-statistics-seeking-advice-on-how-to-model-a-2d-size-distribution-of-clothing-items/71152/25 "2021-11-11T13:35:37Z")

</div>

Those are _excellent_ questions and the answer to both of them is mostly that the example was thrown together quickly/sloppily 😆

> [@tim\_pi](#):
>
> I’m wondering why a `fit` for the prior was chosen (as opposed to just taking `70` from d\_old). And in the same regard: why only use the `prior.μ` in the model, as opposed to using the fitted σ : ie: `μ ~ Normal(prior.μ,prior.σ)`

In this case, `d_old` and `d_new` are the true underlying distributions which I almost never know in practice. I generated samples from `d_old` and sort of assumed that would be the starting point in a real-world scenario (you observe some data) but then you might fit a distribution to that data. In this example, of course, there was no reason to do that since we have the actual distribution. In the model, we are saying we want to draw the parameter `μ` from a normal distribution that has the mean `μ` and, in this case, an arbitrarily chosen standard deviation. I think the choice of standard deviation in the line `μ ~ Normal(prior.μ, 25)` is a reflection of how confident you are that the mean for the new distribution could/should be close to the mean of the old distribution.

> [@tim\_pi](#):
>
> What’s the reasoning behing choosing a `Uniform(0,50)` distribution for σ in the model? I think I understand why it is uniform, but wouldn’t it be a more accurate model if a different interval was chosen? (ie: 15,35)

This was arbitrary. The interpretation of this choice is that I’m saying all values between 0 and 50 are equally likely, i.e., it’s expressing a high degree of uncertainty regarding `σ`. If you wanted to express total uncertainty, you could do `Uniform(0, 100)`, but in this case we can compute the standard deviation of our observed data, so maybe we could have done `σ ~ truncated(Normal(prior.σ, 10), 0, 50)`, or something like that.

Finally, I would highly recommend you check up on all of this (maybe with some of the resources I linked above) because I am nowhere near an expert/authority on this topic, but I’m glad to help out! There are some legitimate subject matter experts that participate regularly in this forum, so if we’re lucky one of them will chime in 😉

---

<div class="post-metadata">

**Author:** ![tim\_pi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim_pi/32/30601_2.png) [@tim\_pi](https://discourse.julialang.org/u/tim_pi)\
**Post date:** [November 11, 2021, 1:47pm UTC](https://discourse.julialang.org/t/learning-statistics-seeking-advice-on-how-to-model-a-2d-size-distribution-of-clothing-items/71152/26 "2021-11-11T13:47:13Z")

</div>

Happy to report that your answers to my questions were in line with what I was intuitively thinking. Nonetheless the example you provided gave me a lot of insight! So thank you once more.

I’ve been wondering how to clip distributions and didn’t know about `truncated`. Nice to see it show up here!

I’ll be checking those links out as well.

[Previous page](https://discourse.julialang.org/t/learning-statistics-seeking-advice-on-how-to-model-a-2d-size-distribution-of-clothing-items/71152.md?page=1)