Hcfyk

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 = bar{y} = frac {sum_{i=1}^ny_i}{n}$




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

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]

$$
begin{array}{cc}
2 & text{NormalDistribution}left[0,frac{sigma }{sqrt{2}}right] \
3 & text{NormalDistribution}left[0,frac{sigma }{sqrt{3}}right] \
4 & text{NormalDistribution}left[0,frac{sigma }{2}right] \
5 & text{NormalDistribution}left[0,frac{sigma }{sqrt{5}}right] \
6 & text{NormalDistribution}left[0,frac{sigma }{sqrt{6}}right] \
7 & text{NormalDistribution}left[0,frac{sigma }{sqrt{7}}right] \
8 & text{NormalDistribution}left[0,frac{sigma }{2 sqrt{2}}right] \
9 & text{NormalDistribution}left[0,frac{sigma }{3}right] \
10 & text{NormalDistribution}left[0,frac{sigma }{sqrt{10}}right] \
end{array}
$$

So we see that the marginal distribution of $bar{y}$ is

NormalDistribution[0, σ/Sqrt[n]]

The joint distribution of $bar{y}$ 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

$$
begin{array}{cc}
2 & text{MultinormalDistribution}left[{0,0},left(
begin{array}{cc}
sigma ^2 & frac{sigma ^2}{2} \
frac{sigma ^2}{2} & frac{sigma ^2}{2} \
end{array}
right)right] \
3 & text{MultinormalDistribution}left[{0,0},left(
begin{array}{cc}
sigma ^2 & frac{sigma ^2}{3} \
frac{sigma ^2}{3} & frac{sigma ^2}{3} \
end{array}
right)right] \
4 & text{MultinormalDistribution}left[{0,0},left(
begin{array}{cc}
sigma ^2 & frac{sigma ^2}{4} \
frac{sigma ^2}{4} & frac{sigma ^2}{4} \
end{array}
right)right] \
5 & text{MultinormalDistribution}left[{0,0},left(
begin{array}{cc}
sigma ^2 & frac{sigma ^2}{5} \
frac{sigma ^2}{5} & frac{sigma ^2}{5} \
end{array}
right)right] \
6 & text{MultinormalDistribution}left[{0,0},left(
begin{array}{cc}
sigma ^2 & frac{sigma ^2}{6} \
frac{sigma ^2}{6} & frac{sigma ^2}{6} \
end{array}
right)right] \
7 & text{MultinormalDistribution}left[{0,0},left(
begin{array}{cc}
sigma ^2 & frac{sigma ^2}{7} \
frac{sigma ^2}{7} & frac{sigma ^2}{7} \
end{array}
right)right] \
8 & text{MultinormalDistribution}left[{0,0},left(
begin{array}{cc}
sigma ^2 & frac{sigma ^2}{8} \
frac{sigma ^2}{8} & frac{sigma ^2}{8} \
end{array}
right)right] \
9 & text{MultinormalDistribution}left[{0,0},left(
begin{array}{cc}
sigma ^2 & frac{sigma ^2}{9} \
frac{sigma ^2}{9} & frac{sigma ^2}{9} \
end{array}
right)right] \
10 & text{MultinormalDistribution}left[{0,0},left(
begin{array}{cc}
sigma ^2 & frac{sigma ^2}{10} \
frac{sigma ^2}{10} & frac{sigma ^2}{10} \
end{array}
right)right] \
end{array}
$$

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}]

$$frac{n e^{-frac{n left(n text{ybar}^2+y^2-2 y text{ybar}right)}{2 (n-1) sigma ^2}}}{2 pi sqrt{n-1} sigma ^2}$$

Here is the distribution of $x=overline{y}$ (Part I of your question):

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 *)

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.