Я пишу имитацию Монте-Карло фотонов, поступающих на детектор. Общее количество обнаруженных фотонов N_detected
соответствует распределению Пуассона (которое я генерирую с использованием scipy.stats.poisson
), а время обнаружения каждого отдельного фотона соответствует заданной функции распределения вероятности (PDF). Чтобы сгенерировать время обнаружения N_detected
фотонов, я генерирую N_detected
случайные числа от 0 до 1 с помощью функции numpy.random.random()
и использую метод обратного преобразования .
Я должен запустить симуляцию несколько раз. Таким образом, чтобы избежать многократного повторения, я бы хотел выполнять каждую симуляцию одновременно, используя массивы. В качестве окончательного результата я хотел бы получить двумерный массив из N_sim
столбцов, где каждый столбец соответствует разной имитации и содержит сгенерированные времена обнаружения. Проблема в том, что каждая симуляция производит разное количество фотонов (поскольку оно случайное), и я не нахожу способ создать двумерный массив со столбцами разной длины. столбцы одинаковой длины (самый большой N_deteceted
) и заполните ненужные элементы NaN
, а остальные случайными числами, которые мне нужны, но я не знаю, смогу ли я это сделать.
Этоэто лучшее, что я смог сделать на данный момент:
import numpy as np
import numpy.ma as ma
from scipy.stats import poisson
beta = 0.0048 # Parameter of the PDF
""" Inverse of the time CDF (needed for the Inverse transform method) """
def inverse_time_cdf(x):
return -np.log( (np.exp(-1000*beta)-1)*x + 1 )/beta
""" Simulation of N_sim experiments through the inverse transfrom method """
T = 1000 # Duration of the experiment [s]
rate = 3.945 # [events/s]
lam = rate*T # Poisson distribution parameter
def photon_arrival_simulation(N_sim):
N_detection = poisson.rvs(lam, size = N_sim) # Number of detections in each experiment
t = np.array([inverse_time_cdf(np.random.random(N)) for N in N_detection])
return N_detection, t
Если возможно, я хотел бы избежать цикла, используемого при понимании списка функции photon_arrival_simulation()
, а также получить двумерный массив вместомассив массивов (поскольку я не могу выполнять операции с массивами, такие как взятие log
для массива массивов).
Я не знаю, должен ли я публиковать этот вопрос здесь или в Physics Stack Exchange , но заранее всем спасибо.