1

I would like to write a python script to generate a uniformly distributed 3D coordinates (e.g., x, y, z) where x, y, and z are float numbers between 0 and 1. For the moment, z can be fixed, thus what I need is a uniform distributed points in a 2D (x-y) plane. I have written a script to do this job and checked both x, and y are uniform numbers. However, I am not sure if these points are uniformly distributed in (x-y) plane.

My code is

1 import matplotlib.pyplot as plt
2 import random
3 import numpy as np
4 import csv
5 nk1=300
6 nk2=300
7 nk3=10
8 kx=[]
9 ky=[]
10 kz=[]
11 for i in range(nk1):
12     for j in range(nk2):
13         for k in range(nk3):
14             xkg1=random.random()
15             xkg2=random.random()
16             xkg3 = float(k)/nk3
17             kx.append(xkg1)
18             ky.append(xkg2)
19             kz.append(xkg3)
20 kx=np.array(kx)
21 count, bins, ignored = plt.hist(kx, normed=True)
22 plt.plot(bins, np.ones_like(bins), linewidth=2, color='r')
23 plt.show()

The plot shows both "kx", and "ky" are uniformly distributed numbers, however, how can I make sure that x-y are uniformly distributed in the 2D plane?

Mad Physicist
  • 107,652
  • 25
  • 181
  • 264
Sky
  • 33
  • 5
  • 1
    Instead of the loops, you could directly use the numpy [`rand`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.rand.html) function: `np.random.rand(nbr_points, nbre_dim)` – xdze2 Aug 22 '18 at 14:47
  • Sorry for the frequent updates to my answer. I'm on mobile and I had to keep checking other stuff. Finished now. – Mad Physicist Aug 22 '18 at 15:27
  • Could you please remove the line numbers from your code? They make it very difficult to copy and paste into an editor or console. – Mad Physicist Aug 22 '18 at 16:19

2 Answers2

2

Just as you used np.histogram1 to check uniformity in 1D, you can use np.histogram2d to do the same thing in 2D, and np.histogramdd in 3D+.

To see an example, let's first fix your loops by making them go away:

kx = np.random.random(nk1 * nk2 * nk3)
ky = np.random.random(nk1 * nk2 * nk3)
kz = np.tile(np.arange(nk3) / nk3, n1 * n2)

hist2d, *_ = np.histogram2d(kx, ky, range=[[0, 1], [0, 1]])

The range parameter ensures that you are binning over [0, 1) in each direction, not over the actual min and max if your data, no matter how close it may be.

Now it's entirely up to you how to visualize the 100 data points in hist2d. One simple way would be to just ravel it and do a bar chart like you did for the 1D case:

plt.bar(np.arange(hist2d.size), hist2d.ravel())
plt.plot([0, hist2d.size], [nk1 * nk2 * nk3 / hist2d.size] * 2)

Another simple way would be to do a heat map:

plt.imshow(hist2d, interpolation='nearest', cmap='hot')

This is actually not as useful as the bar chart, and doesn't generalize to higher dimensions as well.

Your best bet is probably just checking the standard deviation of the raw data.

1 Or rather plt.hist did for you under the hood.

Mad Physicist
  • 107,652
  • 25
  • 181
  • 264
  • @Sky. I probably do. This whole thing was written with no testing on mobile. I'd love to hear what it is and fix it :) – Mad Physicist Aug 22 '18 at 16:01
  • Thank you very much! It is a bit complicated to me about the concept of "histogram2d". I think you have some typo in "hist2d, *_ = np.histogram2d(kx, ky, range=[[0, 1], [0, 1]])" because this function should return 3 values. There is another error in "plt.bar(hist2d.ravel())" because this function needs two inputs. I am now trying to digest your suggestion now. – Sky Aug 22 '18 at 16:14
  • @Sky. If you are using python 3, `*_` will eat up and discard all the additional outputs. Actually, it'll just put them into a two element list called `_`. – Mad Physicist Aug 22 '18 at 16:15
  • @Sky. You're completely right about bar requiring an x input. I've fixed that now – Mad Physicist Aug 22 '18 at 16:18
  • I finally understand the logic of your suggestions. But I modified a bit your script in order to realize the confirmation of uniform distribution in 2D. – Sky Aug 25 '18 at 13:04
  • @Sky. That's excellent if you understand enough to modify it meaningfully, I feel like my job here is done. If the answer helped you, feel free to select it by clicking on the check mark next to it. – Mad Physicist Aug 25 '18 at 13:27
0

With the help of @Mad Physicist, I finally found the way to verify the uniform distribution of random numbers in 2D. Here I post my script, and explain the details:

 1 import numpy as np
 2 import random
 3 import matplotlib.pyplot as plt
 4 import matplotlib
 5 fig = plt.figure()
 6 ax1 = fig.add_subplot(211)
 7 ax2 = fig.add_subplot(212)
 8 nk=100
 9 nk=100
10 nk=1
11 kx1=[]
12 ky1=[]
13 kz1=[]
14 for i in range(nk1):
15     for j in range(nk2):
16         for k in range(nk3):
17             xkg =r andom.random()
18             ykg = random.random()
19             zkg = float(k)/nk3
20             kx.append(xkg)
21             ky.append(ykg)
22             kz.append(zkg)
23 kx=np.array(kx)
24 ky=np.array(ky)
25 kz=np.array(kz)
26 xedges, yedges = np.linspace(0, 1, 6), np.linspace(0, 1, 6)
27 ## count the number of poins in the region definded by (xedges[i], 
    xedges[i+1])
28 ## and (yedges[i], xedges[y+1]). There are in total 10*10 2D 
    squares. 
29 hist, xedges, yedges = np.histogram2d(kx, ky, (xedges, yedges))
30 xidx = np.clip(np.digitize(kx, xedges) - 1, 0, hist.shape[0] - 1)
31 yidx = np.clip(np.digitize(ky, yedges) - 1, 0, hist.shape[1] - 1)
32 ax1.bar(np.arange(hist.size),hist.ravel())
33 ax1.plot([0,hist.size], [nk1 * nk2 * nk3 / hist.size] * 2)
34 c = hist[xidx, yidx]
35 new = ax2.scatter(kx, ky, c=c, cmap='jet') 
36 cax, _ = matplotlib.colorbar.make_axes(ax2)
37 cbar = matplotlib.colorbar.ColorbarBase(cax, cmap='jet')
38 ax2.grid(True)
39 plt.show()
Sky
  • 33
  • 5