Yhyjyk

Given a variable $y_i$, normally distributed with 0 mean and $σ_y$ standard deviation

$y_i$ ~ NormalDistribution[0,$σ_y$ ]

I want to obtain with Mathematica:

1. The distribution of:
   $x = bary = frac sum_i=1^ny_in$




2. The joint distribution of $ (x,y_i )$

Thank you for your helpful comments

I don't know how to get Mathematica to get the joint distribution explicitly for a general value of $n$ but here is how one can easily see the pattern to figure out the general solution.

First the distribution of the mean:

marginalDistribution = TransformedDistribution[Sum[y[i], i, n]/n, 
 Table[y[i] [Distributed] NormalDistribution[0, [Sigma]], i, n], 
 Assumptions -> [Sigma] > 0]
#, marginalDistribution/.n-># &/@Range[2,10]

$$
beginarraycc
2 & textNormalDistributionleft[0,fracsigma sqrt2right] \
3 & textNormalDistributionleft[0,fracsigma sqrt3right] \
4 & textNormalDistributionleft[0,fracsigma 2right] \
5 & textNormalDistributionleft[0,fracsigma sqrt5right] \
6 & textNormalDistributionleft[0,fracsigma sqrt6right] \
7 & textNormalDistributionleft[0,fracsigma sqrt7right] \
8 & textNormalDistributionleft[0,fracsigma 2 sqrt2right] \
9 & textNormalDistributionleft[0,fracsigma 3right] \
10 & textNormalDistributionleft[0,fracsigma sqrt10right] \
endarray
$$

So we see that the marginal distribution of $bary$ is

NormalDistribution[0, σ/Sqrt[n]]

The joint distribution of $bary$ and, say, $y_1$ is given by

jointDistribution = TransformedDistribution[y[1], Sum[y[i], i, n]/n, 
 Table[y[i] [Distributed] NormalDistribution[0, [Sigma]], i, n]]
#, jointDistribution /. n -> # & /@ Range[2, 10] // TableForm

$$
beginarraycc
2 & textMultinormalDistributionleft[0,0,left(
beginarraycc
sigma ^2 & fracsigma ^22 \
fracsigma ^22 & fracsigma ^22 \
endarray
right)right] \
3 & textMultinormalDistributionleft[0,0,left(
beginarraycc
sigma ^2 & fracsigma ^23 \
fracsigma ^23 & fracsigma ^23 \
endarray
right)right] \
4 & textMultinormalDistributionleft[0,0,left(
beginarraycc
sigma ^2 & fracsigma ^24 \
fracsigma ^24 & fracsigma ^24 \
endarray
right)right] \
5 & textMultinormalDistributionleft[0,0,left(
beginarraycc
sigma ^2 & fracsigma ^25 \
fracsigma ^25 & fracsigma ^25 \
endarray
right)right] \
6 & textMultinormalDistributionleft[0,0,left(
beginarraycc
sigma ^2 & fracsigma ^26 \
fracsigma ^26 & fracsigma ^26 \
endarray
right)right] \
7 & textMultinormalDistributionleft[0,0,left(
beginarraycc
sigma ^2 & fracsigma ^27 \
fracsigma ^27 & fracsigma ^27 \
endarray
right)right] \
8 & textMultinormalDistributionleft[0,0,left(
beginarraycc
sigma ^2 & fracsigma ^28 \
fracsigma ^28 & fracsigma ^28 \
endarray
right)right] \
9 & textMultinormalDistributionleft[0,0,left(
beginarraycc
sigma ^2 & fracsigma ^29 \
fracsigma ^29 & fracsigma ^29 \
endarray
right)right] \
10 & textMultinormalDistributionleft[0,0,left(
beginarraycc
sigma ^2 & fracsigma ^210 \
fracsigma ^210 & fracsigma ^210 \
endarray
right)right] \
endarray
$$

So the general distribution is a multivariate normal

MultinormalDistribution[0, 0, σ^2, σ^2/n, σ^2/n, σ^2/n]

The general form of the joint density function can then be found with

FullSimplify[PDF[MultinormalDistribution[0, 0, σ^2, σ^2/n, σ^2/n, σ^2/n], y, ybar],
 Assumptions -> σ > 0, n > 1]

$$fracn e^-fracn left(n textybar^2+y^2-2 y textybarright)2 (n-1) sigma ^22 pi sqrtn-1 sigma ^2$$

Here is the distribution of $x=overliney$ (Part I of your question):

n = 5; (* for example *)
a = Table[y[k] [Distributed] NormalDistribution[0, [Sigma]], k, 1, n]; 
TransformedDistribution[Sum[y[k]/n, k, 5], a]

The result is

NormalDistribution[0, Abs[[Sigma]]/Sqrt[5]]

UPDATE

Okay, here is how to do it with $n$ a variable:

a[n_] := Table[y[k] [Distributed] NormalDistribution[0, [Sigma]], k, 1, n]; 
p[n_] := TransformedDistribution[Sum[y[k]/n, k, n], a[n]];

Now

x [Distributed] p[5] (* n=5, for example *)

Again, the result is

x [Distributed] NormalDistribution[0, Abs[[Sigma]]/Sqrt[5]]

just modified @mjw's answer,

n = 100;(*for example*)ClearAll[y]; 
a = Table[y[k] [Distributed] NormalDistribution[0, [Sigma]], k, 1, n];
meanDist = TransformedDistribution[Sum[y[k]/100, k, 100], a]

JointDistribution can be composed by ProductDistribution,
if these random variables are independent.

if not,you have to use Copula

joint = ProductDistribution[meanDist, 
 Last@*List @@ Part[a, 1]] /. [Sigma] -> 1;
RandomVariate[joint, 100] // Histogram3D

joint = ProductDistribution[meanDist, 
 Last@*List @@ Part[a, 1]] /. [Sigma] -> 1;
m1 = RandomVariate[meanDist /. [Sigma] -> 1, 100000];
m2 = RandomVariate[
 Last@*List @@ Part[a, 1] /. [Sigma] -> 1, 100000];
Correlation[Thread[List[m1, m2]]]
ListPlot[Thread[List[m1, m2]]]

=>

1., -0.00256777, -0.00256777, 1.

I'm not sure about correlation,but it's okay.