6
$\begingroup$

I have tried unsuccessfully to use Mathematica to derive both the expected value (mean) and variance of a particular statistical distribution of interest to me.

The distribution is a discretized version of the Gaussian, and the question is: what is the expected value of $\lfloor X \rfloor$, where $X\sim N(\mu, \sigma)$ and $\lfloor \cdot \rfloor$ denotes the floor function (i.e., in Mathematica, Floor[]) ?

I have tried a direct integration approach (with assumptions) allowing for a general mean and standard deviation:

Integrate[
 Floor[x]*PDF[NormalDistribution[\[Mu], \[Sigma]], 
   x], {x, -\[Infinity], \[Infinity]}, 
 Assumptions -> {\[Mu] \[Element] Reals, \[Sigma] \[Element] 
    PositiveReals}]

and also tried restricting the question to one involving an underlying standard normal distribution:

Integrate[
 Floor[x]*PDF[NormalDistribution[0, 1], 
   x], {x, -\[Infinity], \[Infinity]}]

Equivalently, for the 2nd raw moment:

Integrate[
 (Floor[x])^2 * PDF[NormalDistribution[\[Mu], \[Sigma]], 
   x], {x, -\[Infinity], \[Infinity]}, 
 Assumptions -> {\[Mu] \[Element] Reals, \[Sigma] \[Element] 
    PositiveReals}]

The approach fails in both instances, but the problem is solvable. One solution, for the mean, which uses the Poisson summation formula, is:

$$ \mathbb E\lfloor X \rfloor = \mu - \mathbb E \{X\} = \mu - \frac{1}{2} + \frac{1}{\pi} \sum_{k=1}^\infty e^{-2\pi^2\sigma^2 k^2} \frac{\sin(2\pi k \mu)}{k} $$

and, indeed, it strikes me that using the Poisson summation formula is likely to yield a solution for the question about variance also. However, to use the Poisson summation formula, one needs the Fourier transform of the distribution, but with that I have also been unsuccessful.

What is an appropriate way of using Mathematica to obtain a solution both for the mean, and the variance for my distribution.

Further explanation: I've just seen Bob Hanlon's very instructive numerical-methods solution ... but I'd really like to find a way of reproducing the Poisson summation/Fourier-based analytic solution. One very significant reason for preferring a solution involving the Poisson summation formula (and associated Fourier description) is that it clearly reveals the oscillatory nature of the function, and it is the oscillation that particularly interests me. It's true that I have displayed an analytic description of the mean for the Floor[] function, but there are several other distributions of interest, so I'll like to learn a general Mathematica approach.

$\endgroup$
7
  • 2
    $\begingroup$ I wonder if directly using the probability mass function for the floor of $X$ would help: pmf[x_, \[Mu]_, \[Sigma]_] := CDF[NormalDistribution[\[Mu], \[Sigma]], x + 1] - CDF[NormalDistribution[\[Mu], \[Sigma]], x]. It also seems that if $\mu$ is an integer, then the mean simplifies to $\mu-1/2$ and if $\mu$ and $\sigma$ are integers, then the variance is $\sigma^2+1/12$. $\endgroup$
    – JimB
    Jan 2, 2023 at 2:26
  • $\begingroup$ @Jim8 Yes, I had thought of that, but in terms of obtaining an analytic solution I made no further progress than with the PDF. Numerically, following the solutions offered by user64494 and Bob Hanlon, the mass function almost doubles the time taken as compared with using Integrate[] or NSum[] with the PDF. $\endgroup$ Jan 2, 2023 at 9:50
  • 1
    $\begingroup$ I think your adjustment to $\mu-1/2$ for non-integer values of $\mu$ (which is $ \frac{1}{\pi} \sum _{k=1}^{\infty } e^{-2 \pi ^2 \sigma ^2 k^2} \frac{\sin (2 \pi k \mu )}{k} $) simplifies to $\frac{1}{\pi} e^{-2 \pi ^2 \sigma ^2} \sin (2 \pi \mu )$. $\endgroup$
    – JimB
    Jan 3, 2023 at 4:01
  • $\begingroup$ I was wrong. The simplification is not exact but very, very good for $\sigma>0.3$. $\endgroup$
    – JimB
    Jan 4, 2023 at 1:13
  • $\begingroup$ The approximation just added should make the execution time more tolerable. $\endgroup$
    – JimB
    Jan 4, 2023 at 17:28

3 Answers 3

8
$\begingroup$

Use numeric methods.

Clear["Global`*"]

For the mean,

μ[m_?NumericQ, s_?Positive] := NIntegrate[Floor[x]*
     PDF[NormalDistribution[m, s], x],
    {x, -∞, ∞},
    WorkingPrecision -> 30,
    MinRecursion -> 25,
    MaxRecursion -> 100] //
   N // Quiet

For example,

μ @@@ {{-3, 2}, {0, 1}, {3, 2}}

