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