Модификация кода для улучшения симуляции

Question

Как реализовать реализацию кода, предполагая, что количество баллов составляет 10000 и 100000 повторений.Чтобы моделирование выполнялось правильно?

using Statistics, Random, DataFrames, DataFramesMeta, CSV, PyPlot

macro assertprob(x)
msg = string("wrong ", x, ": ")
:(0≤$(esc(x))≤1 || throw(ArgumentError(string($msg,$(esc(x))))))
end

function simulate(p::Float64, q::Float64)
@assertprob p
@assertprob q
@assertprob p+q

t = 0
while true
    t += 1
    r = rand()
    r < p && return t
    r < p + q && return missing
end
end

function getpoint()
while true
    p, q = rand(), rand()
    p + q ≤ 1 && return (p, q)
end
end

function runsim(points=10^3, reps=10^3)
df = DataFrame(p=Float64[], q=Float64[], rep=Int[],
               sim=Union{Int,Missing}[])
for i in 1:points
    p, q = getpoint()
    for j in 1:reps
        push!(df, (p, q, j, simulate(p, q)))
    end
end
df
end

function analyzesim(df)
@linq df |>
by([:p, :q], msim=mean(collect(skipmissing(:sim)))) |>
transform(mtheory=1 ./ (:p .+ :q)) |>
with(scatter(:msim, :mtheory))
end

Random.seed!(1)
df = runsim()
CSV.write("results.txt", df)
analyzesim(df)

У кого-нибудь есть идея?Заранее спасибо за помощь?

Answer 1

Хотя я полностью согласен с приведенными выше комментариями (этот форум не является местом для обзоров кода, а ваш вопрос плохо сформулирован), для Снежков это тяжелая жизнь.

Я постараюсь ответить на вопрос«реализовать (псевдо?) код, чтобы симуляция выполнялась правильно», но поскольку код выполняется, мы можем рассмотреть другие возможные проблемы.Тем не менее, все зависит от моего понимания того, какое значение вы пытаетесь Монте-Карло

1. simulate возвращает missing, если p > r < p+q < 1.Если это неверный случай, возможно, повторите эксперимент, вернув пропущенный.Тем более, что вы пропускаете пропуски.
2. simulate возвращает либо missing, либо Float64.Это не является стабильным типом и не является предпочтительным.
3. Аналогично, код динамически увеличивается df.Это также снижение производительности и должно быть предварительно выделено.
4. Я не уверен, что DataFrame действительно нужен.Код может хранить p и q каждый в двух векторах npoints длиной, а результаты t из simulate в матрице reps на npoints затем усреднять по первому измерению (по reps)

Наконец, для ускорения работы можно сделать много математических операций.p и q взяты равномерно из первого квадранта шара нормы L1.Иными словами, вероятность любой пары (p, q) равна 1, если она находится в этом регионе.

Небольшая математика дает предельную плотность и кумулятивные пробные функции как: f(p) = 2(1-p) и F(p) = p*(2-p).(Математика была сделана для нормализации).Таким образом, вы можете нарисовать случайный p, сначала нарисовав u = rand() и решив u = 2p - p^2.Затем нарисуйте случайное число q = (1-p)*rand().

Наконец, если я правильно понимаю симуляцию, вы рассчитываете ожидаемое количество попыток, чтобы получить r = rand() < p, но отбрасываете ожидание, если p < r < p+q.

Тем не менее, код учитывает только то, сколько попыток потребовалось для выхода из области [p+q, 1], а не в том случае, если окончательная выборка была в средней области (p, p+q), поскольку код отбрасывает все missing и не выглядитна сколько повторений были равны missing (в среднем q/(p+q)*reps).Это означает, что данные не отсутствуют случайным образом, а вместо этого отсутствуют с вероятностью q/(p+q), что смещает оценку t в сторону от одного из углов (p, q) поддержки.Кроме того, число не пропущенных выборок и, следовательно, дисперсия оценки t s зависит также от q, что вводит дополнительные соображения для выполнения проверки достоверности соответствия или значения p для ваших данных

Я на самом деле очень рад, что сел с ручкой и бумагой, чтобы обдумать это.Математика для вычисления количества образцов, необходимых для посадки в [0, p], но не [p, p+q] (p/(p+q)^2), была забавной, и мне потребовалось некоторое время, чтобы подумать, почему msim было примерно 1/(p+q).

Надеюсь, это поможет!