dimanche 24 novembre 2019

Poisson (npr) Size Alteration Returns ValueError (wrt arbitrary paths and array creation)

If I sample a non-central chi-square distribution using a Poisson distribution, I am unable to alter the size and can only input the mean, "nc / 2" (I must set size = 1 or it also returns the same error):

n = np.random.poisson(nc / 2, 1)  # generates a random variable from the poisson distribution with
# mean: non-centrality parameter / 2
x[t] = c * mp.nsum(lambda i: np.random.standard_normal() ** 2, [0, v + 2 * n])

If I attempt to increase the size to the number of simulations being run

n = np.random.poisson(nc / 2, simulations)

where simulations = 10000, I receive:

"ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()"

Running the code with 1 simulation produces one desired result, and every run produces another random path.

Graph created under 10,000 simulations with size = one

However, it is a necessity to have the graph composed of paths determined by each iteration of the simulation. Under a different condition, the non-central chi-square distribution is determined by the code:

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]))

which does produce the desired result

Graph produced by the line of code above

How can I obtain a different path for x[t] despite not being able to change the size of the Poisson distribution (i.e. not have the same path for each of the 10,000 simulations)

If required:

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

T = 1
beta = 1.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) * sigma ** 2
simulations = 10000
M = 50
dt = T / M


def srd_sampled_nxc2():
    x = np.zeros((M + 1, simulations))
    x[0] = x0
    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]  # The non-centrality parameter lambda
        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:
            n = np.random.poisson(nc / 2, 1)
            x[t] = c * mp.nsum(lambda i: np.random.standard_normal() ** 2, [0, v + 2 * n])
    return x


x1 = srd_sampled_nxc2()


plt.figure(figsize=(10, 6))
plt.plot(x1[:, :10], lw=1)
plt.xlabel('time')
plt.ylabel('index')
plt.show()



Aucun commentaire:

Enregistrer un commentaire