0
ルックアヘッド推定器を使用して区分的確率密度関数の定常状態分布を推定しようとしています。しかし、エラーを取得するQ1:IndexError:サイズ1の軸1に対してインデックス1の範囲外です。
z[condition1] = (1/(sigma*np.sqrt(2*np.pi)))*np.e**(-0.5*((x[condition1]-y[condition1]+Q1-u)/sigma)**2.
IndexError: index 1 is out of bounds for axis 1 with size 1.
誰もが最初の条件が動作していない理由を教えてもらえますか?私はどうだろう何
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt
import sympy as sp
from sympy import *
import numpy as np
from scipy.stats import lognorm, beta
from quantecon import LAE
from sympy import symbols
q= symbols('q')
## == Define parameters == #
a_sigma = 0.4
psi_0 = beta(5, 5, scale=0.5) # Initial distribution
phi = lognorm(a_sigma)
def p(x,y):
u=80
sigma=30
b=0.2
Q=80
Q1=Q*(1-b)
Q2=Q*(1+b)
z = np.zeros_like(x, dtype=float)
# Condition 1 indexes all elements where subformula 1 is valid
condition1 = np.logical_and(y>0.0, x >=y-Q1)
z[condition1] = (1/(sigma*np.sqrt(2*np.pi)))*np.e**(-0.5*((x[condition1]-y[condition1]+Q1-u)/sigma)**2)
condition2 = np.logical_and(y<0.0, x >=y-Q2)
z[condition2] = (1/(sigma*np.sqrt(2*np.pi)))*np.e**(-0.5*((x[condition2]-y[condition2]+Q2-u)/sigma)**2)
condition3 = np.logical_and(y==0.0, x >=-Q1)
#print(-0.5*((k_prime[condition3] + q - u))**2)
j=-0.5*((x[condition3] + q - u))**2
K=[]
for elem in j:
print(elem)
K.append(1/(sigma*sqrt(2*pi))*sp.integrate(sp.exp(elem),(q,Q1,Q2)))
#z[condition3] = K
return z
n = 10000 # Number of observations at each date t
T = 30 # Compute density of k_t at 1,...,T+1
# == Generate matrix s.t. t-th column is n observations of k_t == #
k = np.empty((n, T))
A = phi.rvs((n, T))
k[:, 0] = psi_0.rvs(n) # Draw first column from initial distribution
for t in range(T-1):
k[:, t+1] = k[:, t]+ A[:, t]
# == Generate T instances of LAE using this data, one for each date t == #
laes = [LAE(p, k[:, t]) for t in range(T)]
# == Plot == #
fig, ax = plt.subplots()
ygrid = np.linspace(0.01, 4.0, 200)
greys = [str(g) for g in np.linspace(0.0, 0.8, T)]
greys.reverse()
for psi, g in zip(laes, greys):
ax.plot(ygrid, psi(ygrid), color=g, lw=2, alpha=0.6)
ax.set_xlabel('capital')
title = r'Density of $k_1$ (lighter) to $k_T$ (darker) for $T={}$'
ax.set_title(title.format(T))
plt.show()