Alright, I'm having an insane amount of difficulty for what seems like a simple concept. I need to generate a bunch of basis functions for a curve that underlies some simulated Poisson distributed data and then run Newton's method on it to fit the log-likelihood (Poisson regression). I understand the math behind Newton's method enough to code that up and fit the data. However, I'm having a lot of difficulty just trying to generate a set of basis functions for which the equation $-\mathbf{X}^{\mathrm{T}}\cdot\mathrm{diag}\left(\exp\left(\mathbf{X}\cdot\boldsymbol\theta\right)\right)\cdot\mathbf{X}$ results in a tridiagonal matrix, where $$\mathbf{X}\cdot\boldsymbol\theta= \begin{bmatrix} f_{1}\left(\mathbf{t}_1\right) & f_{2}\left(\mathbf{t}_1\right) & \cdots & f_{d}\left(\mathbf{t}_1\right)\\ f_{1}\left(\mathbf{t}_2\right) & f_{2}\left(\mathbf{t}_2\right) & \cdots & f_{d}\left(\mathbf{t}_2\right) \\ \vdots & \vdots & \ddots & \vdots\\ f_{1}\left(\mathbf{t}_n\right) & f_{2}\left(\mathbf{t}_n\right) & \cdots & f_{d}\left(\mathbf{t}_n\right) \end{bmatrix} \begin{bmatrix} \theta_1\\ \theta_2\\ \vdots\\ \theta_d \end{bmatrix},$$ and $\mathrm{diag}$ results in a matrix with the elements of its argument on the diagonal of an otherwise zero-filled matrix.
Here's what I've got so far:
function ddl = fastnewton(d, varargin)
%% boilerplate for optional function arguments
p = inputParser;
% required arguments
p.addRequired('d', @isposintscalar)
% optional arguments
p.addOptional('n', 1000, @isposintscalar)
% keyword arguments
p.addParamValue('order', 3, @isposintscalar)
p.addParamValue('tminmax', [-10, 10], ...
@(x) (isvector(x) || iscell(x)) && length(x) == 2)
p.addParamValue('weightfunc', @(varargin) abs(randn(varargin{:})), ...
@(x) isa(x, function_handle))
% parse the args and get the results
p.parse(d, varargin{:})
r = p.Results;
d = r.d;
n = r.n;
order = r.order;
if iscell(r.tminmax)
r.tminmax = cell2mat(r.tminmax);
end
tminmax = r.tminmax;
weightfunc = r.weightfunc;
%% make a time vector between tmax and tmin, inclusive
% tminmax == [-10, 10]
tmin = tminmax(1);
tmax = tminmax(2);
t = linspace(tmin, tmax, n); % n == 1000
%% make some random weights using the 'weightfunc' function
% weightfunc == @(varargin) abs(randn(varargin{:}))
theta = weightfunc(1, d);
%% make the bump funcs
breaks = linspace(tmin, tmax, d + order);
coeffs = theta;
pp = spmak(breaks, diag(coeffs));
%% compute the value of the function to be estimated
Fbig = spval(pp, t);
F = sum(Fbig);
%% initialize the sample vector
y = poissrnd(exp(F));
subplot(211)
hold on
plot(t, y, 'LineWidth', 1)
plot(t, Fbig, 'LineWidth', 2)
plot(t, F, 'r', 'LineWidth', 4)
hold off
legend({'y', '', 'sum(F)'})
%% plotting
lw = 3;
opts = {'LineWidth', lw};
subplot(212)
hold on
plot(t, y, 'b', opts{:})
plot(t, F, 'g', opts{:})
hold off
axis tight
X = Fbig';
ddl = sparse(-X' * diag(exp(F)) * X);
figure
hold on
plot(t, y, 'b', 'LineWidth', 1)
plot(t, F, 'r', 'LineWidth', 2);
hold off
figure
spy(ddl)
end
You can run this code as follows: fastnewton(30, 'order', 2). That will give you 30 basis functions of order 2 which will give you a tridiagonal matrix. However, if you run fastnewton(30, 'order', 3) you get nice, smooth basis functions but a matrix with a bandwith of 3 instead, which is not tridiagonal. MATLAB's cute little spy function gives you a graphical representation of the sparsity of a matrix.
All I really need is some way to adjust the basis functions so that they are smooth AND the equation above results in a tridiagonal matrix, then I can do the rest myself. Thanks!
