【问题标题】:Integrate Over Multiple Columns in 1 List to Fill Additional List With Same Number Of Columns整合 1 个列表中的多个列以填充具有相同列数的其他列表
【发布时间】:2020-03-21 12:24:09
【问题描述】:

我打算获取随机变量列表,并通过所述随机变量更改每列中的先前列表。但是,就我的函数而言,每个变量都必须在 Gamma 函数中使用以及积分。

x[t] = c * (1 / (2 ** (v / 2) + test[t - 1]) * (gamma((v / 2) + test[t - 1]))) * integrate.\
                quad(lambda h: np.exp(-h / 2) * h ** ((v / 2) + test[t - 1] - 1), 0, np.inf)

x[ t ] 是一个 np.zeros((x , y)) 列表,而 test[t - 1] 是一个 np.zeros((x - 1, y)) 列表

我已经用适当的随机变量填充了 test[],但是我无法通过这个方程来完成 x 中行 [t] 的列

当我尝试运行我当前的代码时,我收到:

File "C:\Program Files (x86)\Microsoft Visual Studio\Shared\Python37_64\lib\site-packages\scipy\integrate\quadpack.py", line 450, in _quad
    return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
TypeError: only size-1 arrays can be converted to Python scalars

是否有不同的特殊函数允许我使用每列的变量来求解我想要的 x[ t ]?

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import stats
import mpmath as mp
import scipy.integrate as integrate
from scipy.special import gamma

T = 1
beta = 0.5
x0 = 0.05
q = 0
mu = x0 - q
alpha = - (2 - beta) * mu
sigma0 = 0.1
sigma = (2 - beta) * sigma0
b = - ((1 - beta) / (2 * mu) * (sigma0 ** 2))
simulations = 100
M = 50
dt = T / M


def srd_sampled_nxc2():
    x = np.zeros((M + 1, simulations))
    x[0] = x0
    test = np.zeros((M, simulations))
    for t in range(1, M + 1):
        v = 4 * b * alpha / sigma ** 2
        c = (sigma ** 2 * (1 - np.exp(-alpha * dt))) / (4 * alpha)
        nc = np.exp(-alpha * dt) / c * x[t - 1]
        if v > 1:
            x[t] = c * ((np.random.standard_normal(simulations) + nc ** 0.5) ** 2 + mp.nsum(
                lambda i: np.random.standard_normal(simulations) ** 2, [0, v - 1]))
        else:
            max_array = []
            nc_over_2 = [l / 2 for l in nc]
            for p in range(simulations):
                sump = []
                poisson_start = 0
                while poisson_start <= 1:
                    x_i = sum(-np.log(np.random.uniform(0, 1, simulations)) / nc_over_2)
                    sump.append(
                        x_i
                    )
                    poisson_start += x_i
                x_n = max(sump)
                max_array.append(
                    x_n
                )
                sump = []
            test[t - 1] = max_array
            x[t] = c * (1 / (2 ** ((v / 2) + test[t - 1])) * (gamma((v / 2) + test[t - 1]))) * integrate.\
                quad(lambda h: np.exp(-h / 2) * h ** ((v / 2) + test[t - 1] - 1), 0, np.inf)
            max_array = []
    return x

【问题讨论】:

  • 您的问题是test[t - 1] 是一个列表。 quad 需要一个标量。你能建议你真正想做什么吗?你想为每个测试值或其他东西集成吗?
  • 我想通过应用于testt - 1 行每一列的公式来更改x[ ]t 行中的每一列。本质上,使用来自测试每一列的数字来复制x[t]“模拟”次数的公式,从而填充x[]的每一行t行的每一列
  • 我取得了一些进展,但是 while 循环存在问题:对于 t > 2,x_i 变得太小,并且 poisson_start 需要永远达到 1。这是因为nc 变得很大,由于之前的 x 非常大,我猜是因为 x_n 正在增加。如果你愿意,我可以发布代码。
  • 我找到了一种方法来正确存储每一列的积分并将它们称为x[t],但我也刚刚意识到我不应该使用x_n = max(sump);相反,我应该取长度 - 1,所以max_array.append(len(sump) - 1)。我现在似乎面临的问题是x[t] 太小了。我用来采样非中心卡方的公式在大多数情况下为零(鉴于 nc 非常大)我还确保 test[t - 1] 是幂的一部分:((v / 2) + test[t - 1])

标签: python arrays list scipy integration


【解决方案1】:

最终找到了一个易于实施的解决方法:

        else:
            max_array = []
            for p in range(simulations):
                k = nc[t - 1, p]
                lam = k / 2
                poisson_samp = 0
                while poisson_samp <= 1:
                    x_i = -math.log(np.random.uniform(0, 1)) / lam
                    max_array.append(
                        x_i
                    )
                    poisson_samp += x_i
                test[t - 1, p] = len(max_array) - 1
                max_array.clear()
            for f in range(simulations):
                n = test[t - 1, f]
                z = integrate.quad(lambda h: np.exp(-h / 2) * h ** ((v / 2) + n - 1), 0, 1)
                new[t - 1, f] = z[0]
            x[t] = c * (1 / (2 ** ((v / 2) + test[t - 1]) * (gamma((v / 2) + test[t - 1]))) * new[0])

唯一真正的问题是x[t] 的收缩导致除以零——只是一个公式问题。

【讨论】:

    猜你喜欢
    • 2020-12-27
    • 1970-01-01
    • 1970-01-01
    • 2013-05-26
    • 1970-01-01
    • 2016-09-04
    • 1970-01-01
    • 2020-03-04
    • 1970-01-01
    相关资源
    最近更新 更多