2
$\begingroup$

I am trying to generate a synthetic earthquake database where the number of events ($N$) with magnitude ($M$) in the range $[M, M+\delta_M]$ follows: $\log_{10}(N) = a - bM$ where $a$ and $b$ are constants.

I am trying to do this in Python using the "random" module. I know I can (or at least I think I can - as I haven't tried it) use random.expovariate but I thought I could use random.random with a transformation like:

-math.log10(random.random()))

The random function generates numbers between 0 & 1.

I ran this for 2,000,000 samples which I then binned into 0.1 bins and plotted on a log scale.

enter image description here

I'm not worried about the variation above x=4.5. This is due to small number of points and natural randomness. What I am asking about is the very small (at this scale) variation for the points near x=0 - particularly the point at x=0 which appears to be lower than the fitted (by eye) line.

generated data

Given the large number of events in this bin the expected error should be very low. When I look at the difference between the generated and theoretical values I see that points below x ~2.5 start to deviate from linearity.

[edit] Have now changed the red line to represent what my simulation should be generating. The following plot shows the deviation between my simulation and the theoretical points.

diff between observed and theory

To try and work out where my problem lies I wrote code that generated numbers from 0 to 1 with a very fine step. Each number then went through the function noted above. The result was a purely linear result that exactly matched the theoretical values. So my transformation function is fine - so I'm thinking it's just something to do with the random number generator?

Would be grateful for any pointers. Thanks.

$\endgroup$
4
  • $\begingroup$ (1) What exactly does the red line represent? (2) On what basis have you determined the values "start to deviate from linearity"? The fit looks great. $\endgroup$
    – whuber
    Nov 2, 2020 at 23:41
  • $\begingroup$ The red line has the slope of the "b" in the formula and the "a" adjusted (by eye) to fit most of the points. As x decreases, the number of simulations rises exponentially so the fit (to the straight line) should get better and better.The final point (x=0) is clearly off the fitted line. I could add another plot showing the deviation but its shape depends on the choice of the a&b-values. $\endgroup$
    – RustyC
    Nov 3, 2020 at 7:10
  • $\begingroup$ Fitting by eye is not a valid way to assess data. The red line should be drawn to represent the distribution you are simulating from. $\endgroup$
    – whuber
    Nov 3, 2020 at 14:05
  • $\begingroup$ @whuber - Agreed. So have modified the red line to represent the theoretical distribution. $\endgroup$
    – RustyC
    Nov 4, 2020 at 0:07

1 Answer 1

1
$\begingroup$

The plot of residuals (deviations between the observed and theoretical frequencies) was inspired, because it shows what is going wrong.

You have generated $n$ independent random numbers according to some distribution function $F$ (defined as $F(x)$ is the chance each number is $x$ or less). You then binned them into the intervals $(-\infty,0.1],$ $(0.1,0.2],$ and so on, and tallied the counts in each bin. Your plots indicate the bin locations (horizontal axes) and the counts (or related values) on the vertical axes.

It is important to get both parts of the comparison correct: you need to generate the random values in a way that follows the intended distribution $F$ and you need to compute the theoretical counts in the bins.

Your graphics suggest you haven't computed the expected counts correctly, so I will discuss this in detail.

The expected count of values lying in some bin $(a,b]$ is the expected count of all values $b$ or smaller, minus the expected count of all values $a$ or smaller. This will be proportional to $F(b) - F(a).$ In the present case you have generated uniform random values $U$ in the interval $(0,1]$ and taken their common logarithms $X = -\log_{10}U.$ Thus, for any positive number $x,$

$$F(x) = \Pr(X \le x) = \Pr(-\log_{10}U \le x) = \Pr(U \ge 10^{-x}) = 1 - 10^{-x}.$$

Consequently

$$\Pr(X \in (a,b]) = F(b) - F(a) = (1 - 10^{-b}) - (1 - 10^{-a}) = 10^{-a} - 10^{-b}.$$

Because each number $x_i$ in your simulation is generated independently, the expected count of values in a bin $(a,b]$ (with $a \ge 0$) therefore is

$$\begin{aligned} e(a,b) &= E\left[\# x_i\mid a \lt x_i \le b\right] = n\Pr(a \le X \le b) = n(F(b)-F(a))\\ & = n(10^{-a} - 10^{-b}). \end{aligned}$$

For $n=2.2\times 10^6$ and the bins you use, these expected counts are

$$(452477.9, 359416.0, 285494.2, 226776.1, \ldots).$$

They drop below $1$ beginning with the 58th bin $(5.7, 5.8].$

I ran this simulation (in R). Although my random numbers must differ from yours, my results should be qualitatively the same. Here they are, in graphics.

First, here are the counts shown as vertical bars (erected at the right endpoint of each bin) and the expected values shown as a connected red curve:

Figure 1

Next, because on a semi-log plot the expected values all lie on a common line, here is the same plot with a logarithmic vertical axis:

Figure 2

It's hard to tell how good the agreement is because the logarithms compress variation among the highest counts so much. Your residual plot of differences between the counts and their expectations gives a clearer comparison:

Figure 3

(To make this plot I combined all the counts in the rightmost bins because none of them was expected to have more than five values each.)

Sure enough, there's much more variation at the left. But this is to be expected, because (to an incredibly good approximation) the variation in a random count is proportional to the square root of its expectation. This variation is called the standard error of the count. Consequently, to make a fair comparison of these residuals, we "standardize" them by dividing by their standard errors and redo the plot:

Figure 4

Now the residuals randomly oscillate among positive and negative values -- some counts are a little high, some a little low -- but they show no tendency to be larger in any sequence of bins, nor do they show any tendency to have runs of positive or negative deviations. This is strong evidence that the simulation conforms with the expected distribution.

One simple way to assess this evidence is called a chi-squared test. In this application we simply sum the squares of all the $m$ residuals shown in the plot and compare that value to a chi-squared distribution with parameter $m-1.$ In the simulation I have shown, this sum of squares is $42.23$ and $m=51.$ In a $\chi^2(51-1)$ distribution, $22.5\%$ of the probability is less than $42.23$ and the remaining $77.5\%$ is greater: this places the chi-squared statistic squarely in the middle of the distribution. That means there's nothing unusual about the value of $42.23;$ consequently, despite having generated over two million random values, there still is no evidence their distribution differs from the theoretical (that is, intended) distribution.

I suspect that when you perform this analysis with your data, you will arrive at the same conclusion.

$\endgroup$
0

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.