Как распространить неопределенность на пользовательскую функцию, используя python

Question

Я хотел бы распространять неопределенность с помощью Python. Это относительно просто для простых функций через пакет неопределенностей. Тем не менее, это не так очевидно, чтобы достичь того же с пользовательской функцией. Далее следует пример того, что я пытаюсь сделать.

import mcerp as err
import numpy as np


def mult_func(x,xm ,a):
    x[x==0.] = 1e-20
    v = (1.-(xm/x)**a) * (x > xm)
    v[np.isnan(v)] = 0.
    return v


def intg(e,f,cut,s):
    t = mult_func(e,cut,s)
    res = np.trapz(t*f,e)  
    return res  

x=np.linspace(0,1,10000)
y=np.exp(x)

m=0.
mm=0.
N=100000  
for i in range(0,N):
    cut=np.random.normal(0.21,0.02)
    stg=np.random.normal(1.1,0.1)
    v=intg(x,y,cut,stg)
    m=m+v
    mm=mm+v*v



print("avg. %10.5E +/- %10.5E fixed %10.5E"%(m/N,np.sqrt((mm/N-(m/N)**2)),intg(x,y,0.21,1.1)))

То, что сделано выше, это просто случайная выборка двух параметров и вычисление среднего значения и дисперсии. Однако я не уверен, насколько этот метод грубой силы является адекватным. Я мог бы использовать закон больших чисел и попытаться оценить, сколько испытаний N необходимы, чтобы получить определенное значение (P=1-1/(N*k**2)) быть рядом k раз стандартные отклонения вокруг истинного среднего.

В принципе то, что я написал, могло сработать. Тем не менее, я предполагаю, что, будучи настолько гибким языком со многими мощными пакетами, python мог бы выполнять эту задачу гораздо более эффективно. Я думал о uncertainties, mcerp а также pymc, Из-за моего ограниченного опыта использования этих пакетов, я не уверен, что делать дальше.

РЕДАКТИРОВАТЬ: Мой оригинальный пример не был настолько информативным, поэтому я решил сделать новый пример, который на самом деле работает, чтобы проиллюстрировать мою идею.

Answer 1

Numpy поддерживает массивы произвольных числовых типов. Однако не все функции поддерживают произвольные числовые типы.

В этом случае оба numpy.exp а также trapz не поддерживаются

Обратите внимание, что модуль uncertanties содержит unumpy пакет. numpy.exp здесь есть замена: uncertainties.unumpy.exp

Мы определяем trapz как ufunc.

Проверьте это здесь!

a=un.ufloat(0.3,0.01)
b=un.ufloat(1.2,0.071)



def sample_func(a: un.UFloat, b: un.UFloat) -> np.ndarray:
    x=np.linspace(0,a,100)
    y = un.unumpy.exp(x)
    return utrapz(y, x)

def utrapz(y: np.ndarray, x: np.ndarray) -> np.ndarray:
    Δx = x[1:]-x[:-1]
    avg_y = (y[1:]+y[:-1])/2
    return (Δx*avg_y)


print(sample_func(a, b))  

ИЗ:

[0.00026601240063021264+/-nan 0.0005935120815465686+/-6.429403852670308e-06
 0.0006973604419223405+/-3.888235103342809e-06 ...,
 0.002095505706899622+/-6.503985178118233e-05
 0.0021019968633076134+/-6.545802781649068e-05
 0.0021084415802710295+/-6.587387316821736e-05]