<jupyter_start><jupyter_text>Poisson $$f\left(\left.y\right|x_{i}\right)=\frac{\exp\left(-\mu\left(x_{i}\right)\right)\mu\left(x_{i}\right)^{y}}{y!}$$$$\mu\left(X_{i}\right)=\exp\left(X_{i}\theta\right)$$<jupyter_code>import numpy as np
import matplotlib.pylab as plt
import seaborn as sns

np.set_printoptions(precision=4, suppress=True)
sns.set_context('notebook')

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Generate data<jupyter_code># True parameter
theta = .5
# Sample size
n = int(1e2)
# Independent variable, N(0,1)
X = np.random.normal(0, 1, n)
# Sort data for nice plots
X = np.sort(X)

mu = np.exp(X * theta)
# Error term, N(0,1)
Y = np.random.poisson(mu, n)<jupyter_output><empty_output><jupyter_text>Plot the data and the model<jupyter_code>plt.figure(figsize = (8, 8))

plt.scatter(X, Y, label='Observed data')
plt.ylabel(r'$Y$')
plt.xlabel(r'$X$')

plt.show()<jupyter_output><empty_output><jupyter_text>Maximize log-likelihood $$l\left(y|x,\theta\right)=\sum_{i=1}^{n}\log\frac{\exp\left(-\mu\left(x_{i}\right)\right)\mu\left(x_{i}\right)^{y_{i}}}{y_{i}!}$$<jupyter_code>import scipy.optimize as opt
from scipy.stats import poisson

# Define objective function
def f(theta, X, Y):
    Q = - np.sum(np.log(1e-3 + poisson.pmf(Y, np.exp(X * theta))))
    return Q

# Run optimization routine
theta_hat = opt.fmin_bfgs(f, 0., args=(X, Y))

print(theta_hat)<jupyter_output>Optimization terminated successfully.
         Current function value: 130.738440
         Iterations: 5
         Function evaluations: 21
         Gradient evaluations: 7
[ 0.4621]<jupyter_text>Plot objective function, true parameter, and the estimate<jupyter_code># Generate data for objective function plot
th = np.linspace(-3., 3., 1e2)
Q = [f(z, X, Y) for z in th]

# Plot the data
plt.figure(figsize=(8, 4))
plt.plot(th, Q, label='Q')
plt.xlabel(r'$\theta$')
plt.axvline(x=theta_hat, c='red', label='Estimated')
plt.axvline(x=theta, c='black', label='True')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Solve first order conditions<jupyter_code>from scipy.optimize import fsolve

# Define the first order condition
def df(theta, X, Y):
    return - np.sum(X * (Y - np.exp(X * theta)))

# Solve FOC
theta_hat = fsolve(df, 0., args=(X, Y))

print(theta_hat)<jupyter_output>[ 0.4619]<jupyter_text>Plot first order condition<jupyter_code># Generate data for the plot
th = np.linspace(-3., 3., 1e2)
Q = np.array([df(z, X, Y) for z in th])

# Plot the data
plt.figure(figsize=(8, 4))
plt.plot(th, Q, label='Q')
plt.xlabel(r'$\beta$')
plt.axvline(x=theta_hat, c='red', label='Estimated')
plt.axvline(x=theta, c='black', label='True')
plt.axhline(y=0, c='green')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Plot original data and fitted mean<jupyter_code>plt.figure(figsize=(8, 8))

plt.scatter(X, Y, label='Observed data')
plt.plot(X, np.exp(X * theta_hat), label='Fitted mean')
plt.ylabel(r'$Y$')
plt.xlabel(r'$X$')
plt.legend()

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18263175339741464, 'fraction_numerical': 0.018230029830957905, 'mean_word_length': 5.766816143497758, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '631598', 'license': 'MIT', 'n_tokens_mistral': 1292, 'n_tokens_neox': 1211, 'n_words': 240, 'path': 'notebooks/mle_poisson.ipynb', 'repo_name': 'khrapovs/metrix', 'size': '3017'}