I'm trying to choose random positions in a cube, weigthed by a function F(R), where R is the distance to the center[0,0,0].
My idea is to create a discrete list of points in the cube, calculate F(R) for each of theses points and then use random.choices to pick N points. To start out F is simply 1/(R^0.5) so I would expect the majority of the N points to be close to center.
import numpy as np
import matplotlib.pyplot as plt
import random
grid1d = np.linspace(-100,100,100)
grid = np.array([[x,y,z] for x in grid1d for y in grid1d for z in grid1d])
def density(pos):
return np.sqrt(pos[0]**2 + pos[1]**2 + pos[2]**2)** -0.5
weights = np.array([density(pos) for pos in grid])
weights /= np.sum(weights)
positions = random.choices(grid,weights = weights, k = 1000000)
positions_r = [np.sqrt(pos[0]**2 + pos[1]**2 + pos[2]**2) for pos in positions]
plt.hist(positions_r)
plt.show()
max_weight_index = np.argpartition(weights,-10)[-10:]
print(grid[max_weight_index])
'''[[ 1.01010101 -1.01010101 -3.03030303]
[-1.01010101 -1.01010101 3.03030303]
[-1.01010101 -1.01010101 1.01010101]
[-1.01010101 -1.01010101 -1.01010101]
[ 1.01010101 1.01010101 -1.01010101]
[ 1.01010101 -1.01010101 -1.01010101]
[ 1.01010101 -1.01010101 1.01010101]
[-1.01010101 1.01010101 -1.01010101]
[ 1.01010101 1.01010101 1.01010101]
[-1.01010101 1.01010101 1.01010101]]'''
I have no real statistics background so there is a good chance I messed up something while calculating the weights, but I would still expect the majority of the points to be near the center. And at least the highest weights are associated with the points closest to the center. But the histogram shows I'm completly off.
Now if I change F to something steeper like 1/(R^3) the result looks more like what I expected.
I would've expected the first histogram to look somewhat similar just not as steep. Do I misunderstand random.choices? Is my expectation wrong? Is my calculation for R wrong? Do I really have to read a stats book? It feels like I just made a simple mistake or typo somewhere but I just can't find it.
Aucun commentaire:
Enregistrer un commentaire