![]() |
|
The apop_model is intended to provide a consistent expression of any model that (implicitly or explicitly) expresses a likelihood of data given parameters, including traditional linear models, textbook distributions, Bayesian hierarchies, microsimulations, and any combination of the above. The unifying feature is that all of the models act over some data space and some parameter space (in some cases one or both is the empty set), and can assign a likelihood for a fixed pair of parameters and data given the model. This is a very broad requirement, often used in the statistical literature. For discussion of the theoretical structures, see A Useful Algebraic System of Statistical Models (PDF).
This page includes:
Users are encouraged to always use models via the helper functions, like apop_estimate or apop_cdf. The helper functions do some boilerplate error checking, and are where the defaults are called: if your model has a log_likelihood
method but no p
method, then apop_p will use exp(log_likelihood
). If you don't give an estimate
method, then apop_estimate
will call apop_maximum_likelihood.
So the game in writing a new model is to write just enough internal methods to give the helper functions what they need. In the not-uncommon best case, all you need to do is write a log likelihood function.
Here is how one would set up a model that could be estimated using maximum likelihood:
where data
is the input data, and m
is the parametrized model (i.e. your model with a parameters
element set by the caller). This function will return the value of the log likelihood function at the given parameters.
vsize
, msize1
, and msize2
elements specify the shape of the parameter set. For example, if it's three numbers in the vector, then set .vsize=3
and omit the matrix sizes. The default model prep routine will call new_est->parameters = apop_data_alloc(vsize, msize1, msize2)
. dsize
element is the size of one random draw from your model. -1
for vsize
, msize(1|2)
, or dsize
. If the allocation is exceptional in a different way, then you will need to allocate parameters by writing a custom prep
method for the model. You already have more than enough that something like this will work (the dsize
is used for random draws):
Once that baseline works, you can fill in other elements of the apop_model as needed.
For example, if you are using a maximum likelihood method to estimate parameters, you can get much faster estimates and better covariance estimates by specifying the dlog likelihood function (aka the score):
The score has to be registered (see below) using
Your model may need additional settings or auxiliary information to function, which would require associating a model-specific struct with the model.
Before getting into the detail of how to make model-specific groups of settings work, note that there's a lightweight method of storing sundry settings, so in many cases you can bypass all of the following.
The apop_model structure has a void
pointer named more
which you can use to point to a model-specific struct. If more_size
is larger than zero (i.e. you set it to your_model.more_size=sizeof(your_struct)
), then it will be copied via memcpy
by apop_model_copy, and freed by apop_model_free. Apophenia's estimation routines will never impinge on this item, so do what you wish with it.
The remainder of this subsection describes the information you'll have to provide to make use of the conveniences described to this point: initialization of defaults, smarter copying and freeing, and adding to an arbitrarily long list of settings groups attached to a model. You will need four items: a typedef for the structure itself, plus init, copy, and free functions. This is the sort of boilerplate that will be familiar to users of object oriented languages in the style of C++ or Java, but it's really a list of arbitrarily-typed elements, which makes this feel more like LISP. [And being a reimplementation of an existing feature of LISP, this section will be macro-heavy.]
ysg_settings
, with a dataset, its two sizes, and an owner-of-data marker. ysg
stands for Your Settings Group; replace that substring with your preferred name in every instance to follow.The first item is a familiar structure definition. The last line is a macro that declares the three functions below. This is everything you would need in a header file, should you need one. These are just declarations; we'll write the actual init/copy/free functions below.
The structure itself gets the full name, ysg_settings
. Everything else is a macro, and so you need only specify ysg
, and the _settings
part is filled in. Because of these macros, your struct
name must end in _settings
.
If you have an especially simple structure, then you can generate the three functions with these three macros in your .c
file:
These macros generate appropriate functions to do what you'd expect: allocating the main structure, copying one struct to another, freeing the main structure. The spaces after the commas indicate that no special code gets added to the functions that these macros generate.
You'll never call these funtions directly; they are called by Apop_settings_add_group, apop_model_free, and other model or settings-group handling functions.
Now that initializing/copying/freeing of the structure itself is handled, the remainder of this section will be about how to add instructions for the struture internals, like data that is pointed to by the structure elements.
NULL
. In most cases, though, you will need a new line declaring a default for every element in your structure. There is a macro to help with this too. These macros will define for your use a structure named in
, and an output pointer-to-struct named out
. Continuing the above example:Now, Apop_settings_add(a_model, ysg, .size1=100)
would set up a group with a 100-by-10 data set, and set the owner bit to one.
in
for your use. Continuing the example:With those three macros in place and the header as above, Apophenia will treat your settings group like any other, and users can use Apop_settings_add_group to populate it and attach it to any model.
For any given function (e.g., entropy, the dlog likelihood, Bayesian updating), there is probably a special case for well-known models like the Normal distribution. Rather than any procedure that could have a special-case calculation to the apop_model
struct, functions may maintain a registry of models and associated special-case procedures.
This subsection will discuss how to add a function to an existing vtable.
_vtable_add
function to add the function and associate it with the given model. For example, to add a Beta-binomial routine named betabinom
to the registry of Bayesian updating routines, use apop_update_vtable_add(betabinom, apop_beta, apop_binomial)
. apop_update_hash
takes in two models and calculates the hash based on the address of the prior's draw
method and the likelihood's log_likelihood
or p
method. Thus, a vtable lookup for new models that re-use the same methods (at the same addresses in memory) will still find the same special-case function. apop_update_vtable_drop(apop_beta, apop_binomial)
. You can guarantee that a method will not be re-added by following up the _drop
with, e.g., apop_update_vtable_add(NULL, apop_beta, apop_binomial)
. ..._vtable_add
are typically placed in the prep
method of the given model, thus ensuring that the auxiliary functions are registered after the first time the model is sent to apop_estimate.This overview will not go into detail about setting up a new vtable. Briefly:
draw
, log_likelihood
, and methods
of the model. A model where these elements are identical but the name is changed will still match.The remainder of this page covers the detailed expectations regarding the elements of the apop_model structure. I begin with the data (non-function) elements, and then cover the method (function) elements. Some of the following will be requirements for all models and some will be advice to authors; I use the accepted definitions of "must", "shall", "may" and related words.
->matrix
element of the apop_data set sent to model methods. Your likelihood
, p
, cdf
, and estimate
routines must accept data as a single row of the matrix of the apop_data set for such functions to work.ols_shuffle
to convert a matrix where the first column is the dependent variable to a data set with dependent variable in the vector and a column of ones in the first matrix column. By checking for a vector, the prep function knows whether to do the shuffling or not. Most univariate distributions take each scalar element as a separate data point; having one data point per row is a special case.prep
method of the model; see below. Given the model m
and its elements m.vsize
, m.msize1
, m.msize2
, functions that need to allocate a parameter set will do so via apop_data_alloc(m.vsize, m.msize1, m.msize2)
.
.vsize,
.msize1, or
.msize2 to -1
, then the default prep method will set that size to the number of columns in the input data. This is what you want for regression methods, where there is one parameter per independent variable.<info>
is typically a list of scalars. Nothing is guaranteed, but the elements may include:For those elements that require a count of input data, the calculations assume each row in the input apop_data set is a single datum.
Get these via, e.g., apop_data_get(your_model->info, .rowname="log likelihood")
. When writing for any arbitrary function, be prepared to handle NaN
, indicating that the element is not calculated or saved in the info page by the given model.
predict
table. The table has these rows: For OLS-type estimations, each row corresponds to the row in the original data. For filling in of missing data, the elements may appear anywhere, so the row/col indices are essential.
In object-oriented jargon, settings groups are the private elements of the data set, to be pulled out in certain contexts, and ignored in all others. Therefore, there are no rules about internal use. The more
element of the apop_model provides a lightweight means of attaching an arbitrary struct to a model. See Writing new settings groups above for details.
long double your_p_or_ll(apop_data *d, apop_model *params)
. ->parameters
element. NULL
and set the model's error
element to 'p'
if they are. NaN
on errors. If an error in the input model is found, the function may set the input model's error
element to an appropriate char
value. p
methods, it must be the case that log(p(d, m))
equals log_likelihood(d, m)
for all d
and m
.void your_prep(apop_data *data, apop_model *params)
. vsize
, msize1
, or msize2
are -1, then the prep function will set them to the width of the input data. dsize
is -1, then the prep function shall set it to the width of the input data. parameters
element is not allocated, the function shall allocate it via apop_data_alloc(vsize, msize1, msize2)
(or equivalent). data
pointer shall be set to point to the input data. data
. info
element shall be allocated and its title set to "<Info>". void your_estimate(apop_data * data, apop_model *params)
. parmaeters
hold garbage (as in a malloc
without a subsequent assignment to the malloc
-ed space). estimate
function with the prepeped copy. parameters
of the input model. For consistency with other models, the estimate should be the maximum likelihood estimate, unless otherwise documented. <Info>
page may be filled with data. For scalars like log likelihood and AIC, use apop_data_add_named_elmt. estimate
routine; any changes to the data made by estimate
must be documented. error
element to a single character. Documentation should include the list of error characters and their meaning.void your_draw(double out, gsl_rng r, apop_model *params)
paramters
are set, via apop_estimate or apop_model_set_parameters. The author of the draw method should check that parameters
are not NULL
and fill the output with NaNs if necessary parameters are not set. double
of length dsize
; user is expected to make sure that there is adequate space. User also inputs a gsl_rng
, already allocated (probably via apop_rng_alloc). p
method. Data shall be reduced to a single vector via apop_data_pack if it is not already a single vector.long double your_cdf(apop_data *d, apop_model *params)
. paramters
are set, via apop_estimate or apop_model_set_parameters. The author of the CDF method should check that parameters
are not NULL
and return NaN if necessary parameters are not set. matrix
of the input apop_data set (as per a draw produced using the draw
method). May accept other formats. apop_cdf_settings
group may be added to the model. See the apop_cdf function for details of its use.long double your_constraint(apop_data *data, apop_model *params)
. parameters
are set, via apop_estimate, apop_model_set_parameters, or the internals of an MLE search. The author of the constraint method should check that parameters
are not NULL
and return NaN if necessary parameters are not set. parameters
in the input model to a constraint-satisfying value, and (2) return the distance between the input parameters and what you've moved the parameters to. The choice of within-bounds parameters and distance function is left to the author of the constraint function.