6
$\begingroup$

I'm not well versed in statistics so I'm not sure if my question is worded exactly correctly but basically here's the problem I'm trying to solve: imagine you have two equal sized arrays of size n. Each array is filled with random numbers from 0 to 1. The arrays are then subtracted from each other, and the absolute value taken. The average value of this array is then calculated, and this value then subtracted from 1. This process is repeated for T trials, and the values are plotted with a histogram to visualize the distribution. I would like an equation that gives said distribution as well as the mean and standard deviation for the distribution. I have searched extensively online to try and find an answer to no avail. Below you will find some matlab code that does what I'm talking about for an array of size 8 and 100,000 trials:

clc

clear all

n=8;

T=100000;

H=zeros(1,T);

for i=1:T
    A=rand(1,n);
    B=rand(1,n);
    H(i)=mean(1-abs(A-B));
end

histfit(H,100);


Minimum=min(H)
Maximum=max(H)
Mean=mean(H)
Standard_Deviation=std(H)

One can observe that the mean approaches 2/3 with increasing trials, and that the standard deviation decreases with increasing array size n. For n=8 it is about 8.33. It appears to be normal but not quite.

Any help greatly appreciated, thanks.

$\endgroup$
2
  • $\begingroup$ min max mean and std are supposed to be on their own lines $\endgroup$
    – David Webb
    Oct 22 '14 at 7:18
  • $\begingroup$ The standard deviation at n=8 should be about 0.0833 (in fact, it should be exactly 1/12) rather than 8.33. $\endgroup$
    – Glen_b
    Oct 22 '14 at 12:41
10
$\begingroup$

$A-B$ has a symmetric triangular distribution on $(-1,1)$. It has mean 0 and variance $\frac{1}{6}$.

$|A-B|$ has a $\text{beta}(1,2)$ distribution. It has mean $\frac{1}{3}$ and variance $\frac{1}{18}$.

The distribution of a sum of beta random variables is known for $n=2$ (see 1).

Edit:

Actually, this particular beta is so simple we can just do the integral, I don't know why I didn't jut try it before. Let $Y_i=|A_i-B_i|$. We can compute the density of $Z=Y_1+Y_2$ and you could find the distribution of $1-\frac{1}{2}\sum_{i=1}^2 |A_i-B_i|$ from that.

By simple, direct integration (of the convolution integral),

$f_Z(z) = \begin{cases} \frac{2}{3}z\,(z^2-6z+6) &\mbox{for } 0<z\leq 1 \\ \frac{2}{3}(2-z)^3 & \mbox{for } 1<z\leq 2\\ 0 &\mbox{elsewhere }. \end{cases}$

enter image description here

Then the density of $H_2=1-Z/2$ is just a linear rescaling of that density.

[This very rapidly becomes unwieldy. If I've done this right, a sum of 4 $Y_i$'s will have four pieces each consisting of 7th-order polynomials, and 8 would have eight pieces, each 15th order polynomials. If you've got a nice computer algebra system handy, you could certainly calculate them, but I don't think it's really going to be informative.]

--

For moderate $n$ I don't think even this relatively simple case is not known algebraically - but you could do the convolution numerically, fairly simply.

For large $n$ you could make use of the central limit theorem. The mean and variance of $1-\frac{1}{n}\sum_{i=1}^n |A_i-B_i|$ are straightforward - if I haven't made an error, they're $\frac{2}{3}$ and $\frac{1}{18n}$.

While it's not so accurate for $n=8$, by $n=20$ the normal approximation is pretty good:

enter image description here

===

Edit:

What exactly do you mean when you say you could do the convolution numerically? What exactly would I be convolving?

Let $Y_i=|A_i-B_i|$, so $H=1-\frac{1}{n}\sum_{i=1}^n Y_i$.

Aside from a linear rescaling -- which is trivial -- you need the distribution of the sum of $Y_i$.

That's where the convolution comes in; the density of $W+Z$ is the convolution of their pdfs.

Of course in practice, you don't do the convolution integral. The usual statistical approach is to use MGFs or more generally, characteristic functions, but apart from a matter of sign (again, trivial), a CF is just a Fourier transform (in basically the same sense that MGFs are essentially Laplace transforms).

So numerically, you could use FFTs to take care of the convolution. In fact, it proceeds as follows (I'll use the symbol for a characteristic function, but you can freely think Fourier transform):

$\phi(Y_1+Y_2+...+Y_n)= \phi(Y_1)\times \phi(Y_2)\times ...\times \phi(Y_n)$

$\hspace{2cm} = \phi(Y_1)^n$

Then we just convert back by the inverse transform.

In a program where you operate numerically, the usual approach is to suitably discretize the pdf (into some reasonably large power-of-2 pieces; in some cases $2^{8}$ to $2^{10}$ might be enough, in other cases you might perhaps want $2^{16}$ or more), take the FFT, take it to the nth-power, transform back. If you handle the constants correctly (which should be automatically handled by the inverse FFT anyway, you have at the end a discrete approximation to the pdf of the sum of the $n$ iid random variables.

In practice such operations may be a little tedious to get right the first time, but they tend to be quite fast to run.

===

1: Pham-Gia, T. and Turkkan, N. (1994). Reliability of a standby system with beta component lifelength. IEEE Transactions on Reliability, 71–75.

$\endgroup$
8
  • $\begingroup$ We are on the same wavelength :) $\endgroup$
    – wolfies
    Oct 22 '14 at 12:02
  • $\begingroup$ Oh, hey, I didn't even notice your answer there. Very nice. Yes, it really looks like we are thinking along the same lines. $\endgroup$
    – Glen_b
    Oct 22 '14 at 12:47
  • $\begingroup$ You were first :) (crossed paths while typing them up) $\endgroup$
    – wolfies
    Oct 22 '14 at 12:51
  • $\begingroup$ Thank you for your answer. Unfortunate that there's not an exact symbolic solution. What exactly do you mean when you say you could do the convolution numerically? What exactly would I be convolving? $\endgroup$
    – David Webb
    Oct 22 '14 at 23:29
  • $\begingroup$ See the edit to the answer $\endgroup$
    – Glen_b
    Oct 23 '14 at 0:10
10
$\begingroup$

If $A$ and $B$ are standard Uniform and independent, then $(A-B) \sim Triangular(-1,0,1)$.

Then $Z = |A-B|$ will have pdf $f(z)$:

 f = 2 (1 - z);   domain[f] = {z, 0, 1};

Then, the characteristic function (cf) of the sample mean of $Z$ is $\big(E\big[e^{\large i \frac{t}{n} z}\big] \big)^n$:


(source: tri.org.au)

where I am using the Expect function from the mathStatica package for Mathematica to automate the wurly-curlies.

(i) The EXACT pdf of the sample mean of Z

Although the cf does not appear to have a nice tractable form for symbolic inversion, we can invert it numerically, given any arbitrary value of $n$, to yield the pdf of the sample mean of $Z$. This is done below in the plot (see Blue curve).

(ii) CLT APPROXIMATION of the sample mean of Z

We derived above the pdf of $Z$, namely $f(z)$, where $Z$ has mean $E[Z]= \frac13$ and variance $Var(Z) = \frac{1}{18}$.

Then, by the Central Limit Theorem, the sample mean of $Z$ is asymptotically Normal:

$$\bar Z_n \overset{a}{\sim } N\big(E[Z], \frac{Var(Z)}{n}\big) = N\big(\frac13, \frac{1}{18 n}\big)$$

Compare actual to approximate

We can now compare the easy asymptotic Normal approximation to the exact pdf, for any given value of $n$. The following diagram illustrates the:

  • exact solution (derived by inverting the cf), when $n = 10$: BLUE CURVE
  • Normal approximation (CLT), when $n=10$: RED CURVE


(source: tri.org.au)

Even for reasonably small $n=10$, the Normal approximation works rather nicely.

All that remains is the easy transformation from $\bar Z_n$ to $1-\bar Z_n$ which just involves changing your Normal mean from $\frac13$ to $\frac23$, and all done.

$\endgroup$
2
  • $\begingroup$ Thank you for your answer. What program did you use to numerically determine the pdf? Do you know if this is possible in MATLAB? $\endgroup$
    – David Webb
    Oct 22 '14 at 23:27
  • 1
    $\begingroup$ I used Mathematica to numerically invert the cf, to obtain the pdf. This is done by numerical integration. To see how it is done, see Section 2.4 D of Chapter 2 of Rose/Smith (2002) ... a free download of which is available here: mathstatica.com/book/bookcontents.html I am sure you could use the same techniques in MATLAB. $\endgroup$
    – wolfies
    Oct 23 '14 at 7:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.