(* {-3.5, -0.5, 2.5} *)

For the variance,

var[m_?NumericQ, s_?Positive] := 
   NIntegrate[(Floor[x] - μ[m, s])^2*
     PDF[NormalDistribution[m, s], x],
    {x, -∞, ∞},
    WorkingPrecision -> 30,
    MinRecursion -> 25,
    MaxRecursion -> 100] //
   N // Quiet

Examples,

var @@@ {{-3, 2}, {0, 1}, {3, 2}}

(* {4.08333, 1.08333, 4.08333} *)
$\endgroup$
3
  • $\begingroup$ (+1) and best wishes for a happy new year! $\endgroup$
    – bmf
    Jan 2, 2023 at 2:21
  • $\begingroup$ Indeed (+1, Happy New Year), even with my having slightly amended my question. With regard to the numerical solution, could you explain what role the Recursion options are playing here? $\endgroup$ Jan 2, 2023 at 2:26
  • 1
    $\begingroup$ The integrand has numerous discontinuities. Specifying the MinRecursion and MaxRecursion increases the subdivisions to improve the accuracy of the calculations. $\endgroup$
    – Bob Hanlon
    Jan 2, 2023 at 2:32
6
$\begingroup$

A shorter way is as follows.

N[Mean[TransformedDistribution[Floor[x], x \[Distributed] NormalDistribution[-3, 2]]], 3]

-3.50

N[Variance[TransformedDistribution[Floor[x], x \[Distributed] NormalDistribution[-3, 2]]], 3]

4.08

Addition. One more way is as follows.

NSum[n*Integrate[ PDF[NormalDistribution[-3, 2], t], {t, n, n + 1}], {n, -Infinity, Infinity}]

-3.5

NSum[n^2*Integrate[PDF[NormalDistribution[-3, 2], t], {t, n, n + 1}], {n, -Infinity,  Infinity}] -NSum[n*Integrate[PDF[NormalDistribution[-3, 2], t], {t, n, n + 1}], {n, -Infinity, Infinity}]^2

4.083333333333332

$\endgroup$
4
$\begingroup$

This is another numerical approach that seems to be much quicker than the other numerical approaches. This simply starts out acknowledging that you have a discrete distribution over the integers. We have

$$\text{Pr}(\lfloor X \rfloor = x) =\Phi \left(\frac{x+1-\mu}{\sigma }\right)-\Phi \left(\frac{x-\mu }{\sigma }\right)=\frac{1}{2} \text{erf}\left(\frac{\mu -x}{\sqrt{2} \sigma }\right)-\frac{1}{2}\text{erf}\left(\frac{\mu -x-1}{\sqrt{2} \sigma }\right)$$

A function that returns the original mean and standard deviation for $X$ and the estimated mean and variance for $\lfloor X \rfloor$ follows:

stats[μ_, σ_] := Module[{n1, n2, mean, variance},
  n1 = Floor[μ - 6 σ];
  n2 = Ceiling[μ + 6 σ];
  If[IntegerQ[μ], mean = μ - 1/2,
   mean = Sum[x (-(1/2) Erf[(-1 - x + μ)/(Sqrt[2] σ)] + 
     1/2 Erf[(-x + μ)/(Sqrt[2] σ)]), {x, n1, n2}]];
  variance = Sum[(x - mean)^2 (-(1/2) Erf[(-1 - x + μ)/(Sqrt[2] σ)] + 
     1/2 Erf[(-x + μ)/(Sqrt[2] σ)]), {x, n1, n2}];
  {μ, σ, mean // N, variance // N}
]

For example:

stats[-3, 2] // AbsoluteTiming
(* {0.0010298, {-3, 2, -3.5, 4.08333}} *)
stats[0, 1] // AbsoluteTiming
(* {0.0004044, {0, 1, -0.5, 1.08333}} *)
stats[3, 2] // AbsoluteTiming
(* {0.0009363, {3, 2, 2.5, 4.08333}} *)

Note that for $\sigma\geq 1$ even quicker approximations for the mean and variance are, respectively:

$$\mu-1/2$$ $$\sigma^2+1/12$$

$\endgroup$
3
  • $\begingroup$ I know an analytical method is desired and I'd also like to see how Mathematica could do that. $\endgroup$
    – JimB
    Jan 4, 2023 at 17:35
  • $\begingroup$ (+1) for numerics. Definitely the fastest of those so far listed (by a factor of 2) for the mean with $\mu=1/10, \sigma=1/11$,, at least on my ancient machine. $\endgroup$ Jan 5, 2023 at 2:35
  • $\begingroup$ The savings come from the use of the Erf function and that when $\sigma$ is small, there are very few values of $x$ to consider (as the sum of the probabilities of those values of $x$ are almost 1.0). For $\mu=1/10$ and $\sigma=1/11$, the only values of $x$ that really need to be examined are -1 and 0. $\endgroup$
    – JimB
    Jan 5, 2023 at 3:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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