<jupyter_start><jupyter_text>Experimenting with Maximum Likelihood Parameter Estimation (MLE) Problem description:A ball is dropped and accelerates without drag. The position of the ball is measured including some measurement noice. We want to use a well known differential equation to determine the unknown coefficient (gravity) using Maximum Likelihood Parameter Estimation. Constraits:We want to solve the diff. eq. numerically by running a "simulation" when we don't know an analytical solution describing the position of the ball. Solutions:Solution1: Measured acceleration is calculated numerically from the measured positions. The Measured acceleration can then be used to fit the parameters of the diff. eq.Somution2: Estimated Position is calculated by solving the differential equation (simulating) for parameter guess.<jupyter_code>import numpy as np, pandas as pd
from matplotlib import pyplot as plt
import seaborn as sns
from scipy.optimize import minimize
import scipy.stats as stats
import statsmodels.api as sm
from statsmodels.base.model import GenericLikelihoodModel
from scipy.integrate import odeint<jupyter_output><empty_output><jupyter_text>Generate data for a ball drop<jupyter_code>def calculate_force(g, m):
    #return g*m - cd*u**2
    return g*m

def calculate_acceleration(g):
    return calculate_force(g=g, m=1)

def ball_drop_air(states, t, g):
    
    # states:
    # [x,u]
        
    x = states[0]
    u = states[1]
    dxdt = u
    
    dudt = calculate_acceleration(g=g)
    
    d_states_dt = np.array([dxdt,dudt])
    
    return d_states_dt
def simulate(states_0, t, g):
    states = odeint(func = ball_drop_air,y0 = states_0,t = t,args=(g,))
    return states
states_0 = np.array([0,0])
t = np.linspace(0,10,100)
g=9.81
m = 1
cd = 0.01
states = simulate(states_0=states_0,t = t, g=g)
x = states[:,0]
dxdt = states[:,1]
fig,axes = plt.subplots(ncols = 2)
fig.set_size_inches(11,5)
ax = axes[0]
ax.plot(t,x,label = 'Ball position')

ax.set_title('Ball drop')
ax.set_xlabel('time [s]')
ax.set_ylabel('x [m]')
a = ax.legend()

ax = axes[1]
ax.plot(t,dxdt,label = 'Ball speed')

ax.set_title('Ball drop')
ax.set_xlabel('time [s]')
ax.set_ylabel('dx/dt [m/s]')
a = ax.legend()<jupyter_output><empty_output><jupyter_text>Add measurement noice<jupyter_code>ϵ = np.random.normal(loc = 0.0, scale = 4.0, size = len(x))
x_measure = x + ϵ
v_measure = np.gradient(x_measure,t)
a_measure = np.gradient(v_measure,t)
df = pd.DataFrame({'t':t, 'x':x_measure, 'v':v_measure, 'a':a_measure})
fig,ax = plt.subplots()
fig.set_size_inches(11,5)
ax.plot(t,x, label = 'Actual ball position')
df.plot(x='t', y='x', label = 'Measured ball position', ax=ax)
ax.set_title('Ball drop')
ax.set_xlabel('time [s]')
ax.set_ylabel('x [m]')
ax.legend()
fig,ax = plt.subplots()
fig.set_size_inches(11,5)
df.plot(x='t', y='a', ax=ax, label='Measured ball acceleration')<jupyter_output><empty_output><jupyter_text>Solution 1Calculating the acceleration numerically from measured position. This gives a very bad prediction of unknown coefficeient (gravity).  Maximum Likelihood   define likelihood function<jupyter_code>def MLERegression(params):
    
    g, std = params[0], params[1] # inputs are guesses at our parameters
    
    t = None
    
    ahat = calculate_acceleration(g=g)
    
    # next, we flip the Bayesian question
    # compute PDF of observed values normally distributed around mean (yhat)
    # with a standard deviation of sd
    negLL = -np.sum( stats.norm.logpdf(a_measure, loc=ahat, scale=std) )
    # return negative LL
    return(negLL)
# let’s start with some random coefficient guesses and optimize
guess = np.array([5,10])
results = minimize(MLERegression, guess, method = 'Nelder-Mead', 
          options={'disp': True})
results
coefficients = results['x'][0:-1]
g_pred = coefficients[0]
g_pred<jupyter_output><empty_output><jupyter_text>Solution 2 Estimated Position is calculated by solving the differential equation (simulating) for parameter guess.<jupyter_code>def MLERegression2(params):
    
    v_0, a_0, g, std = params[0], params[1], params[2], params[3] # inputs are guesses at our parameters
    
    states_0 = [v_0, a_0]
    states = simulate(states_0, t, g)
    xhat = states[:,0]
    # next, we flip the Bayesian question
    # compute PDF of observed values normally distributed around mean (yhat)
    # with a standard deviation of sd
    negLL = -np.sum( stats.norm.logpdf(x_measure, loc=xhat, scale=std) )
    # return negative LL
    return(negLL)
# let’s start with some random coefficient guesses and optimize
guess = np.array([0, 0, 5, 10])
results = minimize(MLERegression2, guess, method = 'Nelder-Mead', 
          options={'disp': True})
results
coefficients = results['x']
states_0_pred = coefficients[0:2]
g_pred = coefficients[2]
states_pred = simulate(states_0=states_0_pred,t = t, g=g_pred)
x_pred = states_pred[:,0]
v_pred = states_pred[:,1]

fig,ax = plt.subplots()
fig.set_size_inches(11,5)
ax.plot(t,x, label = 'Actual ball position')
df.plot(x='t', y='x', label = 'Measured ball position', ax=ax)
ax.plot(t,x_pred, label = 'Predicted ball position')
ax.set_title('Ball drop')
ax.set_xlabel('time [s]')
ax.set_ylabel('x [m]')
ax.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12190003815337658, 'fraction_numerical': 0.016024418161007248, 'mean_word_length': 4.864653243847875, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 21, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '412449', 'license': 'MIT', 'n_tokens_mistral': 1888, 'n_tokens_neox': 1817, 'n_words': 514, 'path': 'notebooks/01.2_maa_MLE_ball_drop.ipynb', 'repo_name': 'martinlarsalbert/rolldecay', 'size': '5242'}