Skip to content

Instantly share code, notes, and snippets.

@josef-pkt
Created February 25, 2012 23:13
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save josef-pkt/1911544 to your computer and use it in GitHub Desktop.
Save josef-pkt/1911544 to your computer and use it in GitHub Desktop.
Partial linear model - estimated with concentrated least square
# -*- coding: utf-8 -*-
"""Partial Linear Model, with parametric nonlinear part
Created on Sat Feb 25 17:03:13 2012
Author: Josef Perktold
License: BSD-3 (statsmodels)
warning: first draft, not much checked yet,
dangerous: example uses same variable names as class
"""
import numpy as np
from scipy import optimize
from scikits.statsmodels.regression.linear_model import OLS
class PartialLinear(object):
def __init__(self, y, x_nonlin, x_lin, n_params_nonlin):
self.y = y
self.x_nonlin = x_nonlin
self.x_lin = x_lin
self.n_params_nonlin = n_params_nonlin
self.n_params_lin = x_lin.shape[1]
def predict_nonlin(self, p_nonlin, x_nonlin=None):
if x_nonlin is None:
x_nonlin = self.x_nonlin
return 1. / (1 + np.exp(np.dot(x_nonlin, p_nonlin)))
def predict_lin(self, p_lin, x_lin=None):
if x_lin is None:
x_lin = self.x_lin
return np.dot(x_lin, p_lin)
def predict(self, params, x_nonlin=None, x_lin=None):
p1 = params[:self.n_params_nonlin]
p2 = params[-self.n_params_lin:]
ypred = self.predict_nonlin(p1, x_nonlin)
ypred += self.predict_lin(p2, x_lin)
return ypred
def error(self, params):
return self.y - self.predict(params)
def error_concentrated(self, p_nonlin):
ypred_nonlin = self.predict_nonlin(p_nonlin)
ydiff = self.y - ypred_nonlin
self.result_ols = OLS(ydiff, self.x_lin).fit() #attach last result
resid = self.result_ols.resid
return resid
def fit_concentrated(self, start_value=None):
if start_value is None:
start_value = np.ones(self.n_params_nonlin)
return optimize.leastsq(self.error_concentrated, start_value)
def fit_full(self, start_value=None):
if start_value is None:
start_value = np.ones(self.n_params_nonlin + self.n_params_lin)
return optimize.leastsq(self.error, start_value)
nobs = 1000 #use large sample to see if we get close to true params
sige = 0.5
x_lin = np.column_stack((np.ones(nobs), np.random.randn(nobs, 2)))
x_nonlin = np.column_stack((np.ones(nobs), np.random.randn(nobs, 3)))
p_lin_true = np.ones(2+1) * 0.5
p_nonlin_true = np.ones(3+1)
y_true = 1. / (1 + np.exp(np.dot(x_nonlin, p_nonlin_true)))
y_true += np.dot(x_lin, p_lin_true)
y_obs = y_true + sige * np.random.randn(nobs)
mod = PartialLinear(y_obs, x_nonlin, x_lin, x_nonlin.shape[1])
print mod.fit_concentrated(np.array([1.,0,0,0]))
p_c = mod.fit_concentrated()
print p_c
print mod.fit_full(np.array([1., 0, 0, 0, 1, 0, 0])) #constant only
start_value = np.concatenate((p_c[0], mod.result_ols.params))
print mod.fit_full(start_value)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment