Início / coding / Perguntas / 78874684
Accepted
Ishigami
Asked: 2024-08-15 18:33:15 +0800 CST

Python scipy quad e nqaud dando respostas diferentes

  • 772

Não tenho certeza se isso é mais adequado para troca de pilha matemática ou estouro de pilha, mas como as provas matemáticas parecem boas para mim, suspeito que o problema seja o código e, portanto, vou publicá-lo aqui:

insira a descrição da imagem aqui

A primeira fórmula (fórmula 1) é direta por definição: insira a descrição da imagem aqui

e aqui está meu código para a fórmula 1:

import scipy.integrate as integrate
from scipy.integrate import nquad
from scipy.stats import norm
import math

def normalcdf(x):
    return (1+math.erf(x/math.sqrt(2)))/2

def normalpdf(x): 
    return math.exp(-x*x/2)/math.sqrt(2*math.pi)

def integrand12345(x1,x2,x3,x4,x5,theta1,theta2,theta3,theta4,theta5):
  return normalpdf(x1 - theta1) * normalpdf(x2 - theta2) * normalpdf(x3 - theta3) * normalpdf(x4 - theta4) * normalpdf(x5 - theta5)

def range_x1(theta1,theta2,theta3,theta4,theta5):
  return [-np.inf, np.inf]

def range_x2(x1,theta1,theta2,theta3,theta4,theta5):
  return [x1, np.inf]

def range_x3(x2,x1,theta1,theta2,theta3,theta4,theta5):
  return [x2, np.inf] 

def range_x4(x3,x2,x1,theta1,theta2,theta3,theta4,theta5):
  return [x3, np.inf] 

def range_x5(x4,x3,x2,x1,theta1,theta2,theta3,theta4,theta5):
  return [x4, np.inf]

def pi_12345(theta1,theta2,theta3,theta4,theta5):
  return (nquad(integrand12345, [range_x5, range_x4, range_x3, range_x2, range_x1], args=(theta1,theta2,theta3,theta4,theta5)))[0]

A segunda fórmula (fórmula 2) usa integração dupla: insira a descrição da imagem aqui e aqui está o código da fórmula 2:

def integrandforpi_12(x1, x2, t1, t2, *theta): 
  prod = 1
  for ti in theta:
    prod = prod * (1 - normalcdf(x2 - ti))
  return prod * normalpdf(x1 - t1) * normalpdf(x2 - t2)

def range_x1(t1, t2, *theta):
  return [-np.inf, np.inf]

def range_x2(x1, t1, t2, *theta):
  return [x1, np.inf]

def pi_12(t1, t2, *theta):
  return (nquad(integrandforpi_12, [range_x2, range_x1], args=(t1, t2, *theta)))[0]

Minha terceira fórmula (fórmula 3) é baseada no teorema de Bayes: insira a descrição da imagem aqui

e meu código para a fórmula 3 é:

def integrandforpi_i(xi, ti, *theta):
  prod = 1
  for t in theta:
    prod = prod * (1 - normalcdf(xi - t))
  return prod * normalpdf(xi - ti)

def pi_i(ti, *theta):
  return integrate.quad(integrandforpi_i, -np.inf, np.inf, args=(ti, *theta))[0]

Então pi_i calcula a probabilidade de que X_i seja o menor entre os theta_i.

No entanto, quando executo meu código usando as 3 fórmulas diferentes, todas elas dão respostas diferentes e não tenho ideia do porquê. Não tenho certeza se há uma falha na minha derivação da fórmula ou se há um erro no meu código. Qualquer ajuda seria apreciada.

Aqui está um exemplo:

t1,t2,t3,t4,t5 = 0.83720022,0.61704171,1.21121701,0,1.52334078

p12345 = pi_12345(t1,t2,t3,t4,t5)
p12354 = pi_12345(t1,t2,t3,t5,t4)
p12435 = pi_12345(t1,t2,t4,t3,t5)
p12453 = pi_12345(t1,t2,t4,t5,t3)
p12534 = pi_12345(t1,t2,t5,t3,t4)
p12543 = pi_12345(t1,t2,t5,t4,t3)

print('p12345=',p12345)
print('p12354=',p12354)
print('p12435=',p12435)
print('p12453=',p12453)
print('p12534=',p12534)
print('p12543=',p12543)

print('formula 1 gives', p12345+p12354+p12435+p12453+p12534+p12543)

print('formula 2 gives', pi_12(t1,t2,t3,t4,t5))

print('formula 3 gives', pi_i(t2,t3,t4,t5) * pi_i(t1,t2,t3,t4,t5))

e a saída é

p12345= 0.0027679276698449086
p12354= 0.008209750140618218
p12435= 0.0016182955786153714
p12453= 0.001921206801273682
p12534= 0.009675713474375739
p12543= 0.003904872716765966
formula 1 gives 0.028097766381493885
formula 2 gives 0.21897431741874426
formula 3 gives 0.0418669679120933

Observação: A Fórmula 1 é extremamente lenta, leva cerca de 3 horas para rodar no meu pobre e velho laptop. As Fórmulas 2 e 3 são instantâneas.

python

2 respostas

  1. Best Answer

    Como a nquadlista de intervalos está invertida, a Xlista de argumentos do integrando também deve ser invertida:

    def integrand12345(x5,x4,x3,x2,x1,theta1,theta2,theta3,theta4,theta5):

    e

    def integrandforpi_12(x2, x1, t1, t2, *theta):

    Então meus resultados se tornam:

    formula 1 gives 0.04444581969881939
formula 2 gives 0.04444465903549115
formula 3 gives 0.0418669679120933

    Também diminuí significativamente o tempo de cálculo da fórmula 1 passando o optsargumento para nquad:opts={'epsabs':0.001}

    • 3

  2. Isso não responde totalmente à pergunta, mas pode ser útil saber se alguma delas está correta. Para isso, temos simulação:

    import numpy as np
from scipy import stats
rng = np.random.default_rng(43258234589253)

ti = np.asarray([0.83720022, 0.61704171, 1.21121701, 0, 1.52334078])

n_trials = 1000000
n_obs = len(ti)
desired_ranks = np.arange(1, n_obs+1)

Xi = rng.standard_normal(size=(n_trials, n_obs)) + ti
ranks = stats.rankdata(Xi, axis=1)
in_order = np.all(ranks[:, :2] == desired_ranks[:2], axis=1)
p = np.mean(in_order)
print(p)  # 0.044489

    Se este código estiver correto (e corrija-me se entendi mal e corrigirei isso), então nenhuma das integrais está correta.

    Na verdade, também simulei P(x1 é o menor de todos xi) * P(x2 é o menor de x2, x3, x4, x5). Dá o mesmo resultado que a integral 3, mas não o mesmo que a simulação completa ou integrais 1/2 (após correção). Na sua derivação, você está tratando P(x2 é o segundo menor | x1 é o menor) como equivalente a P(x2 é o menor de x2, x3, x4, x5). Essa parte pode não estar certa, porque ignora o fato de que x2 > x1.

    • 1

relate perguntas