Dynamic Display of Changing Posterior in Bayesian Survival
Analysis: The Software
Hani J. Doss
Department of Statistics
The Ohio State
University
Columbus, OH 43210
and
B. Narasimhan
Department of Statistics
Stanford University
Stanford, CA 94305
Revision: 1.34 of Date: 1998/07/02 17:16:33
Please note that this is an imperfect automatic translation of
the file. The math typesetting is most likely clobbered.
Table of Contents
- Copyright
- Introduction
- The Software
- The Master Object Prototype
- The :identifier Method
- The :number-of-markov-chains Method
- The :number-of-points Method
- The :number-of-months Method
- The :number-of-data-values Method
- The :summary-data Method
- The :indicator-counts Method
- The :log-constants-of-proportionality Method
- The :slaves Method
- The :initially-specified-hyperparameter-values Method
- The :hyperparameter-names Method
- The :current-hyperparameter-values Method
- The :number-of-hyperparameters Method
- The :hyperparameters-used-in-markov-chains Method
- The :hyperparameter-ranges Method
- The :hyperparameter-sliders Method
- The :log-mixture-density Method
- The :importance-weights Method
- The :loglik Method
- The :compute-log-hmix Method
- The :calc-weights Method
- The :isnew Method
- The :process-run-file Method
- The :process-frequency-table Method
- The :graphical-interface Method
- The :create-run-file Method
- The :synchronize Method
- The :consolidate-computation Method
- The :reset Method
- The :effective-sample-size Method
- The :print-all-statistics Method
- The :labelled-hyperparameter-values Method
- The :statistics Method
- The :statistics-print-formats Method
- The :statistics-labels Method
- The :toggle-timing Method
- The :superimpose Method
- The :close Method
- Defaults for Master
- The Slave Object Prototype
- The :isnew Method
- The :redraw-background Method
- The :redraw-statistics Method
- The :print-summary Method
- The :close Method
- Defaults for Slave
- The C Programs
- Beta Functions
- Global Variables
- Initialization Routine
- The C Equivalent of :f Method
- The C Equivalent of :logp Method
- The C Equivalent of :loglik Method
- The C Equivalent of :htheta-over-hmix Method
- The C Equivalent of :compute-log-hmix Method
- The C Equivalent of :calc-weights Method
- The C Equivalent of :compute-statistics Method
- The C Equivalent of :compute-law-of-f-of-t Method
- The C Equivalent of :compute-mean-of-fbar-of-t Method
- The :consolidate-computation Method in C
- Installation Information
- Improvements needed
- Acknowledgement
- Index of Code Chunks
- Index of Identifiers
Abstract
We consider the problem of estimating an unknown distribution
function F in the presence of censoring under the conditions
that a parametric model is believed to hold approximately. We
use a Bayesian approach, in which the prior on F is a mixture
of Dirichlet distributions. A hyperparameter of the prior
determines the extent to which this prior concentrates its mass
around the parametric family. A Gibbs sampling algorithm to
estimate the posterior distributions of the parameters of
interest is reviewed. An importance sampling scheme enables us
to use the output of the Gibbs sampler to very quickly
recalculate the posterior when we change the hyperparameters of
the prior. The calculations can be done sufficiently fast to
enable the dynamic display of the changing posterior as the
prior hyperparameters are varied.
This paper provides a literate program completely documenting the
code for performing the dynamic graphics.
[*]
We begin with our usual copyright.
<Copyright>= (U->)
;;;
;;; $Revision: 1.34 $ of $Date: 1998/07/02 17:16:33 $
;;;
;;; Copyright (C) 1994, 1995, 1998. Doss and Narasimhan
;;;
;;; Hani J. Doss (doss@stat.ohio-state.edu) and
;;; B. Narasimhan (naras@stat.stanford.edu)
;;;
;;; This program is free software; you can redistribute it and/or modify
;;; it under the terms of the GNU General Public License as published by
;;; the Free Software Foundation; either version 2 of the License, or
;;; (at your option) any later version.
;;;
;;; This program is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
;;; GNU General Public License for more details.
;;;
;;; You should have received a copy of the GNU General Public License
;;; along with this program; if not, write to the Free Software
;;; Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
;;;
Defines copyright (links are to index).
Introduction
[*]
This document is a literate program implementing the theory
described in our paper [cite doss:nara:1998b]. A literate
program is a program written in a style that makes it easy for
humans to read, understand and modify. For more information on
Literate Programming, see [cite knut:1992]. A quicker introduction
is available on the World Wide Web; see [cite lee:1994]. This
document uses the Noweb [cite rams:1994], [cite rams:1994a],
[cite rams:1995] literate programming tools. Although it is not required,
we recommend that any serious user of this software have the
Noweb tools installed. Noweb tools, besides being
free, are extremely easy to install and require no special expertise
other than basic knowledge about how TeX or LaTeX work. Having
Noweb allows one to take full advantage of our
software---syntax errors will be minimized and the code we have
written can be reused with the user's modifications spliced into
place automatically by Noweb.
We expect any serious user of our software to read the original paper
[cite doss:nara:1998b], a copy of which is included in the software
distribution.
This document is available in three forms: Postscript, PDF and HTML.
All versions are accessible from the web pages of the authors.
We wish to remark that the software only does sensitivity analysis.
No general facility is provided for generating observations from
Markov chains. Indeed, since the range of models for which MCMC
methods are applicable is large and such methods most likely involve
problem-specific issues, it is our opinion that building such a
supertool, if it is at all possible, is a non-trivial task. However,
the Fortran program used in generating the output for our example is
included along with this software and can be used for models similar
to ours. Of course, any appropriate method may be used to generate
the samples as long as the output is available in a form usable by
our software. The requirements on the data that can be used with our
software are spelt out below.
Corresponding to each Markov chain output, there must be two files
with the extensions .in (input file) and .out
(output file). For example, mc1.in and mc1.out.
The input file must have the following structure. The first four items
in the file can be anything, string or number, either on a single line
or any conceivable combination of lines. The next three items
must be the shape of the Gamma distribution on 
, the
scale of the Gamma distribution on 
---the parametrization for
shape a and scale b is proportional to
---and 
. The next
three values values following these quantities can be anything, but
the one following it should be the number of data points, or
sets. Nothing else is read from the input file.
The output file must have the following structure for each data point
generated by the Markov chain. The value of theta must be followed
by the number of distinct values of the data points, which must be
followed by a frequency table of the actual data value and the
corresponding frequency. The layout of the values on lines does not
matter as long as at least a single white space delimits values. If
this structure is violated, errors will result. A peek at the data
files included with this software will help the reader.
A note on performance. The calculations involved in reweighting are
non-trivial and require a reasonably powerful computer for smooth
performance. The efficiency can be improved by dynamically loading C
programs that compute various quantities. The version of software
described here does so by default. An older version that does not use
dynamic loading which is available upon request from the authors.
Without dynamic loading, the performance is very bad indeed.
Dynamic libraries for Windows and Macintosh are provided. Suggestions
for various Unix platforms are also provided.
Section [->] has more details.
It is assumed that a proper installation of Lisp-Stat
described in [cite tier:1990] is available. The version number on
Lisp-Stat should be 3.52.0 or higher since dynamic loading
uses the new shared library mechanism.
Some additions planned for the future are listed in
Section [->].
The Software
[*]
The software consists of the following components.
<*>=
<Copyright>
(require "utility")
(require "call-by-reference")
(defpackage "BSA" (:use "XLISP" "USER" "UTILITY" "CALL-BY-REFERENCE"))
(in-package "BSA")
(import
'(user::rseq
user::standard-deviation))
<The Master Prototype>
<The Slave Prototype>
(export
'(master-proto))
[*]
The master prototype master-proto inherits from dialog-proto of Lisp-Stat
and contains a number of slots. A rich set of methods facilitate
interaction with the master.
<The Master Prototype>= (<-U)
(defproto master-proto
'(identifier number-of-markov-chains
number-of-points
number-of-data-values
number-of-months
data-file-names
summary-data
indicator-counts
log-constants-of-proportionality
slaves
hyperparameters-used-in-markov-chains
initially-specified-hyperparameter-values
current-hyperparameter-values
hyperparameter-names
log-mixture-density
importance-weights
hyperparameter-ranges
hyperparameter-sliders
work-space
shared-library
density-abscissae
density-ordinates
expectation-abscissae
expectation-ordinates
standard-deviation-ordinates
statistics
statistics-print-formats
statistics-labels
number-of-slider-stops
superimpose
timing
timing-button
lazy)
() dialog-proto
"The Master prototype. Creates and manipulates a harem of slaves.")
<Methods for Master Prototype>
<Defaults for Master>
Defines current-hyperparameter-values, data-file-names, density-abscissae, density-ordinates, expectation-abscissae, expectation-ordinates, hyperparameter-names, hyperparameter-ranges, hyperparameter-sliders, hyperparameters-used-in-markov-chains, identifier, importance-weights, indicator-counts, initially-specified-hyperparameter-values, lazy, log-constants-of-proportionality, log-mixture-density, master-proto, number-of-data-values, number-of-markov-chains, number-of-months, number-of-points, number-of-slider-stops, shared-library, slaves, standard-deviation-ordinates, statistics, statistics-labels, statistics-print-formats, summary-data, superimpose, timing, work-space (links are to index).
Throughout, we shall use m for the number of Markov chains, and n for
the number of data points in the output of each Markov chain. Note
that arrays are indexed from 0 so that index i refers to the
(i+1)-th position.
Here is a description of all the slots in master-proto.
- identifier
- holds a string that is used to identify
instances of the object. This string also helps in repeating any
interesting run---all inputs needed for the run are saved in a file
whose name is generated by adding an extension
.run to the
string. In addition, a file with the extension .lsp is also
created so that invoking xlispstat on that file with
automatically recreate everything without annoying questions.
- number-of-markov-chains
- holds the number of Markov
chains (m)to be used in the exploration.
- number-of-points
- holds the number of data points from
each Markov chain (n) to be used in the exploration.
- number-of-data-values
- holds the number of data values,
that is the number of sets.
- number-of-months
- holds the number of months for
which the law of F(t) is computed. Default is given by the variable
*default-number-of-months*.
- data-file-names
- holds a list of the data files
corresponding to the Markov chains, without the extension. Dimension
is m.
- summary-data
- This slot holds the summary data from all
the Markov chains. The data is stored as a 2-d array with each row
holding (d_x, theta), where d_x is the the number of distinct
values of the data points. The first set of n pairs correspond to
the first Markov chain, the next set to the second and so on, which
means that the indices have to be suitably translated to access
values. There are mn pairs in all. Note that we only handle
balanced data, i.e., same n for all Markov chains.
- indicator-counts
- This slot is a 2-d array of size
mn×2 containing the a count of the number of data-values
less than t, the number of months. The number of columns in this
array is indicated by the slot
number-of-months. If this value
is 61, say, then one would be able to calculate the law of F(t)
for t ranging from 0 to 60 months.
- log-constants-of-proportionality
- holds a vector of
estimates of the logs of the constants of proportionality,
(C(kappa_i))^-1, in the paper. The first constant is
implicitly -1 and so the dimension of this slot is m-1.
- slaves
- stores a list of slaves who need to be informed
of changes in the hyperparameter values. Each slave is an instance
of slave-proto.
- hyperparameters-used-in-markov-chains
- stores an m by
3 array holding the values of hyperparameters at which the Markov
chains were run. If a is the array by a, then a[j,i] is the
value of the i-th hyperparameter for the j-th Markov chain, 0
<=j < m and 0 <=i < 3. The three hyperparameters in every
row are alpha(R), the shape and scale of the Gamma prior
on theta respectively.
- initially-specified-hyperparameter-values
- is a
4×1 array of initial values for the hyper-parameters at
which the exploration should begin. These values serve are used in
starting and restarting the exploration. The one extra
hyperparameter is the value of time t in F(t).
- current-hyperparameter-values
- is a 4×1 array
of the current values of the three hyperparameters. These are the
quantities that are changed by the user via sliders.
- hyperparameter-names
- holds a list of names for the
hyper-parameters. Default is Alpha, Theta Shape,
Theta Scale, and Time (Months).
- log-mixture-density
- is a an array holding
logh_mix(x) at the data points. Size is mn. See
[cite geye:1994]. For speed in dynamic graphics, this quantity is
calculated once and saved.
- importance-weights
- holds the importance weights
w[i], 0 <=i < nm, used for the Reweighting Mixtures (RM)
scheme. See [cite doss:nara:1994] and [cite geye:1994].
- hyperparameter-ranges
- stores a list of hyperparameter
ranges to be used in exploration. For internal use only.
- hyperparameter-sliders
- is a list of slider objects for
internal use only.
- work-space
- is a slot used for storing results from
dynamically loaded C routines.
- shared-library
- holds a library for shared handle when
dynamic loading is used. For internal use only.
- density-abscissae
- is a vector of abscissa values
between 0 and 1 at which the density of F(t) will be
plotted. This is an array of length
*default-number-of-plot-stops*.
- density-ordinates
- is a vector of values of the
density of F(t). This is an array of length
*default-number-of-plot-stops*.
- expectation-abscissae
- is a vector of values of
E(F(t)) for t ranging from 0 to
*default-number-of-months* (minus 1, as we start at 0).
- expectation-ordinates
- is a vector of values of
E(F(t)) for t ranging from 0 to
*default-number-of-months* (minus 1, as we start at 0).
- standard-deviation-ordinates
- is a vector of values of
sigma_F(t).
- statistics
- is a vector of 3 values that will hold
E(F(t)), sigma_F(t) and the effective
sample size. The effective sample size is calculated using the
formula
Effective Sample Size = m * n / (1 + cv(W)^2),
where cv(W) is the coefficient of variation of the importance
weights W. See, for example, [cite kong:liu:wong:1994].
- statistics-print-formats
- is a list of lists indicating
the format to be used in printing the statistics values.
- statistics-labels
- is a list of strings (labels).
- number-of-slider-stops
- stores the number of slider
stops for hyperparameters.
- superimpose
- toggles superimposition on and
off. Default is off.
- timing
- signifies if timing is needed.
- timing-button
- is a button for toggling timing on and
off.
- lazy
- is a slot used for efficient synchronization. It
is for internal use by programs and the user shouldn't mess with it.
The methods for master-proto can be broken down as follows.
<Methods for Master Prototype>= (<-U)
<The Master :identifier Method>
<The Master :number-of-markov-chains Method>
<The Master :number-of-points Method>
<The Master :number-of-months Method>
<The Master :number-of-data-values Method>
<The Master :summary-data Method>
<The Master :indicator-counts Method>
<The Master :log-constants-of-proportionality Method>
<The Master :slaves Method>
<The Master :initially-specified-hyperparameter-values Method>
<The Master :hyperparameter-names Method>
<The Master :current-hyperparameter-values Method>
<The Master :number-of-hyperparameters Method>
<The Master :hyperparameters-used-in-markov-chains Method>
<The Master :hyperparameter-ranges Method>
<The Master :hyperparameter-sliders Method>
<The Master :log-mixture-density Method>
<The Master :importance-weights Method>
<The Master :loglik Method>
<The Master :compute-log-hmix Method>
<The Master :calc-weights Method>
<The Master :isnew Method>
<The Master :graphical-interface Method>
<The Master :process-run-file Method>
<The Master :process-frequency-table Method>
<The Master :create-run-file Method>
<The Master :synchronize Method>
<The Master :consolidate-computation Method>
<The Master :reset Method>
<The Master :effective-sample-size Method>
<The Master :print-all-statistics Method>
<The Master :labelled-hyperparameter-values Method>
<The Master :statistics Method>
<The Master :statistics-print-formats Method>
<The Master :statistics-labels Method>
<The Master :superimpose Method>
<The Master :toggle-timing Method>
<The Master :close Method>
Some of these methods are mere accessor and modifier methods for the
slots and we can get them easily out of the way.
The :identifier Method
[*]
<The Master :identifier Method>= (<-U)
(defmeth master-proto :identifier (&optional name)
"Method args: (&optional name)
Sets or retrieves the identifier slot."
(if name
(setf (slot-value 'identifier) name)
(slot-value 'identifier)))
Defines :identifier (links are to index).
The :number-of-markov-chains Method
[*]
<The Master :number-of-markov-chains Method>= (<-U)
(defmeth master-proto :number-of-markov-chains ()
"Method args: ()
Returns m, the number of Markov chains used."
(slot-value 'number-of-markov-chains))
Defines :number-of-markov-chains (links are to index).
The :number-of-points Method
[*]
<The Master :number-of-points Method>= (<-U)
(defmeth master-proto :number-of-points ()
"Method args: ()
Returns n, the number of data points."
(slot-value 'number-of-points))
Defines :number-of-points (links are to index).
The :number-of-months Method
[*]
<The Master :number-of-months Method>= (<-U)
(defmeth master-proto :number-of-months ()
"Method args: ()
Returns n, the number of months for which F(t) is computed."
(slot-value 'number-of-months))
Defines :number-of-months (links are to index).
The :number-of-data-values Method
[*]
<The Master :number-of-data-values Method>= (<-U)
(defmeth master-proto :number-of-data-values ()
"Method args: ()
Returns n, the number of data values, or sets."
(slot-value 'number-of-data-values))
Defines :number-of-data-values (links are to index).
The :summary-data Method
[*]
<The Master :summary-data Method>= (<-U)
(defmeth master-proto :summary-data ()
"Method args: ()
Returns the summary data containing pairs theta and number of distinct
data points. A 2d array of size mn by 2."
(slot-value 'summary-data))
Defines :summary-data (links are to index).
The :indicator-counts Method
[*]
<The Master :indicator-counts Method>= (<-U)
(defmeth master-proto :indicator-counts ()
"Method args: ()
Returns the data values a 2d array of size mn by number-of-months."
(slot-value 'indicator-counts))
Defines :indicator-counts (links are to index).
The :log-constants-of-proportionality Method
[*]
<The Master :log-constants-of-proportionality Method>= (<-U)
(defmeth master-proto :log-constants-of-proportionality ()
"Method args: ()
Returns a list of values of the log constants of
proportionality for the posteriors for each Markov chain."
(slot-value 'log-constants-of-proportionality))
Defines :log-constants-of-proportionality (links are to index).
The :slaves Method
[*]
<The Master :slaves Method>= (<-U)
(defmeth master-proto :slaves ()
"Method args: () Retrieves the slaves."
(slot-value 'slaves))
Defines :slaves (links are to index).
The :initially-specified-hyperparameter-values Method
[*]
<The Master :initially-specified-hyperparameter-values Method>= (<-U)
(defmeth master-proto :initially-specified-hyperparameter-values ()
"Method args: ()
Returns a list of initially specified values of the hyper
parameters. Used mainly for resetting the hyper parameter values
before dynamic exploration."
(slot-value 'initially-specified-hyperparameter-values))
Defines :initially-specified-hyperparameter-values (links are to index).
The :hyperparameter-names Method
[*]
<The Master :hyperparameter-names Method>= (<-U)
(defmeth master-proto :hyperparameter-names (&optional index)
"Method args: (index)
Returns a string identifying the hyperparameter index or the whole slot."
(if index
(select (slot-value 'hyperparameter-names) index)
(slot-value 'hyperparameter-names)))
Defines :hyperparameter-names (links are to index).
The :current-hyperparameter-values Method
[*]
<The Master :current-hyperparameter-values Method>= (<-U)
(defmeth master-proto :current-hyperparameter-values (&optional i x)
"Method args: (&optional i x)
Returns a list of hyper-parameter values. If i is given and a sequence,
then all values in the list are set. The dimensions must match. If i is
not a sequence, returns the i-th value. If x is specified, sets the
i-th hyperparameter value to x. Note that x is ignored if i is a sequence."
(if x
(if (slot-value 'timing)
(time
(progn
(setf (elt (slot-value 'current-hyperparameter-values) i) x)
(send self :synchronize)))
(progn
(setf (elt (slot-value 'current-hyperparameter-values) i) x)
(send self :synchronize)))
(if i
(if (sequencep i)
(progn
(setf (slot-value 'lazy) t)
(dotimes (j (send self :number-of-hyperparameters))
(send (select (slot-value 'hyperparameter-sliders) j)
:value (select i j)))
(setf (slot-value 'lazy) nil)
(send self :synchronize))
(elt (slot-value 'current-hyperparameter-values) i))
(slot-value 'current-hyperparameter-values))))
Defines :current-hyperparameter-values (links are to index).
The
:number-of-hyperparameters Method
[*]
<The Master :number-of-hyperparameters Method>= (<-U)
(defmeth master-proto :number-of-hyperparameters ()
"Method args: ()
Returns the number of hyperparameters."
(length (slot-value 'current-hyperparameter-values)))
Defines :number-of-hyperparameters (links are to index).
The
:hyperparameters-used-in-markov-chains Method
[*]
<The Master :hyperparameters-used-in-markov-chains Method>= (<-U)
(defmeth master-proto :hyperparameters-used-in-markov-chains ()
"Method args: ()
Returns the hyperparameter-values used in running the Markov
chains. 2D array of size number-of-markov-chains by 3."
(slot-value 'hyperparameters-used-in-markov-chains))
Defines :hyperparameters-used-in-markov-chains (links are to index).
The :hyperparameter-ranges Method
[*]
<The Master :hyperparameter-ranges Method>= (<-U)
(defmeth master-proto :hyperparameter-ranges (&optional ranges)
"Method args: (&optional ranges)
Returns the ranges within which the hyperparameters will be
varied. List of lists. This is calculated as the maximum and minimum
of all the values used in the Markov chains if not given."
(if ranges
(setf (slot-value 'hyperparameter-ranges) ranges)
(if (slot-value 'hyperparameter-ranges)
(slot-value 'hyperparameter-ranges)
(let* ((values (column-list
(slot-value 'hyperparameters-used-in-markov-chains)))
(max (mapcar #'max values))
(min (mapcar #'min values)))
(append (mapcar #'(lambda (x y) (list x y)) min max)
(list (list 0 (1- (slot-value 'number-of-months)))))))))
Defines :hyperparameter-ranges (links are to index).
The :hyperparameter-sliders Method
[*]
<The Master :hyperparameter-sliders Method>= (<-U)
(defmeth master-proto :hyperparameter-sliders ()
"Method args: ()
Returns the sliders associated with the hyperparameters."
(slot-value 'hyperparameter-sliders))
Defines :hyperparameter-sliders (links are to index).
The :log-mixture-density Method
[*]
<The Master :log-mixture-density Method>= (<-U)
(defmeth master-proto :log-mixture-density ()
"Method args: ()
Returns the log of the mixture density at each of the data points.
An array of nm values."
(slot-value 'log-mixture-density))
Defines :log-mixture-density (links are to index).
The :importance-weights Method
[*]
<The Master :importance-weights Method>= (<-U)
(defmeth master-proto :importance-weights ()
"Method args: ()
Retrieves the importance sampling weights. An array
of size number-of-markov-chains by number-of-points."
(slot-value 'importance-weights))
Defines :importance-weights (links are to index).
The :loglik Method
[*]
<The Master :loglik Method>= (<-U)
(defmeth master-proto :loglik (n)
"Returns the log-quasi-likelihood as a function of n. The dimension
of n should be one less than the number of Markov chains."
(let ((result (slot-value 'work-space)))
(call-by-reference-oldcfun "logLikelihood"
(slot-value 'shared-library)
(coerce n '(vector c-double)) result)
(aref result 0)))
Defines :loglik (links are to index).
The :compute-log-hmix Method
[*]
<The Master :compute-log-hmix Method>= (<-U)
(defmeth master-proto :compute-log-hmix ()
"Method args: ()
Computes the log of the mixture density at the various data points."
(call-by-reference-oldcfun "computeLogHmix"
(slot-value 'shared-library)))
Defines compute-log-hmix (links are to index).
The :calc-weights Method
[*]
<The Master :calc-weights Method>= (<-U)
(defmeth master-proto :calc-weights ()
"Method args: ()
Calculates the importance weights."
(call-by-reference-oldcfun "calcWeights" (slot-value 'shared-library)))
Defines :calc-weights (links are to index).
The :isnew Method
[*]
The :isnew method for master-proto is a bit involved.
<The Master :isnew Method>= (<-U)
(defmeth master-proto :isnew (&key identifier?
number-of-markov-chains?
number-of-points?
file-names?
log-constants?)
"Method args: (&key identifier?
number-of-markov-chains?
number-of-points
file-names?
log-constants?)
Please refer to the literate program and the paper for a thorough
discussion."
<Set up slot values for master object>
(send self :calc-weights)
<Create slaves of master object>
<Set up dialog and wait for user input>
<Some final touches>)
Defines :isnew (links are to index).
The :synchronize message at the end of this method forces the
computations anyway, so we avoid the computations now. In other words,
initially, we want to be lazy in performing intensive
computations.
<Set up slot values for master object>= (<-U) [D->]
(setf (slot-value 'work-space)
(make-array 1 :initial-element 0.0 :element-type 'c-double))
(setf (slot-value 'lazy) t)
We need to determine m, the number of chains and n, the number of
data points to be used. If the argument identifier? is given, then
all the necessary inputs can be read from the ``run'' file. Otherwise,
we need to prompt the user for everything.
<Set up slot values for master object>+= (<-U) [<-D->]
(setf (slot-value 'hyperparameter-names)
(list "M(R)" "Theta-Shape" "Theta-Scale" "Time (months)"))
(setf (slot-value 'shared-library)
(shlib::shlib-open *default-shared-lib-name*))
(setf (slot-value 'statistics) (make-array 3 :initial-element 0.0
:element-type 'c-double))
(setf (slot-value 'statistics-print-formats)
*default-statistics-print-formats*)
(setf (slot-value 'statistics-labels)
*default-statistics-labels*)
(if identifier?
(progn
(setf (slot-value 'identifier) identifier?)
(send self :process-run-file :number-of-points? number-of-points?))
(send self :graphical-interface :identifier? identifier?
:number-of-markov-chains? number-of-markov-chains?
:number-of-points? number-of-points?
:file-names? file-names?
:log-constants? log-constants?))
Time to get information on how the exploration is to proceed.
<Set up slot values for master object>+= (<-U) [<-D->]
<Set up Hyperparameter-ranges and stops>
(setf (slot-value 'current-hyperparameter-values)
(copy-vector (send self :initially-specified-hyperparameter-values)))
This code chunk sets up a large dialog box for all the various
quantities. In most cases, the user will just continue without change.
<Set up Hyperparameter-ranges and stops>= (<-U)
(let* ((dialog-items ())
(result nil)
(ranges (select (send self :hyperparameter-ranges) (iseq 3)))
(init-vals (if identifier?
(slot-value 'initially-specified-hyperparameter-values)
(mapcar #'(lambda(x) (* 0.5 (sum x))) ranges)))
(stop-vals (if identifier?
(slot-value 'number-of-slider-stops)
(mapcar #'(lambda(x) (1+ (- (second x) (first x))))
ranges)))
(hypernames (select (slot-value 'hyperparameter-names) (iseq 3)))
(prompt-item
(send text-item-proto :new
(format nil
"You can change the interval between which the~%~
hyperparameters may be varied and the number~%~
of stops in the interval (end points included).~%~
Please note that NO ERROR CHECKING is done!~%~%~
The default settings computed from the data files~%~
are shown below. Please click OK when done.")))
(col-1 (send text-item-proto :new "Hyperparameter"))
(col-2 (send text-item-proto :new "Minimum"))
(col-3 (send text-item-proto :new "Maximum"))
(col-4 (send text-item-proto :new "Stops"))
(col-5 (send text-item-proto :new "Initial"))
(headings (list col-1 col-2 col-3 col-4 col-5))
(widths (map-elements #'send headings :width))
(hnames (mapcar #'(lambda(x)
(send text-item-proto :new x)) hypernames))
(mins (mapcar #'(lambda(x)
(send edit-text-item-proto :new
(format nil "~g" (select x 0))
:text-length 10))
ranges))
(maxs (mapcar #'(lambda(x)
(send edit-text-item-proto :new
(format nil "~g" (select x 1))
:text-length 10))
ranges))
(stops (mapcar #'(lambda(x)
(send edit-text-item-proto :new
(format nil "~g" x)
:text-length 10))
stop-vals))
(inits (mapcar #'(lambda(x)
(send edit-text-item-proto :new
(format nil "~g" x)
:text-length 10))
init-vals))
(abort (send modal-button-proto :new "Abort" :action #'top-level))
(dialog nil)
(ok (send modal-button-proto :new "OK"
:action
#'(lambda()
(list
(map-elements #'send mins :text)
(map-elements #'send maxs :text)
(map-elements #'send stops :text)
(map-elements #'send inits :text))))))
(dotimes (i (length hnames))
(send (select hnames i) :width (select widths 0))
(send (select mins i) :width (select widths 1))
(send (select maxs i) :width (select widths 2))
(send (select stops i) :width (select widths 3))
(send (select inits i) :width (select widths 4)))
(setf dialog-items
(list prompt-item
headings
(list (mapcar #'(lambda (x y z w u) (list x y z w u))
hnames mins maxs stops inits))
(list ok abort)))
(setf dialog (send modal-dialog-proto :new dialog-items))
(send dialog :title "Hyperparameter Ranges and Stops")
(setf result (mapcar #'convert-to-numbers (send dialog :modal-dialog)))
(send self :hyperparameter-ranges
(append
(mapcar #'(lambda(x y) (list x y))
(select result 0) (select result 1))
(list (list 0 (1- (slot-value 'number-of-months))))))
(setf (slot-value 'initially-specified-hyperparameter-values)
(make-array 4 :element-type 'c-double
:initial-contents (append (select result 3)
(list
(truncate
(* 0.5
(slot-value 'number-of-months)))))))
(setf (slot-value 'number-of-slider-stops)
(append (select result 2) (list (slot-value 'number-of-months)))))
Now, make the initial arrays for calling the C functions.
<Set up slot values for master object>+= (<-U) [<-D->]
(setf (slot-value 'density-abscissae)
(make-array *default-number-of-plot-stops*
:initial-contents (rseq 0 1 *default-number-of-plot-stops*)
:element-type 'c-double))
(setf (slot-value 'density-ordinates)
(make-array *default-number-of-plot-stops*
:initial-element 0 :element-type 'c-double))
(let ((number-of-stops (send self :number-of-months)))
(setf (slot-value 'expectation-abscissae)
(make-array number-of-stops
:initial-contents (iseq 0 (1- number-of-stops))
:element-type 'c-double))
(setf (slot-value 'expectation-ordinates)
(make-array number-of-stops
:initial-element 0 :element-type 'c-double))
(setf (slot-value 'standard-deviation-ordinates)
(make-array number-of-stops
:initial-element 0 :element-type 'c-double)))
The current values of the hyperparameters must be those specified
initially.
<Set up slot values for master object>+= (<-U) [<-D->]
(let ((tmp-array (send self :initially-specified-hyperparameter-values)))
(setf (slot-value 'current-hyperparameter-values)
(make-array (length tmp-array) :initial-contents tmp-array
:element-type 'c-double)))
We set up the number of points to be used in Reverse Logistic
Regression if necessary. After that, all the data areas will have been
set up and so we pass addresses of the data-arrays to the C routines.
<Set up slot values for master object>+= (<-U) [<-D->]
(let* ((n (slot-value 'number-of-points))
(nc (min *default-number-of-points* n)))
(unless (or log-constants? identifier?
(= (slot-value 'number-of-markov-chains) 1))
(setf nc
(select (get-tested-value-dialog
(format nil "How many points should I use from each chain for~%~
estimating the constants of proportionality?~%~
Unless you are prepared to wait for a long time, you~%~
should go with the default or less!~%")
:initial nc
:test #'(lambda(x) (and (numberp x) (> x 0) (<= x n)))
:error-message "Invalid entry. Please try again!")
0)))
<Pass data array addresses to C Routines>)
<Pass data array addresses to C Routines>= (<-U)
(call-by-reference-oldcfun "initializeAddress"
(slot-value 'shared-library) (slot-value 'summary-data)
(slot-value 'hyperparameters-used-in-markov-chains)
(slot-value 'current-hyperparameter-values)
(slot-value 'importance-Weights)
(slot-value 'log-mixture-density)
(slot-value 'log-constants-of-proportionality)
(slot-value 'number-of-markov-chains)
(slot-value 'number-of-points)
(slot-value 'number-of-months)
(slot-value 'number-of-data-values)
(slot-value 'indicator-counts)
(slot-value 'statistics)
nc)
If the constants of proportionality have to be estimates, this is the
time to do it after duly notifying the user.
<Set up slot values for master object>+= (<-U) [<-D]
(unless (or log-constants? identifier?)
(unless (= (slot-value 'number-of-markov-chains) 1)
(message-dialog "The Constants of proportionality will now be~%~
be estimated by Reverse Logistic Regression.~%~
This will take a while and only happens once.~%~%~
You will see the results of iterations on the~%~
Console as the maximization is done.~%~%
Click OK to continue.")
<Perform Reverse Logistic Regression>))
(unless identifier?
(send self :compute-log-hmix)
(send self :create-run-file))
<Perform Reverse Logistic Regression>= (<-U)
(flet ((loglik (n) (send self :loglik n)))
(setf initial-guess
(nelmeadmax #'loglik (slot-value 'log-constants-of-proportionality)
:epsilon *default-maximization-tolerance*
:count-limit 50000))
(dotimes (i (length initial-guess))
(setf (aref (slot-value 'log-constants-of-proportionality) i)
(select initial-guess i))))
We are ready to create the slave objects, a plot of the density of
F(t) as well as a plot of E(S(t)) for various values of t.
<Create slaves of master object>= (<-U)
(setf (slot-value 'slaves)
(list
(send slave-proto :new self)
(send slave-proto :new self :survival-plot t)))
(send (second (slot-value 'slaves)) :title
(strcat "Plot of E(S(t))-" (slot-value 'identifier)))
The sliders for controlling the hyperparameters need to be set up.
There are three dialog items that correspond to every hyperparameter:
a text-item where the hyperparameter name will be displayed, a
value-text-item where the value of the hyperparameter will be
displayed and underneath the first two, a slider showing the slider
stop. Figure [->] shows a typical slider. For nice looks,
the width of the hyperparameter name string and the value string
should add up to the width of the slider.
A typical slider in a dialog
[*]
<Set up dialog and wait for user input>= (<-U)
(let* ((hyperparameter-labels (send self :hyperparameter-names))
(hyperparameter-ranges (send self :hyperparameter-ranges))
(triples <Make triples for sliders>)
(slider-and-dialog-items <Make sliders with triples>)
(dialog-items (select slider-and-dialog-items 0))
(sliders-alone (select slider-and-dialog-items 1))
(reset-button (send button-item-proto :new "Reset"
:action #'(lambda() (send self :reset))))
(timing-button (send button-item-proto :new "Timing:OFF"
:action #'(lambda() (send self :toggle-timing))))
(stats-button (send button-item-proto :new "Statistics"
:action #'(lambda()
(send self :print-all-statistics)))))
(setf (slot-value 'hyperparameter-sliders) sliders-alone)
(setf (slot-value 'timing-button) timing-button)
(call-next-method (list (list reset-button stats-button timing-button)
(list dialog-items))))
Let us first create the triples we need for the hyperparameter sliders
for use with the function make-sliders. The triples are (a) the
label for the hyperparameter, (b) the range within which the
hyperparameter will be varied, and (c) the action that is to be taken
when the slider is pressed as a function.
<Make triples for sliders>= (<-U)
(mapcar
#'(lambda(x y z)
(list x y #'(lambda(w)
(send self :current-hyperparameter-values z w))))
hyperparameter-labels hyperparameter-ranges
(iseq (send self :number-of-hyperparameters)))
We are now ready to create the sliders and associated text and value
items. The function make-sliders returns a multiple-values list
of slider-items using triples with a specified
layout along with a list of the appropriate scroll-items. A
slider-item is a list of two elements, the first element being a
list containing a text-item and a value-item and the second is
a scroll-item.
<Make sliders with triples>= (<-U)
(multiple-value-list
(make-sliders triples
:layout (list (list t) (list t t) (list t))
:formats *default-hyperparameter-print-format*
:no-of-slider-stops (slot-value 'number-of-slider-stops)))
Finally, we create a decent title, become alert once again as opposed
to being lazy and then synchronize everything and show those windows.
<Some final touches>= (<-U)
(send self :title (strcat (send self :identifier) "-Master"))
(let ((hd (send self :number-of-hyperparameters))
(hyper-sliders (slot-value 'hyperparameter-sliders))
(hyper-vals (send self :current-hyperparameter-values)))
(setf hyper-vals (coerce hyper-vals 'list))
(send (select hyper-sliders 1) :display-value)
(dotimes (l hd)
(unless (= l 1)
(send (select hyper-sliders l) :value (select hyper-vals l)))))
(setf (slot-value 'lazy) nil)
(send self :synchronize)
The :process-run-file Method
[*]
In handling the ``run'' file, note the assumptions made about the file
structure and how some irrelevant headings are expected and skipped.
<The Master :process-run-file Method>= (<-U)
(defmeth master-proto :process-run-file (&key number-of-points?)
"Method args: (&key number-of-points?)
Processes a run file to get all inputs."
(format t "~%Processing Run file ... ")
(let ((fh (open (strcat (slot-value 'identifier) ".run") :direction :input))
(file-names nil)
(alpha-shape-scale nil)
(constants nil)
(m nil)
(n nil)
(mn nil)
(index 0)
(indicator-counts nil)
(summary-data nil)
(log-hmix nil))
<Read number of Markov chains>
<Read number of points>
(setf mn (* m n))
<Read number of sets>
<Read number of months>
(setf (slot-value 'importance-weights)
(make-array mn :initial-element 0.0 :element-type 'c-double))
<Process the table containing info on Markov chains>
<Process Hyperparameter Ranges etc.>
<Read summary data, log mixture density, indicator counts>
(setf (slot-value 'indicator-counts) indicator-counts)
(setf (slot-value 'summary-data) summary-data)
(setf (slot-value 'log-mixture-density) log-hmix)
(close fh))
(format t "done~%"))
<Read number of Markov chains>= (<-U)
;; Ignore the information in first line---just some info.
(read fh nil)
;; Ignore the heading "Number of Markov chains"
(read fh nil)
(setf m (read fh nil))
(setf (slot-value 'number-of-markov-chains) m)
<Read number of points>= (<-U)
;; Ignore the heading "Number of Points"
(read fh nil)
(setf n (read fh nil))
(if number-of-points?
(setf n number-of-points?))
(setf (slot-value 'number-of-points) n)
<Read number of sets>= (<-U)
;; Ignore the heading "Number of Data Values"
(read fh nil)
(setf (slot-value 'number-of-data-values) (read fh nil))
<Read number of months>= (<-U)
;; Ignore the heading "Number of Months"
(read fh nil)
(setf (slot-value 'number-of-months) (read fh nil))
<Process Hyperparameter Ranges etc.>= (<-U)
;; Ignore the leading line
(read fh nil)
(let ((range nil)
(r (iseq 3))
(s (iseq 3))
(i (iseq 3)))
(dotimes (j 3)
(read fh nil) (read fh nil)
(setf range (cons (read fh nil) nil))
(read fh nil)
(setf range (cons (read fh nil) range))
(setf (select r j) (reverse range))
(read fh nil)
(setf (select i j) (read fh nil))
(read fh nil)
(setf (select s j) (read fh nil)))
(send self :hyperparameter-ranges
(append r
(list (list 0 (1- (slot-value 'number-of-months))))))
(setf (slot-value 'number-of-slider-stops)
(append s (list (slot-value 'number-of-months))))
(setf (slot-value 'initially-specified-hyperparameter-values)
(append i (list (* 0.5 (slot-value 'number-of-months))))))
<Process the table containing info on Markov chains>= (<-U)
;; Now read the table of file name, alpha,
;; shape, scale and log constant after ignoring the heading.
(read fh nil)
(dotimes (i m)
(setf file-names (cons (read fh nil) file-names))
(setf alpha-shape-scale (cons (read fh nil) alpha-shape-scale))
(setf alpha-shape-scale (cons (read fh nil) alpha-shape-scale))
(setf alpha-shape-scale (cons (read fh nil) alpha-shape-scale))
(setf constants (cons (read fh nil) constants)))
(setf (slot-value 'data-file-names) (reverse file-names))
(setf (slot-value 'log-constants-of-proportionality)
(if (= m 1)
(make-array 1 :initial-element 0
:element-type 'c-double)
(make-array (1- m) :initial-contents (rest (reverse constants))
:element-type 'c-double)))
(setf (slot-value 'hyperparameters-used-in-markov-chains)
(make-array (list m 3)
:initial-contents (reverse alpha-shape-scale)
:element-type 'c-double))
<Read summary data, log mixture density, indicator counts>= (<-U)
;; Ignore heading.
(read fh nil)
(setf indicator-counts
(make-array (list mn (slot-value 'number-of-months))
:element-type 'c-long))
(setf summary-data
(make-array (list mn 2) :initial-element 0.0
:element-type 'c-double))
(setf log-hmix
(make-array mn :initial-element 0.0 :element-type 'c-double))
(dotimes (j n)
(dotimes (k m)
(setf index (+ (* n k) j))
(setf (aref summary-data index 0) (read fh nil))
(setf (aref summary-data index 1) (read fh nil))
(setf (aref log-hmix index) (read fh nil))
(dotimes (l (slot-value 'number-of-months))
(setf (aref indicator-counts index l) (truncate (read fh nil))))))
Defines :process-run-file (links are to index).
The :process-frequency-table Method
[*]
<The Master :process-frequency-table Method>= (<-U)
(defmeth master-proto :process-frequency-table (fh nd ind)
"Args: (fh nd ind)
This is a convenience method used in reading the data files. It reads
the frequency table of size nd from fh, the file stream, and fills
the indicator-counts slot at index ind."
(let ((x-vals (make-array nd))
(freq (make-array nd))
(sort-index nil)
(indicator-counts (slot-value 'indicator-counts))
(number-of-months (slot-value 'number-of-months))
(j 0)
(k 0)
(freq-sum 0)
(smallest-month 0))
<Read frequency table and sort values>
(dotimes (i nd)
(setf smallest-month
<Find smallest month not less than i-th ordered x-val>)
(when smallest-month
<Fill table entries with current frequency sum>
<Bump starting index for next stretch>
<Update frequency sum>))
<Handle situation when x-values run out before month values>))
Defines :process-frequency-table (links are to index).
<Read frequency table and sort values>= (<-U)
(dotimes (l nd)
(setf (aref x-vals l) (read fh nil))
(setf (aref freq l) (read fh nil)))
(setf sort-index (make-array nd :initial-contents (order x-vals)))
In the code snippet below, nil is returned if no such month is
found.
<Find smallest month not less than i-th ordered x-val>= (<-U)
(loop
(if (>= j number-of-months)
(return nil))
(if (>= j (elt x-vals (elt sort-index i)))
(return j))
(setf j (1+ j)))
<Fill table entries with current frequency sum>= (<-U)
(dotimes (x (- smallest-month k))
(setf (aref indicator-counts ind (+ x k)) freq-sum))
<Bump starting index for next stretch>= (<-U)
(setf k smallest-month)
<Update frequency sum>= (<-U)
(setf freq-sum (+ freq-sum (elt freq (elt sort-index i))))
<Handle situation when x-values run out before month values>= (<-U)
(dotimes (x (- number-of-months k))
(setf (aref indicator-counts ind (+ x k)) freq-sum))
The :graphical-interface Method
[*]
<The Master :graphical-interface Method>= (<-U)
(defmeth master-proto :graphical-interface
(&key identifier?
number-of-markov-chains?
number-of-points?
file-names?
log-constants?)
"Method args: (&key identifier?
number-of-chains?
number-of-points
file-names?
log-constants?)
Sets up a graphical interface for some inputs if not specified."
(let* ((id (if identifier?
identifier?
<Get an identifier>))
(m (if number-of-markov-chains?
number-of-markov-chains?
<Get the number of Markov chains>))
(file-names (if file-names?
file-names?
<Get data file names>))
(n (if number-of-points?
number-of-points?
<Get the number of points>))
(initial-guess (if log-constants?
log-constants?
(if (> m 1)
<Get initial guess>))))
(setf (slot-value 'identifier) id)
(setf (slot-value 'data-file-names) file-names)
(setf (slot-value 'number-of-markov-chains) m)
(setf (slot-value 'number-of-points) n)
(setf (slot-value 'number-of-months) *default-number-of-months*)
(setf (slot-value 'indicator-counts)
(make-array (list (* m n) (slot-value 'number-of-months))
:element-type 'c-long))
(setf (slot-value 'importance-weights)
(make-array (* m n) :initial-element 0 :element-type 'c-double))
(setf (slot-value 'log-mixture-density)
(make-array (* m n) :initial-element 0
:element-type 'c-double))
(setf (slot-value 'log-constants-of-proportionality)
(if (= m 1)
(make-array 1
:initial-element 0
:element-type 'c-double)
(make-array (1- m)
:initial-contents initial-guess
:element-type 'c-double)))
(let* ((alphas nil)
(shapes nil)
(scales nil)
(summary-data (make-array (list (* m n) 2) :initial-element 0
:element-type 'c-double )))
<Read in data files and set up data>)))
*
<Get an identifier>= (<-U)
(get-nonempty-string-dialog
(format nil "Please enter a unique short descriptive name~%~
for this exploration.~%~
Ex: CancerData") :initial "BreastCancer")
<Get the number of Markov chains>= (<-U)
(select (get-tested-value-dialog
(format nil "How many Markov chains for ~a?" id)
:initial *default-number-of-markov-chains*
:test #'(lambda(x) (and (numberp x) (> x 0)))
:error-message "No. of Markov chains must be >= 1!")
0)
<Get data file names>= (<-U)
(let* ((dialog-items ())
(filenames nil)
(prompt-item
(send text-item-proto :new
(format nil "Please enter names of all data files without~%~
the extensions for run ~a and click OK.~%" id)))
(col-1 (send text-item-proto :new "MC number"))
(col-2 (send text-item-proto :new " Data Filename "))
(headings (list col-1 col-2))
(width-a (send col-1 :width))
(width-b (send col-2 :width))
(widths (map-elements #'send headings :width))
(file-name-items (map-elements
#'send edit-text-item-proto :new
(repeat "" m)
:text-length 30))
(mc-numbers (mapcar #'(lambda(x)
(send text-item-proto :new
(format nil "~5d" x)))
(1+ (iseq m))))
(abort (send modal-button-proto :new "Abort" :action #'top-level))
(dialog nil)
(ok (send modal-button-proto :new "OK"
:action
#'(lambda()
(map-elements #'send file-name-items :text)))))
(dolist (x mc-numbers)
(send x :width width-a))
(dolist (x file-name-items)
(send x :width width-b))
(setf dialog-items
(list prompt-item
headings
(list (mapcar #'(lambda (x y) (list x y))
mc-numbers file-name-items))
(list ok abort)))
(loop
(setf dialog (send modal-dialog-proto :new dialog-items))
(setf file-names (send dialog :modal-dialog))
(if (or (some-files-dont-exist (map-elements #'strcat file-names ".in"))
(some-files-dont-exist (map-elements #'strcat file-names ".out")))
(message-dialog "Some files don't exist. Please try again!")
(return file-names))))
<Get the number of points>= (<-U)
(select (get-tested-value-dialog
(format nil "How many points to use for ~a in reweighting?" id)
:initial *default-number-of-points*
:test #'(lambda(x) (and (numberp x) (> x 0)))
:error-message "No. of points must be > 0!")
0)
<Get initial guess>= (<-U)
(select (get-nonnil-value-dialog
(format nil "Enter an initial guess for estimating the ~% ~
constants in the format shown below.~% ~
Dimension should be ~d!~%" (1- m))
:initial (1+ (iseq (1- m)))
:test #'(lambda(x)
(let ((val (select x 0)))
(and val (listp val) (= (length val (1- m))))))
:error-message "Improper guess!")
0)
To process the data, we use both the input and output file used in
running the Markov chains.
<Read in data files and set up data>= (<-U) [D->]
(dotimes (j m)
;;; First process the input file parameters
(let ((fh (open (strcat (select file-names j) ".in") :direction :input)))
;; discard the first four values
(dotimes (i 4)
(read fh nil))
(setf shapes (cons (read fh nil) shapes))
(setf scales (cons (read fh nil) scales))
(setf alphas (cons (read fh nil) alphas))
;; discard next three values (warmup, iterations and gap)
;; Assumes all input files are consistent.
(read fh nil) (read fh nil) (read fh nil)
(setf (slot-value 'number-of-data-values) (read fh nil))
(close fh))
;;; Now process the corresponding output file.
(let ((fh (open (strcat (select file-names j) ".out") :direction :input))
(ind (* j n))
(nd nil))
(dotimes (i n)
(setf (aref summary-data ind 1) (read fh nil))
(setf nd (read fh nil))
(setf (aref summary-data ind 0) nd)
(send self :process-frequency-table fh nd ind)
(setf ind (1+ ind)))
(close fh)))
Now that all the Markov chain files have been processed, the
slot-values can be set up.
<Read in data files and set up data>+= (<-U) [<-D]
(setf (slot-value 'hyperparameters-used-in-markov-chains)
(make-array (list m 3)
:initial-contents (bind-columns (reverse alphas)
(reverse shapes)
(reverse scales))
:element-type 'c-double))
(setf (slot-value 'summary-data) summary-data)
The :create-run-file Method
[*]
<The Master :create-run-file Method>= (<-U)
(defmeth master-proto :create-run-file ()
"Method args: ()
Creates a data and a lisp file for subsequent runs."
(format t "~%Creating Run file for subsequent runs ... ")
(let ((mn nil)
(index 0)
(m (slot-value 'number-of-markov-chains))
(n (slot-value 'number-of-points))
(p (slot-value 'number-of-data-values))
(file-names (slot-value 'data-file-names))
(summary-data (slot-value 'summary-data))
(indicator-counts (slot-value 'indicator-counts))
(hypers-used (slot-value 'hyperparameters-used-in-markov-chains))
(log-constants (slot-value 'log-constants-of-proportionality))
(log-hmix (slot-value 'log-mixture-density))
(fh (open (strcat (slot-value 'identifier) ".run")
:direction :output)))
(setf mn (* m n))
<Write number of Markov chains, etc.>
<Write table of info on Markov chains>
<Write info on hyperparameters>
<Write summary data, log mixture density, indicator counts>
(close fh))
<Create lisp file for subsequent runs>
(format t "done~%")
<Show informative message>)
*
<Write number of Markov chains, etc.>= (<-U)
(format fh "~s~%~%" "Automatically generated file. Do not edit unless~
you know what you are doing!")
(format fh "~s ~g~%" "Number of Markov Chains" m)
(format fh "~s ~g~%" "Number of Points" n)
(format fh "~s ~g~%" "Number of Data Values"
(slot-value 'number-of-data-values))
(format fh "~s ~g~%" "Number of Months used"
(slot-value 'number-of-months))
<Write info on hyperparameters>= (<-U)
(let ((r (send self :hyperparameter-ranges))
(s (slot-value 'number-of-slider-stops))
(i (slot-value 'initially-specified-hyperparameter-values)))
(format fh "~%~s~%"
"Hyperparameter Exploration Range, Initial Value, stops, etc.")
(format fh "~%~s ~s ~g ~s ~g ~s ~g ~s ~g~%"
"M(R)" "Min" (select (select r 0) 0)
"Max" (select (select r 0) 1)
"Initial" (select i 0)
"Stops" (select s 0))
(format fh "~%~s ~s ~g ~s ~g ~s ~g ~s ~g~%"
"Theta-Shape" "Min" (select (select r 1) 0)
"Max" (select (select r 1) 1)
"Initial" (select i 1)
"Stops" (select s 1))
(format fh "~%~s ~s ~g ~s ~g ~s ~g ~s ~g~%"
"Theta-Scale" "Min" (select (select r 2) 0)
"Max" (select (select r 2) 1)
"Initial" (select i 2)
"Stops" (select s 2)))
<Write table of info on Markov chains>= (<-U)
(format fh "~%~%~s~%~%"
"Table of file name, alpha, shape, scale, log constants")
(dotimes (j m)
(format fh "~s ~g ~g ~g ~g~%"
(select file-names j)
(aref hypers-used j 0)
(aref hypers-used j 1)
(aref hypers-used j 2)
(if (= j 0) -1 (elt log-constants (1- j)))))
<Write summary data, log mixture density, indicator counts>= (<-U)
(format fh "~%~%~s~%~%"
"Table of no of distinct values, theta, log-hmix, Delta-X-t.")
(dotimes (j n)
(dotimes (k m)
(setf index (+ (* k n) j))
(format fh "~g ~g " (aref summary-data index 0)
(aref summary-data index 1))
(format fh "~g " (aref log-hmix index))
(dotimes (l (slot-value 'number-of-months))
(format fh "~d " (aref indicator-counts index l)))
(format fh "~%")))
<Create lisp file for subsequent runs>= (<-U)
(let* ((id (slot-value 'identifier))
(fh (open (strcat id ".lsp") :direction :output)))
(format fh ";;;Automatically generated file. Do not edit~
unless you know what you are doing.~%")
(format fh "(require ~s)~%" "bsa")
(format fh "(use-package ~s)~%" "BSA")
(format fh "(defvar ~a-master (send master-proto :new
:identifier? ~s))~%" id id)
(close fh))
<Show informative message>= (<-U)
(let ((id (slot-value 'identifier)))
(message-dialog
(format nil "For your information: Two files~%~
~a.run and ~a.lsp~%were created.~%~
To repeat this run quickly the next time around~%~
you only need to load the file ~a.lsp into~%~
xlispstat." id id id)))
The :synchronize Method
[*]
The :synchronize method is responsible for synchronizing the
slaves so that the density estimates they display is for the current
values of the hyperparameters. Thus, if any hyperparameter value is
changed, the :synchronize method should be invoked.
<The Master :synchronize Method>= (<-U)
(defmeth master-proto :synchronize ()
"Method args: ()
Synchronizes all slaves."
(when (not (slot-value 'lazy))
(send self :consolidate-computation)
(let* ((t-format (select *default-hyperparameter-print-format* 3))
(slave1 (first (slot-value 'slaves)))
(slave2 (second (slot-value 'slaves)))
(t-value (send self :current-hyperparameter-values 3)))
(send slave1
:title (format nil "Law of F(~v,vf)-~a"
(first t-format)
(second t-format) t-value
(slot-value 'identifier)))
(setf (select (slot-value 'statistics-labels) 0)
(format nil "E(S(~g))" t-value))
(setf (select (slot-value 'statistics-labels) 1)
(format nil "SD(S(~g))" t-value))
(send slave1 :start-buffering)
(send slave1 :clear-lines :draw nil)
(send slave1 :add-lines
(slot-value 'density-abscissae)
(slot-value 'density-ordinates))
(send slave1 :adjust-to-data)
(send slave1 :buffer-to-screen)
(send slave2 :start-buffering)
(send slave2 :clear-lines :draw nil)
(send slave2 :add-lines
(slot-value 'expectation-abscissae)
(slot-value 'expectation-ordinates))
(send slave2 :adjust-to-data)
(send slave2 :buffer-to-screen))))
Defines :synchronize (links are to index).
The :consolidate-computation Method
[*]
In this method we compute all the relevant statistics such as
E(S(t)) and sigma_S(t) for various values of t, as well as
the effective sample size. This method was added to consolidate all
computations for acceptable performance and makes several other
methods above superfluous.
<The Master :consolidate-computation Method>= (<-U)
(defmeth master-proto :consolidate-computation ()
"Method args: ()
Consolidates all computations for sake of speed. This method
recalculates the importance weights, computes relevant statistics and
saves them in the slots."
(call-by-reference-oldcfun "consolidateComputation"
(slot-value 'shared-library) *default-number-of-plot-stops*
(slot-value 'statistics)
(slot-value 'density-abscissae) (slot-value 'density-ordinates)
(slot-value 'expectation-abscissae)
(slot-value 'expectation-ordinates)
(slot-value 'standard-deviation-ordinates)))
Defines :consolidate-computation (links are to index).
The :reset Method
[*]
<The Master :reset Method>= (<-U)
(defmeth master-proto :reset ()
"Method args: ()
Resets the state of all objects."
(send self :current-hyperparameter-values
(send self :initially-specified-hyperparameter-values)))
Defines :reset (links are to index).
The :effective-sample-size Method
[*]
The actual computation of the effective sample size is now done in the
C code that computes the importance weights (see
section [->]). We note that that the formula for
computing effective sample size was originally used in estimating
standard deviations in stratified sampling. It can be used to give
ballpark figures, but should not be taken too literally.
<The Master :effective-sample-size Method>= (<-U)
(defmeth master-proto :effective-sample-size ()
(select (slot-value 'statistics) 2))
Defines :effective-sample-size (links are to index).
The :print-all-statistics Method
[*]
<The Master :print-all-statistics Method>= (<-U)
(defmeth master-proto :print-all-statistics ()
(let* ((hyperparameter-names (slot-value 'hyperparameter-names))
(hyper-strings (send self :labelled-hyperparameter-values))
(max-name-len (max (mapcar #'length hyperparameter-names))))
(format t "~%~%*** Statistics for ~a ***~%~%" (send self :identifier))
(format t "Hyperparameter Settings:~%" )
(dolist (x hyper-strings)
(format t " ~a~%" x))
(format t "Total Sample Size = ~5d.~%"
(* (slot-value 'number-of-markov-chains)
(slot-value 'number-of-points)))
(format t "Effective Sample Size = ~5d.~%"
(select (slot-value 'statistics) 2)))
(dolist (x (slot-value 'slaves))
(send x :print-summary))
(format t "~%*** End of Statistics ***~%~%"))
Defines :print-all-statistics (links are to index).
The :labelled-hyperparameter-values Method
[*]
<The Master :labelled-hyperparameter-values Method>= (<-U)
(defmeth master-proto :labelled-hyperparameter-values ()
(let* ((hyper-names (slot-value 'hyperparameter-names))
(max-name-len (max (mapcar #'length hyper-names))))
(mapcar
#'(lambda(a b)
(let ((y (select *default-hyperparameter-print-format* a)))
(if (listp y)
(format nil "~va = ~v,vf" max-name-len
b (first y) (second y)
(send self :current-hyperparameter-values a))
(format nil (strcat "~va = " y) max-name-len b
(send self :current-hyperparameter-values a)))))
(iseq (send self :number-of-hyperparameters)) hyper-names)))
Defines :labelled-hyperparameter-values (links are to index).
The :statistics Method
[*]
The :statistics method for the master returns the value of
the slot statistics.
<The Master :statistics Method>= (<-U)
(defmeth master-proto :statistics ()
"Method args: ()
Returns the slot value statistics."
(slot-value 'statistics))
Defines :statistics (links are to index).
The :statistics-print-formats Method
[*]
The :statistics-print-formats method for the master returns
the value of the slot
statistics-print-formats.
<The Master :statistics-print-formats Method>= (<-U)
(defmeth master-proto :statistics-print-formats ()
"Method args: ()
Returns the slot value statistics-print-formats."
(slot-value 'statistics-print-formats))
Defines :statistics-print-formats (links are to index).
The :statistics-labels Method
[*]
The :statistics-labels method for the master returns the
value of the slot statistics-labels.
<The Master :statistics-labels Method>= (<-U)
(defmeth master-proto :statistics-labels ()
"Method args: ()
Returns the slot value statistics-labels."
(slot-value 'statistics-labels))
Defines :statistics-labels (links are to index).
The :toggle-timing Method
[*]
This method toggles timing on and off.
<The Master :toggle-timing Method>= (<-U)
(defmeth master-proto :toggle-timing ()
"Method args: ()
Toggles timing on and off."
(send self :hide-window)
(setf (slot-value 'timing) (not (slot-value 'timing)))
(if (slot-value 'timing)
(send (slot-value 'timing-button)
:slot-value 'text "Timing:ON")
(send (slot-value 'timing-button) :slot-value 'text
"Timing:OFF"))
(send self :show-window))
Defines :toggle-timing (links are to index).
The :superimpose Method
[*]
This method returns the value of the superimpose slot.
<The Master :superimpose Method>= (<-U)
(defmeth master-proto :superimpose ()
"Method args: ()
Returns the value of slot superimpose."
(slot-value 'superimpose))
Defines :superimpose (links are to index).
The :close Method
[*]
The :close method for the master must close the slaves window
that are active. Finally, it must commit hara-kiri.
<The Master :close Method>= (<-U)
(defmeth master-proto :close ()
"Method args: ()
Kills all subordinate slave and commits suicide."
(dolist (x (slot-value 'slaves))
(send x :remove))
(call-next-method))
Defines :close (links are to index).
[*]
We shall use the prefix ``BSA'' for Bayesian Sensitivity Analysis.
<Defaults for Master>= (<-U)
(defvar *default-master-object-prefix* "BSA-")
(defvar *default-number-of-markov-chains* 8)
(defvar *default-number-of-points* 50)
(defvar *default-number-of-plot-stops* 51)
(defvar *default-no-of-slider-stops* '(64 51 121 61))
(defvar *default-hyperparameter-print-format*
'((7 2) (7 2) (7 2) (7 2)))
(defvar *default-number-of-months* 61)
(defvar *default-maximization-tolerance* 1e-5)
(defvar *default-shared-lib-name*
(if (member ':msdos xlisp::*features*)
"win/bsa.dll"
"./libbsa@SHLIB_SUFFIX@"))
Defines *default-hyperparameter-print-format*, *default-master-object-prefix*, *default-maximization-tolerance*, *default-no-of-slider-stops*, *default-number-of-markov-chains*, *default-number-of-months*, *default-number-of-points*, *default-shared-lib-name (links are to index).
[*]
The Slave prototype inherits from scatterplot-proto.
<The Slave Prototype>= (<-U)
(defproto slave-proto '(master survival-plot color-index)
() scatterplot-proto
"The Slave prototype. Master is its master upon whom the slave
relies for all data. Survival-plot is nil or t indicating what is
to be plotted. Color-index is used in superimposition.")
<Methods for Slave Prototype>
<Defaults for Slave>
Defines slave-proto (links are to index).
The methods for slave-proto follow.
<Methods for Slave Prototype>= (<-U)
<The Slave :isnew Method>
<The Slave :redraw-background Method>
<The Slave :redraw-statistics Method>
<The Slave :print-summary Method>
<The Slave :close Method>
The :isnew Method
[*]
<The Slave :isnew Method>= (<-U)
(defmeth slave-proto :isnew (master &key (go-away t) (survival-plot nil))
"Method args: master &rest args
Creates a new instance of the slave-proto object."
(setf (slot-value 'master) master)
(setf (slot-value 'survival-plot) survival-plot)
(setf (slot-value 'color-index) 0)
(call-next-method 2 :go-away go-away :draw nil)
(if (not (slot-value 'survival-plot))
(let ((skip (+ (send self :text-ascent) (send self :text-descent)))
(len (length (send master :statistics-print-formats))))
(send self :margin
0 (+ (* len skip) (send self :text-descent)) 0 0
:draw nil)))
(send self :redraw))
Defines :isnew (links are to index).
The :redraw-background Method
[*]
<The Slave :redraw-background Method>= (<-U)
(defmeth slave-proto :redraw-background ()
"Method args: ()
Redraws the background of the screen"
(send self :start-buffering)
(call-next-method)
(unless (slot-value 'survival-plot)
(send self :redraw-statistics))
(send self :buffer-to-screen))
Defines :redraw-background (links are to index).
The :redraw-statistics Method
[*]
<The Slave :redraw-statistics Method>= (<-U)
(defmeth slave-proto :redraw-statistics ()
"Method args: ()
Redraws the statistics on the screen"
(let* ((master (slot-value 'master))
(ascent (send self :text-ascent))
(descent (send self :text-descent))
(skip (+ ascent descent))
(em (send self :text-width "m"))
(en (send self :text-width "n"))
(canvas-width (send self :canvas-width))
(y 0)
(stats-labels (send master :statistics-labels))
(name-width (max (mapcar #'(lambda(x) (send self :text-width x))
stats-labels)))
(value-strings (mapcar #'(lambda(x y)
(format nil x y))
(send master :statistics-print-formats)
(coerce (send master :statistics) 'list)))
(value-width (max (mapcar #'(lambda(x)
(send self :text-width x))
value-strings)))
(x2 (- canvas-width value-width em))
(x1 (- x2 name-width em)))
(dotimes (i (length stats-labels))
(let ((sl (select stats-labels i))
(sv (select value-strings i)))
(setf y (+ y skip))
(send self :draw-text sl x1 y 0 0)
(send self :draw-text sv x2 y 0 0)))))
Defines :redraw-statistics (links are to index).
The :print-summary Method
[*]
This method basically calculates the mean and variance of beta
distributions.
<The Slave :print-summary Method>= (<-U)
(defmeth slave-proto :print-summary ()
"Method args: ()
Prints the mean and variance of the Law of F(t)."
(if (slot-value 'survival-plot)
(let* ((master (slot-value 'master))
(x (send master :slot-value 'expectation-abscissae))
(y (send master :slot-value 'expectation-ordinates))
(stddevs (send master :slot-value
'standard-deviation-ordinates)))
(format t "~%Table of E(S(t))~%")
(format t "------------------------------------~%")
(format t "Time t (months) E(S(t)) SD(S(t))~%")
(format t "------------------------------------~%")
(dotimes (i (length x))
(format t "~14d ~5,3f ~5,3f~%" (aref x i)
(aref y i) (aref stddevs i))))
(let* ((master (slot-value 'master))
(stats (send master :slot-value 'statistics))
(stat-print-formats (send master :slot-value
'statistics-print-formats))
(hypers (send master :current-hyperparameter-values))
(alpha (aref hypers 0))
(t-value-index (truncate (aref hypers 3))))
(format t (strcat "Mean of S(~d): "
(first stat-print-formats) "~%")
t-value-index (elt stats 0))
(format t (strcat "Std. Dev. of S(~d): "
(second stat-print-formats) "~%")
t-value-index (elt stats 1)))))
Defines :print-summary (links are to index).
The :close Method
[*]
<The Slave :close Method>= (<-U)
(defmeth slave-proto :close ()
(ok-or-cancel-dialog "Please use the master to quit"))
[*]
<Defaults for Slave>= (<-U)
(defvar *default-slave-plot-size* '(300 265))
(defvar *default-slave-plot-stops* 51)
(defvar *default-statistics-labels*
'("E(S(t))" "SD(S(t))" "Eff. Sample Size"))
(defvar *default-statistics-print-formats*
'("~5,3f" "~5,3f" "~5,0f"))
Defines *default-slave-plot-size*, *default-slave-plot-stops*, *default-statistics-labels*, *default-statistics-print-formats* (links are to index).
The C Programs
[*]
In designing the C programs, we assume one thing---that no compaction
is done during the lifetime of the master object. This allows us to
refrain from passing pointers to data arrays each time a routine is
called. Instead, all required pointers could be passed once when the
master object is created and they remain fixed for the life of the
master. Current versions of Lisp-Stat don't do memory
compaction or move objects around.
Here's our copyright for C programs.
<C Copyright>= (U->)
/**
*** $Revision: 1.34 $ of $Date: 1998/07/02 17:16:33 $
***
*** Copyright (C) 1994, 1995, 1998. Doss and Narasimhan
***
*** Hani J. Doss (doss@stat.ohio-state.edu) and
*** B. Narasimhan (naras@stat.stanford.edu)
***
*** This program is free software; you can redistribute it and/or modify
*** it under the terms of the GNU General Public License as published by
*** the Free Software Foundation; either version 2 of the License, or
*** (at your option) any later version.
***
*** This program is distributed in the hope that it will be useful,
*** but WITHOUT ANY WARRANTY; without even the implied warranty of
*** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
*** GNU General Public License for more details.
***
*** You should have received a copy of the GNU General Public License
*** along with this program; if not, write to the Free Software
*** Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
**/
Our C routines follow.
<C Routines>=
<C Copyright>
<Header files>
<Beta functions>
<Global Variables>
<Address Initialization Routine>
<The Density Routine>
<The LogProbability Routine>
<The LogLikelihood Routine>
<The hthetaOverHmix Routine>
<The compute-log-hmix Routine>
<The calcWeights Routine>
<The computeStatistics Method>
<The computeLawOfFOfT Method>
<The computeMeanOfFBarOfT Method>
<The consolidateComputation Routine>
Here are the headers files we will need. Note the declaration for the
prototype of the gamma function that comes with
Lisp-Stat. This function is not declared static in the
Lisp-Stat source and is therefore available to us.
<Header files>= (<-U)
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
[*]
We need two other functions, betdens and logbeta are indeed
declared static in the Lisp-Stat source. Therefore, a
copy of these routines is needed. In Lisp-Stat, the lnGamma
function is known as mygamma.
<Beta functions>= (<-U)
extern double mygamma(double a);
static double logbeta (double a, double b)
{
static double da = 0.0, db = 0.0, lbeta = 0.0;
if (a != da || b != db) { /* cache most recent call */
da = a; db = b;
lbeta = mygamma(da) + mygamma(db) - mygamma(da + db);
}
return(lbeta);
}
static double betadens (double x, double a, double b)
{
if (x <= 0.0 || x >= 1.0)
return 0.0;
return exp(log(x) * (a - 1) + log(1 - x) * (b - 1) - logbeta(a, b));
}
Defines betadens, logbeta, mygamma (links are to index).
[*]
We begin with some global variables for all these C routines. Note
that these are declared static.
<Global Variables>= (<-U)
static double *summaryData, *hyperValuesUsedInMarkovChains,
*currentHyperValues, *importanceWeights,
*logMixtureDensity,
*logConstantsOfProportionality;
static double *statistics;
static int numberOfMarkovChains, numberOfPoints,
numberOfPointsForConstants,
numberOfDataValues, numberOfMonths, numberOfDataValues,
*indicatorCounts;
#define LOW_LOG_PROBABILITY -1e10;
Defines currentHyperValues, hyperValuesUsedinMarkovChains, importanceWeights, indicatorCounts, logConstantsOfProportionality, logMixtureDensity, LOW_LOG_PROBABILITY, numberOfDataValues, numberOfMarkovChains, numberOfMonths, numberOfPoints, numberOfPointsForConstants, statistics, summaryData (links are to index).
[*]
This initialization routine gathers addresses and stores them a static
storage area. On the Windows platform, we also need to export the
routines defined here.
<Address Initialization Routine>= (<-U)
#ifdef _Windows
__declspec(dllexport) void initializeAddress (double *a, double *b,
double *c, double *d, double *e, double *ff, int *g, int *h,
int *i, int *j, int *k, double *l, int *m);
#endif
/*
* Initialization routine, must be called as soon as the slots
* are created in the master :isnew method.
*/
void initializeAddress (argSummaryData, argHyperValuesUsedinMarkovChains,
argCurrentHyperValues, argImportanceWeights,
argLogMixtureDensity, argLogConstantsOfProportionality,
argNumberOfMarkovChains, argNumberOfPoints,
argNumberOfMonths, argNumberOfDataValues,
argIndicatorCounts, argStatistics, argNumberOfPointsForConstants)
double *argSummaryData, *argHyperValuesUsedinMarkovChains,
*argCurrentHyperValues, *argImportanceWeights,
*argLogMixtureDensity, *argLogConstantsOfProportionality;
double *argStatistics;
int *argNumberOfMarkovChains, *argNumberOfPoints,
*argNumberOfMonths, *argNumberOfDataValues, *argIndicatorCounts,
*argNumberOfPointsForConstants;
{
summaryData = argSummaryData;
hyperValuesUsedInMarkovChains = argHyperValuesUsedinMarkovChains;
currentHyperValues = argCurrentHyperValues;
importanceWeights = argImportanceWeights;
logMixtureDensity = argLogMixtureDensity;
logConstantsOfProportionality = argLogConstantsOfProportionality;
indicatorCounts = argIndicatorCounts;
numberOfMarkovChains = *argNumberOfMarkovChains;
numberOfPoints = *argNumberOfPoints;
numberOfMonths = *argNumberOfMonths;
numberOfDataValues = *argNumberOfDataValues;
statistics = argStatistics;
numberOfPointsForConstants = *argNumberOfPointsForConstants;
}
Defines initializeAddress (links are to index).
[*]
This routine calculates the j-th posterior density at a point. See
:f method for master-proto. The formulas used for this and the
next three methods all follow the development in [cite geye:1994].
<The Density Routine>= (<-U)
double f (argJ, argEta, argIndex)
int argJ, argIndex;
double *argEta;
{
double nj, numberOfDistinctXS = summaryData[2 * argIndex],
theta = summaryData[2 * argIndex + 1], alpha, shape, scale;
int index = argJ * 3;
if (argJ == 0)
nj = -1.0;
else
nj = argEta[argJ - 1];
alpha = hyperValuesUsedInMarkovChains[index];
shape = hyperValuesUsedInMarkovChains[index + 1];
scale = hyperValuesUsedInMarkovChains[index + 2];
return exp(numberOfDistinctXS * log(alpha) +
shape * log(scale) + (shape - 1.0) * log(theta) -
scale * theta - mygamma(shape) + nj);
}
Defines f (links are to index).
[*]
This routine calculates the probability that a given point comes
from a particular Markov chain. See :logp method for
master-proto.
<The LogProbability Routine>= (<-U)
double logp (argJ, argEta, argIndex)
int argJ, argIndex;
double *argEta;
{
double tmp, tmp1, sum = 0;
int k;
for (k = 0; k < numberOfMarkovChains; k++) {
tmp = f(k, argEta, argIndex);
if (k == argJ)
tmp1 = tmp;
sum += tmp;
}
if (tmp1 > 0.0)
return log(tmp1 / sum);
else
return LOW_LOG_PROBABILITY;
}
Defines logp (links are to index).
[*]
The routine that calculates the log-likelihood. See :loglik method
for master-proto.
<The LogLikelihood Routine>= (<-U)
#ifdef _Windows
__declspec(dllexport) void logLikelihood (double *a, double *b);
#endif
void logLikelihood (argEta, result)
double *argEta, *result;
{
double sum = 0.0;
int i, j, index = 0;
for (j = 0; j < numberOfMarkovChains; j++) {
for (i = 0; i < numberOfPointsForConstants; i++) {
sum += logp(j, argEta, index);
index++;
}
}
*result = sum;
}
Defines loglik (links are to index).
[*]
See :htheta-over-hmix method for master-proto.
<The hthetaOverHmix Routine>= (<-U)
double hthetaOverHmix (argIndex)
int argIndex;
{
double numberOfDistinctXS, theta, alpha, shape, scale;
numberOfDistinctXS = summaryData[2 * argIndex];
theta = summaryData[2 * argIndex + 1];
alpha = currentHyperValues[0];
shape = currentHyperValues[1];
scale = currentHyperValues[2];
return exp(numberOfDistinctXS * log(alpha) +
shape * log(scale) + (shape - 1.0) * log(theta) -
scale * theta - mygamma(shape)
- logMixtureDensity[argIndex]);
}
Defines hthetaOverHmix (links are to index).
[*]
See :compute-log-hmix method for master-proto.
<The compute-log-hmix Routine>= (<-U)
#ifdef _Windows
__declspec(dllexport) void computeLogHmix ();
#endif
void computeLogHmix ()
{
int i, k, mn;
double sum;
mn = numberOfMarkovChains * numberOfPoints;
for (i = 0; i < mn; i++) {
sum = 0.0;
for (k = 0; k < numberOfMarkovChains; k++) {
sum += f(k, logConstantsOfProportionality, i);
}
logMixtureDensity[i] = log(sum);
}
}
Defines computeLogHmix (links are to index).
[*]
See :calc-weights method for master-proto.
<The calcWeights Routine>= (<-U)
#ifdef _Windows
__declspec(dllexport) void calcWeights ();
#endif
void calcWeights ()
{
int i, j, index = 0, sampleSize;
double tmp, sum = 0.0, meanWeight, sumCenteredWeightsSquared;
for (j = 0; j < numberOfMarkovChains; j++) {
for (i = 0; i < numberOfPoints; i++) {
tmp = hthetaOverHmix(index);
importanceWeights[index] = tmp;
sum += tmp;
index++;
}
}
index = 0;
sampleSize = numberOfMarkovChains * numberOfPoints;
meanWeight = 1.0 / sampleSize;
sumCenteredWeightsSquared = 0.0;
for (j = 0; j < numberOfMarkovChains; j++) {
for (i = 0; i < numberOfPoints; i++) {
importanceWeights[index] /= sum;
tmp = importanceWeights[index] - meanWeight;
sumCenteredWeightsSquared += tmp * tmp;
index++;
}
}
tmp = sumCenteredWeightsSquared / (sampleSize - 1);
statistics[2] = floor(sampleSize /
(1 + sampleSize * sampleSize * tmp));
}
Defines calcWeights (links are to index).
[*]
See :compute-statistics method for slave-proto. In computing
the mean and variance of the mixture of Beta densities, we have to use
the facts E(X) = E(E(X|W)) and V(X) = E(V(X|W)) + V(E(X|W)). In
our case, F(t) is a Beta random variable X_i with probability
w_i, the normalized importance weight.
<The computeStatistics Method>= (<-U)
#ifdef _Windows
__declspec(dllexport) void computeStatistics (double *a);
#endif
void computeStatistics (argStatistics)
double *argStatistics;
{
double alpha = currentHyperValues[0],
expectation = 0.0,
expectationOfConditionalSquaredMean = 0.0,
expectationOfVariance = 0.0,
tmp,
sumCenteredWeightsSquared = 0,
betaMean, betaVariance,
betaParam1, betaParam2;
int i, j, index = 0, tValueIndex = (int) currentHyperValues[3];
<Compute Mean and Variance in C>
argStatistics[0] = 1.0 - expectation;
argStatistics[1] = sqrt(expectationOfConditionalSquaredMean -
expectation * expectation);
}
Defines computeStatistics (links are to index).
<Compute Mean and Variance in C>= (<-U U->)
for (j = 0; j < numberOfMarkovChains; j ++) {
for (i = 0; i < numberOfPoints; i++) {
betaParam1 = <C Expression for first beta parameter>;
betaParam2 = <C Expression for second beta parameter>;
betaMean = betaParam1 / (betaParam1 + betaParam2);
tmp = importanceWeights[index] * betaMean;
betaVariance = betaMean * betaParam2 /
((betaParam1 + betaParam2) * (betaParam1 + betaParam2 + 1));
expectation += tmp;
expectationOfVariance += importanceWeights[index] * betaVariance;
expectationOfConditionalSquaredMean += tmp * betaMean;
index++;
}
}
The C Equivalent of :compute-law-of-f-of-t Method
[*]
See :compute-law-of-f-of-t method for slave-proto. Note the
difference in that some parameters need to be passed in this C routine
unlike the Lisp routine.
This method basically calculates for equally spaced values of u in [0,1] the weighted sum
where 
denotes the beta density with parameters a
and b evaluated at u.
<The computeLawOfFOfT Method>= (<-U)
#ifdef _Windows
__declspec(dllexport) void computeLawOfFOfT (int *n, double *a, double *b);
#endif
void computeLawOfFOfT (int *n, double *x, double *y)
{
double alpha, sum, betaParam1, betaParam2;
int i, j, k, index, tValueIndex;
alpha = currentHyperValues[0];
tValueIndex = (int) currentHyperValues[3];
for (k = 0; k < *n; k++) {
sum = 0.0;
index = 0;
for (j = 0; j < numberOfMarkovChains; j++) {
for (i = 0; i < numberOfPoints; i++) {
betaParam1 = <C Expression for first beta parameter>;
betaParam2 = <C Expression for second beta parameter>;
sum += importanceWeights[index] *
betadens(x[k], betaParam1, betaParam2);
index++;
}
}
y[k] = sum;
}
}
Defines computeLawOfFOfT (links are to index).
Here are the formulas for the two beta density parameters. The first
one is of course 
with the understanding that 
indicates degree of concentration around the
exponential family for theta.
The expression for the first beta parameter follows.
<C Expression for first beta parameter>= (<-U <-U)
(indicatorCounts[index * numberOfMonths + tValueIndex] +
alpha * (1 - exp(-tValueIndex * summaryData[2*index + 1])))
The second one is
which in C is given by the snippet below.
<C Expression for second beta parameter>= (<-U <-U)
(alpha + numberOfDataValues - betaParam1)
[*]
See :compute-mean-of-fbar-of-t method for slave-proto. This
method basically calculates for several values of t, the quantity
E(S(t)). The values of t for which this is done is the same as the
list of values used in the slider for t.
<The computeMeanOfFBarOfT Method>= (<-U)
#ifdef _Windows
__declspec(dllexport) void computeMeanOfFBarOfT (int *n, double *a,
double *b, double *c);
#endif
void computeMeanOfFBarOfT (n, x, y, standardDeviations)
double *x, *y, *standardDeviations;
int *n;
{
double alpha = currentHyperValues[0],
expectation = 0.0,
expectationOfConditionalSquaredMean = 0.0,
expectationOfVariance = 0.0,
tmp,
sumCenteredWeightsSquared = 0,
betaMean, betaVariance,
betaParam1, betaParam2;
int i, j, k, index = 0, tValueIndex;
y[0] = 1.0;
standardDeviations[0] = 0.0;
for (k = 1; k < *n; k++) {
index = 0;
expectation = 0.0;
expectationOfConditionalSquaredMean = 0.0;
expectationOfVariance = 0.0;
tValueIndex = (int) (x[k]);
<Compute Mean and Variance in C>
y[k] = 1 - expectation;
standardDeviations[k] = sqrt(expectationOfConditionalSquaredMean -
expectation * expectation);
}
}
Defines computeMeanOfFBarOfT (links are to index).
[*]
This method consolidates all computation by first computing importance
weights, then computing all the statistics in one shot.
<The consolidateComputation Routine>= (<-U)
#ifdef _Windows
__declspec(dllexport) void consolidateComputation (int *n,
double *a, double *b, double *c, double *d, double *e,
double *ff);
#endif
void consolidateComputation (n, argStatistics, argDensityAbscissae,
argDensityOrdinates,
argExpectationAbscissae, argExpectationOrdinates,
argStandardDeviationOrdinates)
double *argStatistics, *argDensityAbscissae, *argDensityOrdinates,
*argExpectationAbscissae, *argExpectationOrdinates,
*argStandardDeviationOrdinates;
int *n;
{
int k;
calcWeights();
computeLawOfFOfT(n, argDensityAbscissae, argDensityOrdinates);
computeMeanOfFBarOfT(&numberOfMonths,
argExpectationAbscissae, argExpectationOrdinates,
argStandardDeviationOrdinates);
k = (int) currentHyperValues[3];
statistics[0] = argExpectationOrdinates[k];
statistics[1] = argStandardDeviationOrdinates[k];
}
Defines consolidateComputation (links are to index).
[*]
We present some information on how to compile and use the software on
various platforms. This is presented as a README file.
<README file>=
Copyright (C) Hani J. Doss and B. Narasimhan
--------------------------------------------
This file provides information on installing and using
BSA (Bayesian Sensitivity Analysis) software. We expect any
serious user of the software to read our paper that has now
appeared in
@InBook{doss:nara:1998b,
author = {Hani J. Doss and B. Narasimhan},
title = {Practical Nonparametric and Semiparametric
Bayesian Statistics},
chapter = {Dynamic Display of Changing Posterior in Bayesian
Survival Analysis},
publisher = {Springer},
year = {1998},
editor = {Dipak Dey, Peter M\"uller, Debajyoti Sinha},
number = {133},
series = {Lecture Notes in Statistics},
pages = {63--87}
}
We have included a copy of this paper (ddg.ps.gz) in the subdirectory ddg.
Requirements:
A) XLisp-Stat 3.53-5 or higher. Freely available from
ftp://ftp.stat.umn.edu/pub/xlispstat/current
B) Windows or Unix. Windows includes 95 and NT.
(Mac version is in development.)
C) On Unix, you will also need a C compiler.
Step 1.
------
You probably received the entire package as a compressed archive
named bsa.tgz. On Unix, the contents of the archive may be
extracted into a directory called bsa by executing the commands:
gunzip -c bsa.tgz | tar -xvf -
If you are using GNU tar, this can be done in one shot via:
tar -xvzf bsa.tgz
To extract the files on Windows, you need to use an extractor
like WinZip. See http://www.winzip.com for more details.
Step 2
------
On Windows, you skip this step. Macs are not supported yet because
the authors don't know how to create a dynamic shared library.
On Unix, you need to configure the software to your environment.
Change to the bsa directory and type
./configure
make
This will compile the lisp files and create a shared library for your
platform.
Step 3
------
Start using the program. In Unix, the command
xlispstat BreastCancerRadiationOnly
or
xlispstat BreastCancerRadiationChemo
will fire up the readymade examples.
On Windows, you fire up Lisp-Stat and load the files
BreastCancerRadiationOnly.lsp
or
BreastCancerRadiationChemo.lsp
to proceed.
If you are inclined to use commands instead of mouse-clicks, you can
send ``messages'' to the master object in the listener window as
shown in the following examples:
(send BreastCancerRadiationOnly-master :current-hyperparameter-values
'(1 4 150 12))
(send BreastCancerRadiationOnly-master :print-all-statistics)
(send BreastCancerRadiationOnly-master :slot-value
'bsa::importance-weights)
Note how the master-objects names are related to the lisp file names and
how internal slot-names have to be prefixed by the package name.
For dealing with a new problem, we provide a few points regarding the
software. A number of inputs are required for running the
program. These are discussed in detail in the literate program
(bsa.ps) under the section titled ``Introduction.'' For convenience we
repeat the details here. This excerpted part is indented two spaces
for easy reference.
First note that the software only does sensitivity analysis.
No general facility is provided for generating observations from
Markov chains. Indeed, since the range of models for which MCMC
methods are applicable are large and such methods most likely involve
problem-specific issues, it is our opinion that building such a
supertool, if it is at all possible, is a non-trivial task. However,
the Fortran program used in generating the output for our example is
included along with this software and can be used for models similar
to ours. Of course, any appropriate method may be used to generate
the samples as long as the output is available in a form usable by our
software. The requirements on the data that can be used with our
software are spelt out below.
Corresponding to each Markov chain output, there must be two files
with the extensions ".in" (input file) and ".out" (output file).
For example, "mc1.in" and "mc1.out".
The input file must have the following structure. The first four items
in the file can be anything, string or number, either on a single line
or any conceivable combination of lines. The next three items
*must* be the shape of the Gamma distribution on theta, the
scale of the Gamma distribution on theta---the parametrization for
shape a and scale b is proportional to x^{a-1} exp(-bx)---and
M(R). The next three values values following these quantities can
be anything, but the one following it should be the number of data points,
that is, the number of sets or intervals. In the Fortran program we use
-99 is used to denote infinity. Nothing else is read from the input file.
The output file must have the following structure for each data point
generated by the Markov chain. The value of theta must be followed
by the number of distinct values of the data points, which must be
followed by a frequency table of the actual data value and the
corresponding frequency. The layout of the values on lines does not
matter as long as at least a single white space delimits values. If
this structure is violated, errors will result. A peek at the data
files included with this software will help the reader.
It is assumed that a proper installation of XLisp-Stat is available.
For a new problem, you probably have several Markov chain output files
although even one should work. (In the latter case, reweighting reduces
to simple Importance Sampling.)
a) It is best to create a new directory for your problem and have
your data files there. For example, the directory "BreastCancer"
contains relevant data files for our Breast Cancer data.
b) The only files you actually need to run the program are:
1) Either one of bsa.fsl or bsa.lsp
2) Either one of utility.fsl or utility.lsp
3) Either one of call-by-reference.fsl or call-by-reference.lsp
and
4) the shared library libbsa.so or libbsa.sl as the case
may be. On Windows, instead of the shared library,
we need the whole subdirectory "win".
Copy these files/directories to where you have the data
files and work there.
5) The file new-problem.lsp.
Invoke Lisp-Stat and load the file named "new-problem.lsp".
The first time (and first time only), the following inputs will be
needed.
Inputs
------
1) An indentifier for uniquely identifying the run. Use a meaningful
name here. Let us assume this is BreastCancer (the default)
in the discussion below.
2) The number of Markov chain outputs that you want to use for
reweighting. Must be >= 1, with 1 denoting straight Importance
Sampling.
3) The names of the files containing output from Markov chains,
*without the extensions*. The software will automatically
tag on the extensions .in and .out when looking for files.
4) The number of points per chain to use in the dynamic reweighting.
Thus if you specify 50 and have 8 chains, then 50 points from each
of the eight chains (= 400) will be used.
5) An initial guess for maximizing the log-quasilikelihood which
will provide an estimate of the constants of proportionality.
6) The range between which you want to vary the hyperparameters. If
you use only one chain, then you *must* specify the range. Otherwise,
the range will be a single point. If you specified many chains,
the default settings for each hyperparameter will be the minimum and
maximum values from values used in all Markov chains. The number of
stops should be an odd number if you want to hit the middle of the
interval.
7) The number of points to use in estimating the constants of
proportionality. If you use all of the data, the estimation can
take a while. It is almost always better to go with the default
or less. (If you are really interested in using more points,
then start off with 10, and use the estimates thus obtained
to start your larger optimization. This will save you a lot of
time.)
Once you specify this, the maximization will take place. This is a
good point to go refill your coffee cup.
After the estimation, two files are created so that you are not
bombarded with questions in subsequent explorations. For example,
BreastCancer.lsp (and)
BreastCancer.run
For repeating the exploration next time, you only need to load
the file BreastCancer.lsp into XLisp-Stat. This will bypass
all the inputs we discussed above except for the question about
ranges. The file BreastCancer.run contains pre-processed
information for faster loading and will be used when
BreastCancer.lsp is loaded.
All files are text files and can be viewed with a text viewer.
Examples for Breast Cancer Data
-------------------------------
The two files
BreastCancerRadiationOnly.lsp (and)
BreastCancerRadiationChemo.lsp
and the corresponding run files are provided for experimentation.
These exist in the main directory "bsa" itself and concern the
dataset on two treatments described in the paper
Dynamic Display of Changing Posterior in
Bayesian Survival Analysis
by
Hani J. Doss and B. Narasimhan
By default, they use 50 points each and 8 Markov chains. We wish to
note that an earlier version of the software was used to produce the
results in the paper and subsequently a bug was found. This does not
change any of the conclusions of the paper but the numbers shown in
table 1.1 in the paper are off from the actual values obtained using
the software. A replacement is provided in the file newtable.tex
and shows that the agreement between estimates obtained by reweighting
and those obtained by actual runs of Markov chains is, if anything,
better that what the original table indicated.
Just load the lisp files into XLisp-Stat to do the dynamic exploration.
If you have a sufficiently fast machine, you can use more points.
To completely reproduce our work from scratch, you need to use the
data files in the subdirectory "BreastCancer".
The authors may be contacted via email at:
doss@stat.ohio-state.edu (Hani J. Doss)
naras@stat.Stanford.EDU (B. Narasimhan)
Note on the program itself
--------------------------
The programs are written in a literate style using the Noweb literate
programming tools. We provide two utility packages that one can use
independently of the program: utility.lsp and call-by reference.lsp.
The former contains functions we have found useful in writing Lispstat
programs; the latter implements a call-by-reference glue between
Lispstat and C.
Enjoy!
Improvements needed
[*]
Acknowledgement
[*]
This research was supported by Air Force Office of Scientific Research
Grant F49620-94-1-0028.
We are indebted to Fred W. Huffer for the Fortran program dirichlet
that comes with this software. All we did was modify it slightly to
suit our data.
Index of Code Chunks
[*]
This list is generated automatically. The numeral is that of the first
definition of the chunk.
- <*>: D1
- <Address Initialization Routine>: U1, D2
- <Beta functions>: U1, D2
- <Bump starting index for next stretch>: U1, D2
- <C Copyright>: D1, U2
- <C Expression for first beta parameter>: U1, U2, D3
- <C Expression for second beta parameter>: U1, U2, D3
- <C Routines>: D1
- <Compute Mean and Variance in C>: U1, D2, U3
- <Copyright>: D1, U2
- <Create lisp file for subsequent runs>: U1, D2
- <Create slaves of master object>: U1, D2
- <Defaults for Master>: U1, D2
- <Defaults for Slave>: U1, D2
- <Fill table entries with current frequency sum>: U1, D2
- <Find smallest month not less than i-th ordered x-val>: U1, D2
- <Get an identifier>: U1, D2
- <Get data file names>: U1, D2
- <Get initial guess>: U1, D2
- <Get the number of Markov chains>: U1, D2
- <Get the number of points>: U1, D2
- <Global Variables>: U1, D2
- <Handle situation when x-values run out before month values>: U1, D2
- <Header files>: U1, D2
- <Make sliders with triples>: U1, D2
- <Make triples for sliders>: U1, D2
- <Methods for Master Prototype>: U1, D2
- <Methods for Slave Prototype>: U1, D2
- <Pass data array addresses to C Routines>: U1, D2
- <Perform Reverse Logistic Regression>: U1, D2
- <Process Hyperparameter Ranges etc.>: U1, D2
- <Process the table containing info on Markov chains>: U1, D2
- <Read frequency table and sort values>: U1, D2
- <Read in data files and set up data>: U1, D2, D3
- <Read number of Markov chains>: U1, D2
- <Read number of months>: U1, D2
- <Read number of points>: U1, D2
- <Read number of sets>: U1, D2
- <Read summary data, log mixture density, indicator counts>: U1, D2
- <README file>: D1
- <Set up dialog and wait for user input>: U1, D2
- <Set up Hyperparameter-ranges and stops>: U1, D2
- <Set up slot values for master object>: U1, D2, D3, D4, D5, D6, D7, D8
- <Show informative message>: U1, D2
- <Some final touches>: U1, D2
- <The calcWeights Routine>: U1, D2
- <The computeLawOfFOfT Method>: U1, D2
- <The compute-log-hmix Routine>: U1, D2
- <The computeMeanOfFBarOfT Method>: U1, D2
- <The computeStatistics Method>: U1, D2
- <The consolidateComputation Routine>: U1, D2
- <The Density Routine>: U1, D2
- <The hthetaOverHmix Routine>: U1, D2
- <The LogLikelihood Routine>: U1, D2
- <The LogProbability Routine>: U1, D2
- <The Master :calc-weights Method>: U1, D2
- <The Master :close Method>: U1, D2
- <The Master :compute-log-hmix Method>: U1, D2
- <The Master :consolidate-computation Method>: U1, D2
- <The Master :create-run-file Method>: U1, D2
- <The Master :current-hyperparameter-values Method>: U1, D2
- <The Master :effective-sample-size Method>: U1, D2
- <The Master :graphical-interface Method>: U1, D2
- <The Master :hyperparameter-names Method>: U1, D2
- <The Master :hyperparameter-ranges Method>: U1, D2
- <The Master :hyperparameter-sliders Method>: U1, D2
- <The Master :hyperparameters-used-in-markov-chains Method>: U1, D2
- <The Master :identifier Method>: U1, D2
- <The Master :importance-weights Method>: U1, D2
- <The Master :indicator-counts Method>: U1, D2
- <The Master :initially-specified-hyperparameter-values Method>: U1, D2
- <The Master :isnew Method>: U1, D2
- <The Master :labelled-hyperparameter-values Method>: U1, D2
- <The Master :log-constants-of-proportionality Method>: U1, D2
- <The Master :loglik Method>: U1, D2
- <The Master :log-mixture-density Method>: U1, D2
- <The Master :number-of-data-values Method>: U1, D2
- <The Master :number-of-hyperparameters Method>: U1, D2
- <The Master :number-of-markov-chains Method>: U1, D2
- <The Master :number-of-months Method>: U1, D2
- <The Master :number-of-points Method>: U1, D2
- <The Master :print-all-statistics Method>: U1, D2
- <The Master :process-frequency-table Method>: U1, D2
- <The Master :process-run-file Method>: U1, D2
- <The Master Prototype>: U1, D2
- <The Master :reset Method>: U1, D2
- <The Master :slaves Method>: U1, D2
- <The Master :statistics Method>: U1, D2
- <The Master :statistics-labels Method>: U1, D2
- <The Master :statistics-print-formats Method>: U1, D2
- <The Master :summary-data Method>: U1, D2
- <The Master :superimpose Method>: U1, D2
- <The Master :synchronize Method>: U1, D2
- <The Master :toggle-timing Method>: U1, D2
- <The Slave :close Method>: U1, D2
- <The Slave :isnew Method>: U1, D2
- <The Slave :print-summary Method>: U1, D2
- <The Slave Prototype>: U1, D2
- <The Slave :redraw-background Method>: U1, D2
- <The Slave :redraw-statistics Method>: U1, D2
- <Update frequency sum>: U1, D2
- <Write info on hyperparameters>: U1, D2
- <Write number of Markov chains, etc.>: U1, D2
- <Write summary data, log mixture density, indicator counts>: U1, D2
- <Write table of info on Markov chains>: U1, D2
Index of Identifiers
[*]
Here is a list of the identifiers used, and where they appear.
Underlined entries indicate the place of definition.
This index is generated automatically.
- betadens: D1, U2
- :calc-weights: D1, U2
- calcWeights: U1, D2, U3
- :close: D1, U2
- computeLawOfFOfT: D1, U2
- compute-log-hmix: D1, U2
- computeLogHmix: U1, D2
- computeMeanOfFBarOfT: D1, U2
- computeStatistics: D1
- :consolidate-computation: U1, D2
- consolidateComputation: U1, D2
- copyright: D1
- :current-hyperparameter-values: D1, U2, U3, U4, U5, U6, U7, U8
- current-hyperparameter-values: D1, U2, U3, U4, U5, U6, U7, U8, U9
- currentHyperValues: D1, U2, U3, U4, U5, U6, U7
- data-file-names: D1
- *default-hyperparameter-print-format*: U1, U2, U3, D4
- *default-master-object-prefix*: D1
- *default-maximization-tolerance*: U1, D2
- *default-no-of-slider-stops*: D1
- *default-number-of-markov-chains*: U1, D2
- *default-number-of-months*: U1, D2
- *default-number-of-points*: U1, U2, D3
- *default-shared-lib-name: U1, D2
- *default-slave-plot-size*: D1
- *default-slave-plot-stops*: D1
- *default-statistics-labels*: U1, D2
- *default-statistics-print-formats*: U1, D2
- density-abscissae: D1
- density-ordinates: D1
- :effective-sample-size: D1
- expectation-abscissae: D1
- expectation-ordinates: D1
- f: D1, U2, U3
- hthetaOverHmix: D1, U2
- :hyperparameter-names: D1, U2
- hyperparameter-names: D1, U2, U3, U4
- :hyperparameter-ranges: D1, U2, U3, U4, U5
- hyperparameter-ranges: D1, U2, U3, U4, U5, U6, U7
- :hyperparameter-sliders: D1
- hyperparameter-sliders: D1, U2
- :hyperparameters-used-in-markov-chains: D1
- hyperparameters-used-in-markov-chains: D1, U2
- hyperValuesUsedinMarkovChains: D1
- :identifier: D1, U2, U3, U4, U5
- identifier: D1, U2, U3, U4, U5, U6, U7, U8, U9, U10, U11
- :importance-weights: D1
- importance-weights: D1, U2, U3
- importanceWeights: D1, U2, U3, U4, U5
- :indicator-counts: D1
- indicator-counts: D1, U2, U3, U4, U5, U6, U7, U8, U9
- indicatorCounts: D1, U2, U3
- initializeAddress: U1, D2
- :initially-specified-hyperparameter-values: D1, U2, U3, U4
- initially-specified-hyperparameter-values: D1, U2, U3, U4, U5
- :isnew: D1, D2, U3
- :labelled-hyperparameter-values: U1, D2
- lazy: D1
- logbeta: D1
- :log-constants-of-proportionality: D1
- log-constants-of-proportionality: D1, U2
- logConstantsOfProportionality: D1, U2, U3
- :loglik: D1, U2
- loglik: U1, U2, D3
- :log-mixture-density: D1
- log-mixture-density: D1, U2
- logMixtureDensity: D1, U2, U3, U4
- logp: D1, U2
- LOW_LOG_PROBABILITY: D1, U2
- master-proto: U1, D2, U3, U4, U5, U6, U7, U8, U9, U10, U11, U12, U13, U14, U15, U16, U17, U18, U19, U20, U21, U22, U23, U24, U25, U26, U27, U28, U29, U30, U31, U32, U33, U34, U35, U36, U37, U38, U39, U40, U41
- mygamma: D1, U2, U3
- :number-of-data-values: D1
- number-of-data-values: D1, U2
- numberOfDataValues: D1, U2, U3
- :number-of-hyperparameters: U1, D2, U3, U4, U5
- :number-of-markov-chains: D1, U2
- number-of-markov-chains: D1, U2, U3, U4, U5, U6, U7, U8, U9
- numberOfMarkovChains: D1, U2, U3, U4, U5, U6, U7, U8
- :number-of-months: D1, U2
- number-of-months: D1, U2, U3, U4, U5, U6, U7, U8, U9
- numberOfMonths: D1, U2, U3, U4
- :number-of-points: D1, U2
- number-of-points: D1, U2, U3, U4, U5, U6, U7, U8, U9, U10, U11
- numberOfPoints: D1, U2, U3, U4, U5, U6
- numberOfPointsForConstants: D1, U2, U3
- number-of-slider-stops: D1
- :print-all-statistics: U1, D2, U3
- :print-summary: U1, D2
- :process-frequency-table: D1, U2
- :process-run-file: U1, U2, D3
- :redraw-background: D1
- :redraw-statistics: U1, D2
- :reset: U1, D2
- shared-library: D1
- slave-proto: U1, D2, U3, U4, U5, U6, U7
- :slaves: D1
- slaves: D1, U2, U3
- standard-deviation-ordinates: D1
- :statistics: D1, U2, U3, U4, U5
- statistics: D1, U2, U3, U4, U5, U6, U7, U8, U9, U10, U11, U12, D13, U14, U15, U16, U17
- :statistics-labels: D1, U2
- statistics-labels: D1, U2, U3, U4, U5
- :statistics-print-formats: D1, U2, U3
- statistics-print-formats: D1, U2, U3, U4, U5, U6
- :summary-data: D1
- summary-data: D1, U2, U3, U4, U5, U6, U7, U8, U9
- summaryData: D1, U2, U3, U4, U5
- :superimpose: D1
- superimpose: D1, U2
- :synchronize: U1, U2, D3
- timing: D1, U2, U3
- :toggle-timing: U1, D2
- work-space: D1
References
[1]
Hani J. Doss and B. Narasimhan.
Bayesian Poisson regression: Sensitivity analysis through dynamic
graphics.
Technical report, Penn State Erie, The Behrend College, 1994.
[2]
Hani J. Doss and B. Narasimhan.
Dynamic display of changing posterior in Bayesian survival
analysis.
In Practical Nonparametric and Semiparametric Bayesian
Statistics, D. Dey, P. Müller, and D. Sinha, eds, 1998.
Springer-Verlag, N. Y.
[3]
C. J. Geyer.
Estimating normalizing constants and reweighting mixtures in Markov
chain Monte Carlo.
Technical report, School of Statistics, University of Minnesota,
1994.
[4]
Donald E. Knuth.
Literate Programming.
Center for Study of Language and Information, Stanford University,
1992.
[5]
A. Kong, J. S. Liu, and W. H. Wong.
Sequential imputations and Bayesian missing data problems.
J. Amer. Statist. Assoc., 1994.
[6]
Christopher Lee.
Literate programming---propaganda and tools.
World Wide Web, 1994.
Available from
http://www.ius.cs.cmu.edu/help/Programming/literate.html.
[7]
Norman Ramsey.
Literate programming simplified.
IEEE Software, 1994.
[8]
Norman Ramsey.
The Noweb Hacker's Guide, 1994.
Manual for Noweb.
[9]
Norman Ramsey.
The noweb home page.
World Wide Web, 1995.
On the Web at http://www.cs.virginia.edu/ nr/noweb/.
[10]
Luke Tierney.
LISP-STAT: An Object-Oriented Environment for Statistical
Computing and Dynamic Graphics.
John Wiley &Sons, 1990.