In [1]:
# Import relevant libraries
from __future__ import division, print_function
import numpy as np
import scipy.stats as stats
from scipy.stats.mstats import mquantiles
from scipy.optimize import minimize
The following is a very compact description of the method. See the authors' paper for a much more in-depth discussion.
The key idea is that we define two objects:
This skips a lot of details. For example:
That said, I will use the minimal notation needed to describe the process, and use an example to explore and illustate the details of the process.
Assume we have two things:
We will explore this in more detail in the example below.
To fit the theoretical model to the empirical data, choose the parameter $\theta$ to minimize the following quadratic objective function:
$$ \hat{\theta} = \underset{\theta \in \Theta}{\textrm{argmin}} \; \left(\hat{\mathbf{\Phi}} - \mathbf{\Phi}_{\theta}\right)^{\textrm{T}} \mathbf{W}_{\theta} \left(\hat{\mathbf{\Phi}} - \mathbf{\Phi}_{\theta}\right). $$Here $\mathbf{W}_{\theta}$ is a symmetric positive definite weighting matrix.
In addition, the bolded $\Phi$ values are defined as the "stacked" vectors, for each of the $J$ random variables:
$$ \mathbf{\Phi} = \left(\Phi_{1}^{\textrm{T}}, \Phi_{2}^{\textrm{T}}, ..., \Phi_{j}^{\textrm{T}}, ..., \Phi_{J}^{\textrm{T}} \right)^{\textrm{T}} $$for both $\hat{\Phi}$ and $\Phi_{\theta}$.
We'll use the example of the $\alpha$-stable distribution to demonstrate the estimation method.
As discussed in the notebook "An Aside on Alpha Stable Distributions," the $\alpha$-stable distribution is a four-parameter distribution denoted $S(\alpha, \beta, \mu, \sigma)$ where:
There are three special cases of the family of stable distributions:
Importantly, the $\alpha$ parameter governs whether moments of the distribution exist. For $X \sim S(\alpha, \beta, \mu, \sigma)$:
$$\mathbb{E} \left[ X^{p}\right] \lt \infty \; \forall p \lt \alpha .$$We use a quantiles-based estimator defined by McCulloch (1986).
Let's define the two functions of parameters:
Dominicy and Verdas (2013) follow McCulloch (1986) and standardize the empirical data which comprises $\hat{\Phi}$, and thus additionally standardize the theoretical measurements.
The theoretical function of parameters $\Phi_{\theta}$ is obtained by drawing a simulated sample of data from the distribution implied by $\theta$, which we will denote as $\overset{\sim}{\Phi}_{\theta}$.
Finally, to represent $\Phi_{\theta}$ we draw $R$ simulation paths and find the average value:
$$ \overset{\sim}{\Phi}_{\theta}^{R} = \frac{1}{R}\sum_{r=1}^{R} \overset{\sim}{\Phi}_{\theta}^{r}. $$We find the $\hat{\theta}$ by minimizing the quadratic objective:
$$ \hat{\theta} = \underset{\theta \in \Theta}{\textrm{argmin}} \; \left(\hat{\mathbf{\Phi}} - \mathbf{\overset{\sim}{\Phi}}_{\theta}^{R} \right)^{\textrm{T}} \mathbf{W}_{\theta} \left(\hat{\mathbf{\Phi}} - \mathbf{\overset{\sim}{\Phi}}_{\theta}^{R}\right). $$There are a few more details which we cover as we proceed.
In [2]:
# Let's define some functions:
def h(q):
# Assume qhat = qhat_{.05}, qhat_{.25}, qhat_{.50}, qhat_{.75}, qhat_{.95}
return np.array([ q[-1] - q[0],
q[-1] - 2 * q[2] + q[0],
q[-2] - q[1],
q[2] ])
def g(q):
# Assume qhat = qhat_star_{.05}, qhat_star_{.25}, qhat_star_{.50}, qhat_star_{.75}, qhat_star_{.95}
return np.array([ 1.0 / (q[-2] - q[1]),
1.0 / (q[-1] - q[0]),
1.0,
1.0 ])
# Test this:
def Phi_hat_alt(q):
return h(q) * g(q)
def Phi_hat(q):
'''
q is a vector of quantiles in ascending order. The function will
error-check for this.
'''
# Check if q is in order:
#assert np.all(np.diff(q) >= 0), "Quantiles q are not in ascending order: q="+str(q)
return np.array([ (q[-1] - q[0]) / (q[-2] - q[1]),
(q[-1] - 2 * q[2] + q[0]) / (q[-1] - q[0]),
q[-2] - q[1],
q[2] ])
def Phi_theta_plain(q, mu, sigma):
'''
q is a vector of quantiles in ascending order. The function will
error-check for this.
mu and sigma are float values for the mu, sigma values in the function.
'''
# Assert that q is in ascending order!
#assert np.all(np.diff(q) >= 0), "Quantiles q are not in ascending order: q="+str(q)
return np.array( [ (q[-1] - q[0]) / (q[-2] - q[1]),
(q[-1] - 2 * q[2] + q[0]) / (q[-1] - q[0]),
(q[-2] - q[1]) * sigma,
q[2] * sigma + mu] )
def Phi_theta(q, theta):
'''
q is a vector of quantiles in ascending order. The function will
error-check for this.
mu and sigma are float values for the mu, sigma values in the function.
'''
# Assert that q is in ascending order!
#assert np.all(np.diff(q) >= 0), "Quantiles q are not in ascending order: q="+str(q)
# Recall:
#theta = {'alpha':theta_vec[0],
# 'beta':theta_vec[1],
# 'mu':theta_vec[2],
# 'sigma':theta_vec[3]}
if theta[0] != 1.0:
zeta = theta[2] + theta[1]*theta[3]*np.tan(np.pi * theta[0] / 2.0)
# zeta= mu + beta * sigma * tan( pi * alpha / 2)
else:
zeta = theta[2] # mu
return np.array( [ (q[-1] - q[0]) / (q[-2] - q[1]),
(q[-1] - 2 * q[2] + q[0]) / (q[-1] - q[0]),
(q[-2] - q[1]) * theta[3],
q[2] * theta[3] + zeta] )
def Phi_theta_dict(q, theta):
'''
q is a vector of quantiles in ascending order. The function will
error-check for this.
mu and sigma are float values for the mu, sigma values in the function.
'''
# Assert that q is in ascending order!
#assert np.all(np.diff(q) >= 0), "Quantiles q are not in ascending order: q="+str(q)
# Recall:
#theta = {'alpha':theta_vec[0],
# 'beta':theta_vec[1],
# 'mu':theta_vec[2],
# 'sigma':theta_vec[3]}
if theta['alpha'] != 1.0:
zeta = theta['mu'] + theta['beta']*theta['sigma']*np.tan(np.pi * theta['alpha'] / 2.0)
# zeta= mu + beta * sigma * tan( pi * alpha / 2)
else:
zeta = theta['mu'] # mu
return np.array( [ (q[-1] - q[0]) / (q[-2] - q[1]),
(q[-1] - 2 * q[2] + q[0]) / (q[-1] - q[0]),
(q[-2] - q[1]) * theta['sigma'],
q[2] * theta['sigma'] + zeta] )
# Generate data
def generate_data(theta, N, R, rng):
'''
The parameter theta contains the alpha-stable distribution parameters.
theta must be a dict with keys: alpha, beta, mu, sigma.
N is length of vector of IID draws, R is number of vector.
rng is a numpy RandomState instance.
seed is for a random number generator.'''
#assert theta['alpha'] > 0 and theta['alpha'] <= 2, "alpha not in (0,1]: alpha = "+str(theta['alpha'])
#assert theta['beta'] >= -1 and theta['beta'] <= 1, "beta not in [-1,1]: beta = "+str(theta['beta'])
#assert theta['sigma'] >= 0 , "sigma not >= 0: sigma = "+str(theta['sigma'])
# Generate the data
return stats.levy_stable.rvs(alpha=theta['alpha'], beta=theta['beta'],
loc=theta['mu'], scale=theta['sigma'],
size=(N, R), random_state=rng)
def generate_data_vec(theta, N, R, rng, testing=False):
'''
Same as "generate_data" but with theta as an ordered numpy vector.
The parameter theta contains the alpha-stable distribution parameters.
theta values must be in the order: alpha, beta, mu, sigma.
N is length of vector of IID draws, R is number of vector.
rng is a numpy RandomState instance.
seed is for a random number generator.'''
#if testing:
# assert theta[0] > 0 and theta[0] <= 2, "alpha not in (0,1]: alpha = "+str(theta[0])
# assert theta[1] >= -1 and theta[1] <= 1, "beta not in [-1,1]: beta = "+str(theta[1])
# assert theta[3] >= 0 , "sigma not >= 0: sigma = "+str(theta[3])
# Generate the data
return stats.levy_stable.rvs(alpha=theta[0], beta=theta[1],
loc=theta[2], scale=theta[3],
size=(N, R), random_state=rng)
def find_emperical_quantiles(data, q):
'''
data is a 1D or 2D vector of data,
q are the quantiles.
'''
# Draw quantiles from the data. NOTE: Assume that data is (N,R)
# shape. Note that returned values will have shape: (len(q), R)
# NOTE2: this is only for exposition; just directly use mquantiles.
return mquantiles(a=data, prob=q, alphap=0.4, betap=0.4, axis=0)
def find_theoretical_quantiles_squiggle_R(theta, N, R, q, rng):
'''
Construct the thereotical quantiles by simulation. Given
theta, N, R, q, and rng, return a vector of the five quantiles
associated with this distribution.
Parameters:
-----------
theta: dict of stable distribution parameters. Must be a dict
with keys: alpha, beta, mu, sigma
N: int, number of observations per simulation
R: int, number of simulations to run
q: vector of floats; the quantiles values. Assumed in ascending order.
rng: a NumPy RandomState object. Allows reproducibility.
Returns:
-----------
q_hat: vector of floats, the average of quantilels over R simulations
'''
# Assert that q is in ascending order!
#assert np.all(np.diff(q) > 0), "Quantiles q are not in ascending order: q="+str(q)
# Generate data:
# NOTE: data is (N,R) shape:
data = generate_data(theta, N, R, rng)
# Find the quantiles:
# Draw quantiles from the data. NOTE: Assume that data is (N,R)
# shape. Note that returned values will have shape: (len(q), R)
quantiles_R = mquantiles(a=data, prob=q, alphap=0.4, betap=0.4, axis=0)
# Average over each quantile; the resulting vector will be
#in ascending order:
return np.apply_along_axis(func1d=np.mean, axis=1, arr=quantiles_R) #theory_quantiles
def find_theoretical_quantiles_squiggle_R_vec(theta_vec, N, R, q, rng, testing=False):
'''
Construct the thereotical quantiles by simulation. Given
theta, N, R, q, and rng, return a vector of the five quantiles
associated with this distribution.
Parameters:
-----------
theta: vector of stable distribution parameters. Must be a numpy array
with floats in this order: alpha, beta, mu, sigma
N: int, number of observations per simulation
R: int, number of simulations to run
q: vector of floats; the quantiles values. Assumed in ascending order.
rng: a NumPy RandomState object. Allows reproducibility.
Returns:
-----------
q_hat: vector of floats, the average of quantilels over R simulations
'''
# Assert that q is in ascending order!
#if testing:
# assert np.all(np.diff(q) > 0), "Quantiles q are not in ascending order: q="+str(q)
# Generate data:
# NOTE: data is (N,R) shape:
data = generate_data_vec(theta_vec, N, R, rng)
# Find the quantiles:
# Draw quantiles from the data. NOTE: Assume that data is (N,R)
# shape. Note that returned values will have shape: (len(q), R)
quantiles_R = mquantiles(a=data, prob=q, alphap=0.4, betap=0.4, axis=0)
# Average over each quantile; the resulting vector will be
#in ascending order:
return np.apply_along_axis(func1d=np.mean, axis=1, arr=quantiles_R) #theory_quantiles =
def sparsity_function():
pass
def G_hat():
# Diagonal matrix
# Diagonal elements are g( q_hat^{*}_{j} )
# Full thing:
#
# g(q_hat^{*}) = ( g(q_hat^{*}_{1})^T, g(q_hat^{*}_2)^T, ... , g(q_hat^{*}_{J})^T )^T
#
pass
In [3]:
# Next steps:
# Generate a sample of data from a "mystery dist"
# Choose starting theta
# run estimation *step 1 only*.
q = [0.05, 0.25, 0.5, 0.75, 0.95]
mystery_theta = {'alpha':1.5,
'beta':-0.5,
'mu':0.5,
'sigma':1.2}
# Cauchy: alpha=1, beta=0, loc=mu=0.0, scale=sigma/2.0 = 1.0/2.0
# Using raw data:
# array([ 1.33024131, -0.57463142, -0.16851961, 0.96667289])
# Using standardized sample, with sigma=std(data):
#array([ 1.33023827, -0.57449001, 0.06295755, 0.21148216])
'''
Optimization terminated successfully. (Exit mode 0)
Current function value: 3.93901180299e-08
Iterations: 15
Function evaluations: 95
Gradient evaluations: 15
Took:7.11473720074 min.
Out[7]:
fun: 3.9390118029863406e-08
jac: array([-0.00065797, 0.00016368, 0.00025707, 0.00022284, 0. ])
message: 'Optimization terminated successfully.'
nfev: 95
nit: 15
njev: 15
status: 0
success: True
x: array([ 1.33023827, -0.57449001, 0.06295755, 0.21148216])
'''
# SO NOTE: THIS APPEARS WRONG! using std.
# Next:
# - try with 'use sigma/mu to stdize sample'
# - update W_theta
# - ...email
# RNG setup:
seed0 = 567891234
rng = np.random.RandomState(seed0)
# Draw sample of mystery data:
Msample = 200
empirical_data = generate_data(theta=mystery_theta, N=Msample, R=1, rng=rng)
#N=10000
#R=200
N=7500
R=150
# Standardize the sample?
# Need to ask authors about this!
z_empirical_data = (empirical_data - np.mean(empirical_data) ) / np.std(empirical_data)
# NOTE: confusion over what "standardize the sample" means. Here are there options:
# 1. don't standardize
# 2. standardize using mean, std.
# 3. for each theta_check, after generatefs all the samples,
# figure out quantiles, and then a,b,mu, sig, use *that* mu andd sig to normalize.
# Still not clear why would do this (3) one. But whatevs.
# SO, start and try out the (2) version...yes?
# Process:
# Form the quantiles over the data
# Form W0
# Choose a theta0
# Generate theta_squiggle_R from theta0
# Find the quandratic value
# Iterate
# Get data quantiles:
W0 = np.eye(4)
#empirical_q = mquantiles(a=z_empirical_data, prob=q, alphap=0.4, betap=0.4, axis=0).compressed() # Remove masking
empirical_q = mquantiles(a=empirical_data, prob=q, alphap=0.4, betap=0.4, axis=0).compressed() # Remove masking
theta_hat = Phi_hat(empirical_q)
theta0 = {'alpha':1.2,
'beta':-0.25,
'mu':0.5,
'sigma':0.5} # Levy: alpha=0.5, beta=1.0, loc=mu, scale=sigma/2.0
# Note: *internal* to the function, unpack these!
theta_vec = np.array([theta0['alpha'], theta0['beta'], theta0['mu'], theta0['sigma'], ])
def quadratic_objective_old_dict_style(theta, W, theta_hat, seed, N, R, qvals):
#unpack:
#theta = {'alpha':theta_vec[0],
# 'beta':theta_vec[1],
# 'mu':theta_vec[2],
# 'sigma':theta_vec[3]}
# Generate theta_squiggle_R from theta0:
rng = np.random.RandomState(seed)
theta_squiggle_R_q = find_theoretical_quantiles_squiggle_R(theta=theta, N=N, R=R, q=qvals, rng=rng)
theta_squiggle_R = Phi_theta_dict(q=theta_squiggle_R_q, theta=theta) # mu=theta['mu'], sigma=theta['sigma'])
theta_minus_theta = theta_hat - theta_squiggle_R
return np.dot(theta_minus_theta, W0).dot(theta_minus_theta)
def quadratic_objective(theta_vec, W, theta_hat, seed, N, R, qvals):
# Generate theta_squiggle_R from theta0:
# NOTE: this takes in the "theta_hat" already calculated. Need to think
# about whether to *always* recalcuate the the theta_hat based on the
# empirical data, **standardized by the theoretical vavlues**.
# I *presume* that will be more expensive.
rng = np.random.RandomState(seed)
theta_squiggle_R_q = find_theoretical_quantiles_squiggle_R_vec(theta_vec=theta_vec, N=N, R=R, q=qvals, rng=rng)
theta_squiggle_R = Phi_theta(q=theta_squiggle_R_q, theta=theta_vec) # mu=theta_vec[2], sigma=theta_vec[3])
theta_minus_theta = theta_hat - theta_squiggle_R
return np.dot(theta_minus_theta, W0).dot(theta_minus_theta)
def quadratic_objective_vec_recalculate_theta_hat(theta_vec, W, seed, N, R, qvals, verbose=False, the_data=empirical_data):
# Generate theta_squiggle_R from theta0:
# This version takes in the empirical data and *always* recalcuates
# the the theta_hat based on the
# empirical data, **standardized by the theoretical vavlues**.
# I *presume* that will be more expensive.
# Note that I am *binding* the empirical data to the function. Will see how that goes...
# Recall:
# theta = {'alpha':theta_vec[0],
# 'beta':theta_vec[1],
# 'mu':theta_vec[2],
# 'sigma':theta_vec[3]}
if verbose:
print("mu="+str(theta_vec[2]))
print("sigma="+str(theta_vec[3]))
# First find theta-hat:
standardized_data = (the_data - theta_vec[2]) / theta_vec[3]
# Now find the quantiles:
empirical_q = mquantiles(a=standardized_data, prob=q, alphap=0.4, betap=0.4, axis=0).compressed() # Remove masking
# *now* find the theta-hat:
theta_hat = Phi_hat(empirical_q)
# Now rest of the run:
rng = np.random.RandomState(seed)
theta_squiggle_R_q = find_theoretical_quantiles_squiggle_R_vec(theta_vec=theta_vec, N=N, R=R, q=qvals, rng=rng)
theta_squiggle_R = Phi_theta(q=theta_squiggle_R_q, theta=theta_vec) #mu=theta_vec[2], sigma=theta_vec[3])
theta_minus_theta = theta_hat - theta_squiggle_R
return np.dot(theta_minus_theta, W0).dot(theta_minus_theta)
seed1=1236789
quadratic_objective(theta_vec, W0, theta_hat, seed1, N, R, q)
# Let's examine:
#theta_vec_true = np.array([theta0['alpha'], theta0['beta'], theta0['mu'], theta0['sigma'], ])
print(quadratic_objective(theta_vec, W0, theta_hat, seed1, N, R, q))
print(quadratic_objective_old_dict_style(theta0, W0, theta_hat, seed1, N, R, q))
print(quadratic_objective_vec_recalculate_theta_hat(theta_vec, W0, seed1, N, R, q))
In [4]:
theta_hat
def G_theta(q):
return np.diag(g(q))
def Omega_theta(q):
return np.diag(h(q))
def Wstar_theta(q):
pass
g(empirical_q)
h(empirical_q)
#np.diag(g(empirical_q))
#np.diag(h(empirical_q))
#W0
#G_theta(empirical_q)
Omega_theta(empirical_q)
Out[4]:
Denote by $\hat{G}_j$ a M x M diagonal matrix with diagonal elements $g(\hat{q}_{j}^{*})$.
We gather all these vectors into $$g( \hat{q}^{*}) = \left( g(\hat{q}_{1}^{*})^{T}, ..., g(\hat{q}_{J}^{*})^{T} \right)^{T} $$
with the corresponding block diagonal JM x JM matrix
$$\hat{G} = diag(diag(g(\hat{q}_{1}^{*}), ..., diag(g(\hat{q}_{J}^{*}))).$$Similarly, let $G_{\theta}$ be a JM × JM diagonal matrix composed of J diagonal blocks, each of size M :
$$G_{\theta} = diag(diag( g(q_{\theta,1}^{*} )), . . . , diag( g(q_{\theta,M}^{*} ))).$$A result: $\hat{G} \rightarrow G_{\theta}$ as __.
Major Q: what is the convergence iteration here? How do we go from $\hat{q}_{k}^{*} \rightarrow \hat{q}_{k+1}^{*}?$
Are we doing the following?
Updated the Phi_theta function to include the zeta value.
BUT getting the wrong sign on mu. Need to think about -- see if can swtichc back and get "old"results.
Also need to think more about iteration to $W_{\theta}$. Need to be careful about what exactly the iteration steps are.
Need also to think aa
Also need to email Qs.
In [ ]:
# NEW TRIAL: just trying out a "simple" two-step estimator. Key note: see Sheppard 6.4.1
# Recall:
def quadratic_objective(theta_vec, W, theta_hat, seed, N, R, qvals):
# Generate theta_squiggle_R from theta0:
# NOTE: this takes in the "theta_hat" already calculated. Need to think
# about whether to *always* recalcuate the the theta_hat based on the
# empirical data, **standardized by the theoretical vavlues**.
# I *presume* that will be more expensive.
rng = np.random.RandomState(seed)
theta_squiggle_R_q = find_theoretical_quantiles_squiggle_R_vec(theta_vec=theta_vec, N=N, R=R, q=qvals, rng=rng)
theta_squiggle_R = Phi_theta(q=theta_squiggle_R_q, theta=theta_vec) # mu=theta_vec[2], sigma=theta_vec[3])
theta_minus_theta = theta_hat - theta_squiggle_R
return np.dot(theta_minus_theta, W0).dot(theta_minus_theta)
from time import time
t0=time()
near_zero = 1.0e-8
res = minimize(fun=quadratic_objective_vec_recalculate_theta_hat, x0=theta_vec, args=(W0, seed1, N, R, q), method='SLSQP',
jac=None, bounds=[(near_zero, 2.0), (-1.0+near_zero, 1.0-near_zero), (None, None), (near_zero, None)], options={"disp":True} )
t1=time()
print("Took:"+str( (t1-t0)/60. )+ " min.")
try:
print(res.x)
print("printed res.x -- res is an object, x is attribute")
except:
print(res['x'])
print("printed res['x'] -- res is a dict, x is key")
# Next step:
# Need to implement the 2-step estimator, with the recommended update of the W matrix.
In [ ]:
# Usign old version:
from time import time
t0=time()
#quadratic_objective_vec
near_zero = 1.0e-8
res = minimize(fun=quadratic_objective, x0=theta_vec, args=(W0, theta_hat, seed1, N, R, q), method='SLSQP',
jac=None, bounds=[(near_zero, 2.0), (-1.0+near_zero, 1.0-near_zero), (None, None), (near_zero, None)], options={"disp":True} )
#res = minimize(fun=quadratic_objective, x0=theta_vec, args=(W0, theta_hat, seed1, N, R, q), method='SLSQP',
# jac=None, bounds=[(0.51, 2.0), (-0.99999999, 0.9999999), (None, None), (0.0, None)], options={"disp":True} )
t1=time()
print("Took:"+str( (t1-t0)/60. )+ " min.")
res
In [ ]:
# USING THE NEW -ALWAYS-RECALCUALTE VERSION:
from time import time
t0=time()
near_zero = 1.0e-12
res = minimize(fun=quadratic_objective_vec_recalculate_theta_hat, x0=theta_vec, args=(W0, seed1, N, R, q), method='SLSQP',
jac=None, bounds=[(near_zero, 2.0), (-1.0+near_zero, 1.0-near_zero), (None, None), (near_zero, None)], options={"disp":True} )
t1=time()
print("Took:"+str( (t1-t0)/60. )+ " min.")
res
In [ ]:
# REsults with standardized-my-theory-values version:
'''
Optimization terminated successfully. (Exit mode 0)
Current function value: 6.39496039618e-07
Iterations: 9
Function evaluations: 60
Gradient evaluations: 9
Took:0.344144268831 min.
Out[12]:
fun: 6.3949603961751021e-07
jac: array([-0.00169004, -0.00146682, 0.001104 , 0.00482058, 0. ])
message: 'Optimization terminated successfully.'
nfev: 60
nit: 9
njev: 9
status: 0
success: True
x: array([ 1.35001977, -0.58735602, -0.09799538, 0.97669182])
'''
# [ 1.3334765 , -0.58013206, 0.034585 , 1.03876482])
# NEW:
# 1.33501652, -0.5832553 , -0.31480513, 1.03965954
#TRUE:
#
# mystery_theta = {'alpha':1.5,
# 'beta':-0.5,
# 'mu':0.5, # So note that mu is not correct.
# 'sigma':1.2} # Still not gotten great.
# =======================================================
# REsults with data *not* standardized version:
'''
Optimization terminated successfully. (Exit mode 0)
Current function value: 5.65339794279e-08
Iterations: 10
Function evaluations: 63
Gradient evaluations: 10
Took:0.227895700932 min.
Out[20]:
fun: 5.6533979427881372e-08
jac: array([ -8.66875749e-04, 3.26018435e-04, -6.81478371e-05,
-9.61340689e-04, 0.00000000e+00])
message: 'Optimization terminated successfully.'
nfev: 63
nit: 10
njev: 10
status: 0
success: True
x: array([ 1.35006237, -0.58585084, -0.14527886, 0.9652678 ])
'''
# REsults with data standardized by mean, stdev, version:
'''
Optimization terminated successfully. (Exit mode 0)
Current function value: 9.45908866437e-09
Iterations: 18
Function evaluations: 112
Gradient evaluations: 18
Took:0.31707701683 min.
Out[16]:
fun: 9.4590886643706467e-09
jac: array([ 4.14401579e-04, -5.15464206e-05, 1.52895583e-04,
1.10592089e-04, 0.00000000e+00])
message: 'Optimization terminated successfully.'
nfev: 112
nit: 18
njev: 18
status: 0
success: True
x: array([ 1.35011253, -0.5861334 , 0.06476111, 0.21112074])
'''
In [ ]:
theta_minus_theta
In [ ]:
# Find the quandratic value:
val = np.dot(theta_minus_theta, W0)
val = np.dot(val, theta_minus_theta)
val
In [ ]: