Как я могу использовать метод Монте-Карло, чтобы найти неопределенности значения?

Question

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

У меня есть около 5 значений, одно времяи у меня есть случайные неопределенности (ошибки) для каждого из этих значений.Таким образом, масса равна (10 + - 0,1) кг, где погрешность составляет 0,1 кг

Как я на самом деле могу найти распределение измерений, если я, например, выполнил этот эксперимент 5000 раз?

Я знаю, что могу сделать 2 массива ошибок и, возможно, поместить их в функцию.Но что мне тогда делать?Должен ли я положить ошибки в уравнение, а затем добавить ответ в массивы, а затем положить измененные значения массива в уравнение и повторить это тысячу раз.Или я действительно вычисляю реальное значение и добавляю его в массив.

Пожалуйста, не могли бы вы помочь мне понять это.

Редактировать:

У меня проблема в основном в сфере плотности ds, которая падает на расстояние l во времени tчерез жидкость плотности dl это вписывается в уравнение для вязкости, и мне нужно найти распределение измерений вязкости.

Уравнение не должно иметь значения при любом уравнении, которое я могу использовать.метод, подобный этому, чтобы найти распределение измерений.Погода Я бросаю мяч в окно или что-то еще.

Answer 1

Хорошо, давайте попробуем на простом примере - у вас есть пневматический пистолет, который стреляет шариками с массой m и скоростью v.Вы должны измерить кинетическую энергию

E = m * v ² / 2

Существует распределение скорости - гауссовское со средним значением 10 и стандартным отклонением 1.

Существует распределение масс - но мы не можем сделать гауссовский, давайте предположим, что оно усечено нормально, с нижним пределом 1, так что нет отрицательных значений, с loc, равным 5, и масштабом, равным 3.

Итак, что мы будем делать - скорость выборки, масса образца, использовать их, чтобы найти кинетическую энергию, сделать это несколько раз, построить распределение энергии, получить среднее значение, получить стандартное отклонение, нарисовать графики и т. Д.

Простой код Python

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import truncnorm

def sampleMass(n, low, high, mu, sigma):
    """
    Sample n mass values from truncated normal
    """
    tn = truncnorm(low, high, loc=mu, scale=sigma)
    return tn.rvs(n)

def sampleVelocity(n, mu, sigma):
    return np.random.normal(loc = mu, scale = sigma, size=n)


mass_low  = 1.
mass_high = 1000.
mass_mu   = 5.
mass_sigma = 3.0

vel_mu    = 10.0
vel_sigma = 1.0

nof_trials = 100000

mass = sampleMass(nof_trials, mass_low, mass_high, mass_mu, mass_sigma) # get samples of mass
vel  = sampleVelocity(nof_trials, vel_mu, vel_sigma) # get samples of velocity

kinenergy = 0.5 * mass * vel*vel # distribution of kinetic energy

print("Mean value and stddev of the final distribution")
print(np.mean(kinenergy))
print(np.std(kinenergy))

print("Min/max values of the final distribution")
print(np.min(kinenergy))
print(np.max(kinenergy))


# print histogram of the distribution
n, bins, patches = plt.hist(kinenergy, 100, density=True, facecolor='green', alpha=0.75)

plt.xlabel('Energy')
plt.ylabel('Probability')
plt.title('Kinetic energy distribution')
plt.grid(True)

plt.show()

с выводом типа

Mean value and stddev of the final distribution
483.8162951263243
118.34049421853899
Min/max values of the final distribution
128.86671038372
1391.400187563612

Answer 2

Basic Monte Carlo очень прост.Следующее может начать:

import random,statistics,math

#The following function generates a
#random observation of f(x) where
#x is a vector of independent normal variables
#whose means are given by the vector mus
#and whose standard deviations are given by sigmas

def sample(f,mus,sigmas):
    x = (random.gauss(m,s) for m,s in zip(mus,sigmas))
    return f(*x)

#do n times, returning the sample mean and standard deviation:

def monte_carlo(f,mus,sigmas,n):
    samples = [sample(f,mus,sigmas) for _ in range(n)]
    return (statistics.mean(samples), statistics.stdev(samples))

#for testing purposes:

def V(r,h):
    return math.pi*r**2*h

print(monte_carlo(V,(2,4),(0.02, 0.01),1000))

С выводом:

(50.2497301631037, 1.0215188736786902)