[r-cran-coda] 01/60: Imported Upstream version 0.9-1
Andreas Tille
tille at debian.org
Fri Dec 16 12:11:21 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-coda.
commit 00b5a369723f43c8a661286ee3d8698c44bcecd8
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 12:09:53 2016 +0100
Imported Upstream version 0.9-1
---
AUTHORS | 37 ++
CHANGELOG | 261 ++++++++++++
COPYING | 340 +++++++++++++++
DESCRIPTION | 11 +
INDEX | 44 ++
NAMESPACE | 95 +++++
R/codamenu.R | 963 +++++++++++++++++++++++++++++++++++++++++++
R/cumuplot.R | 40 ++
R/gelman.R | 236 +++++++++++
R/geweke.R | 81 ++++
R/heidel.R | 200 +++++++++
R/jags.R | 22 +
R/mcextractor.R | 94 +++++
R/mcmc.R | 310 ++++++++++++++
R/mcmclist.R | 191 +++++++++
R/output.R | 386 +++++++++++++++++
R/raftery.R | 97 +++++
R/thin.R | 0
R/util.R | 255 ++++++++++++
README | 37 ++
data/line.R | 230 +++++++++++
man/Cramer.Rd | 28 ++
man/as.ts.mcmc.Rd | 24 ++
man/autocorr.Rd | 30 ++
man/autocorr.plot.Rd | 24 ++
man/bugs2jags.Rd | 33 ++
man/coda.options.Rd | 77 ++++
man/codamenu.Rd | 15 +
man/crosscorr.Rd | 25 ++
man/crosscorr.plot.Rd | 27 ++
man/cumuplot.Rd | 30 ++
man/densplot.Rd | 37 ++
man/effectiveSize.Rd | 28 ++
man/gelman.diag.Rd | 96 +++++
man/gelman.plot.Rd | 60 +++
man/geweke.diag.Rd | 48 +++
man/geweke.plot.Rd | 47 +++
man/heidel.diag.Rd | 70 ++++
man/linepost.Rd | 16 +
man/mcmc.Rd | 56 +++
man/mcmc.convert.Rd | 46 +++
man/mcmc.list.Rd | 56 +++
man/mcmc.subset.Rd | 33 ++
man/mcmcUpgrade.Rd | 29 ++
man/mcpar.Rd | 26 ++
man/multi.menu.Rd | 37 ++
man/nchain.Rd | 36 ++
man/plot.mcmc.Rd | 32 ++
man/raftery.diag.Rd | 88 ++++
man/read.and.check.Rd | 44 ++
man/read.coda.Rd | 52 +++
man/read.coda.interactive.Rd | 39 ++
man/read.openbugs.Rd | 44 ++
man/spectrum0.Rd | 63 +++
man/spectrum0.ar.Rd | 34 ++
man/summary.mcmc.Rd | 36 ++
man/thin.Rd | 28 ++
man/time.mcmc.Rd | 36 ++
man/traceplot.Rd | 29 ++
man/varnames.Rd | 37 ++
man/window.mcmc.Rd | 36 ++
misc/change.tfoption.Rd | 29 ++
misc/codamenu2.Rd | 40 ++
misc/read.yesno.Rd | 29 ++
misc/set.mfrow.Rd | 26 ++
65 files changed, 5686 insertions(+)
diff --git a/AUTHORS b/AUTHORS
new file mode 100644
index 0000000..8403eee
--- /dev/null
+++ b/AUTHORS
@@ -0,0 +1,37 @@
+The main authors of CODA
+------------------------
+
+Mary Kathryn Cowles University of Nebraska Medical Center, USA
+ <kcowles at stat.uiowa.edu>
+Nicky Best Imperial college, London, UK
+ <n.best at ic.ac.uk>
+Karen Vines Open University, Milton Keynes, UK
+ <S.K.Vines at open.ac.uk>
+Martyn Plummer International Agency for Research on Cancer,
+ Lyon, France
+ <plummer at iarc.fr>
+
+CODA was conceived and motivated by Mary Kathryn Cowles, who developed the
+original program as part of her PhD thesis concerning practical issues
+in implementing the Gibbs sampler. It was modified by Nicky Best to
+provide software suitable for general release, including contributions
+by Karen Vines. Martyn Plummer ported CODA to R and modified the code
+so that the user is not bound by the menu-driven interface.
+
+All authors, and the MRC Biostatistic Unit, have given permission for
+the license terms to be changed to GPL.
+
+Further credits
+---------------
+
+The function "gelman.diag" is based on the "itsim" function
+contributed to the Statlib archive by Andrew Gelman
+<gelman at wald.stat.columbia.edu>.
+
+The function "raftery.diag" is based on the "gibbsit" function
+contributed to the Statlib archive by Steven Lewis
+<slewis at stat.washington.edu>.
+
+Andrew Martin <admartin at wustl.edu> added namespace facilities and
+fixed a large number of bugs that this change uncovered (see CHANGELOG)
+
diff --git a/CHANGELOG b/CHANGELOG
new file mode 100644
index 0000000..0e429e9
--- /dev/null
+++ b/CHANGELOG
@@ -0,0 +1,261 @@
+0.9-1
+- spectrum0 function now has default max.length argument of 200.
+ This means that the output will be batched to a length between
+ 100 and 200 before fitting the glm to the spectrogram. This
+ should improve robustness for chains with high autocorrelation,
+ or markedly non-gaussian distributions.
+- The read.bugs() function has been removed; read.coda() has been
+ modified to allow specification of both output and index files;
+ read.bugs.interactive() has been modified in the same way;
+ read.openbugs() is a new wrapper function around read.coda() for
+ OpenBUGS output.
+
+0.8-3
+- Fixed documentation errors in coda.menu.Rd and linepost.Rd
+
+0.8-2
+- Added generic function thin to list of exported functions
+
+0.8-1
+Continuing problems with namespace:
+- Addition of namespace requires new version number for coda,
+ as saved workspaces are not backwards compatible.
+- Ensured that mcmc attributes are not assumed for objects returned by
+ as.matrix.mcmc, as this now no longer returns an mcmc object.
+ This occurred in, for example, effectiveSize (which returned NA)
+ and gelman.diag (which dropped variable names).
+- Changed "[.mcmc.list" and "[.mcmc" so that they return an
+ mcmc.list/mcmc object respectively when subsetting columns.
+ plot.mcmc() and plot.mcmc.list() now work again when there is
+ only one variable.
+- Imported the required time series generics from package stats.
+ Failure to do this may result in a saved workspace that cannot be reloaded.
+Other problems
+- Removed S compatibility (statements conditional on is.R and wrapper
+ function coda.global.assign). Note that S compatibility never worked
+ at all, and I now have no intention of supporting it.
+- Fixed (old) bug in "[.mcmc" which made column subsets of mcmc objects
+ return invisibly.
+- Changed the plotting functions so it is no longer necessary to press
+ return to see the first page of plots.
+
+0.7-3 (changes done by Andrew Martin <admartin at wustl.edu>)
+- Added NAMESPACE
+- export only functions that I think should have been exported (based on
+ the documentation and common sense).
+- Registered all S3 methods.
+- Fixed documentation for non-exported functions.
+- Patched the "mcmc" function to deal with really big thinning intervals.
+- Patched the "as.matrix.mcmc" function so it really returns matrices.
+- Fixed plot method, by having plot.mcmc pass an mcmc object
+ rather than a matrix, and by fixing the [.mcmc.list method so it
+ returns an mcmc.list not a matrix.
+- Fixed "varnames<-" which was broken when the as.matrix method was
+ fixed.
+- Fixed documentation mcmc.convert.Rd such that the usage is consistent
+ with S3 class definition (this fixed an error thrown in the QC tools).
+
+0.7-2
+- The spectrum0() function now returns zero when it is given a constant
+ vector. summary.mcmc() (which calls spectrum0) now works correctly.
+
+0.7
+- Modified to run on R 1.9.0 with new organization of base library
+
+0.6-2
+- Fixed documentation bug in raftery.diag
+
+0.6-1
+- spectrum0.ar no longer crashes when the chain is a linear function
+ of the iteration number
+- codamenu automatically drops variables that are linear functions
+ of the iteration number
+- read.bugs renamed to read.coda. read.bugs exists as an alias.
+- initial support for JAGS: bugs2jags function converts WinBUGS
+ data to R dump format used by JAGS.
+- added cumuplot function (not yet in codamenu system)
+
+0.5-14
+- Fixed reporting of sample size in densplot(), for variables that are
+ in the range [0,1] or [0,Inf). Thanks to Roy Levy.
+
+0.5-13
+- mcmc() now works with data frames (provided that they contain only
+ numeric values).
+
+0.5-12
+- Documentation errors in coda.options.Rd and nchain.Rd fixed.
+ Thanks to Kurt Hornik.
+
+0.5-11
+- Fixed bug in mcmc() function that causes problems in Geweke diagnostic
+- Fixed geweke.plot so that it works with samples running from iteration
+ N to 2N. These were previously mis-diagnosed as being too short.
+ Thanks to Vasco Leemans for finding both bugs.
+- In autocorr.plot, the `ask' argument was not used. Thanks to
+ Dennis A Wolf.
+
+0.5-10
+- Eliminated use of "=" for assignment operator. This is a syntax error
+ for R < 1.4.0
+
+0.5-9
+- Further documentation bugs removed
+
+0.5-8
+- Removed further documentation bugs found by "R CMD check coda" using
+ R-1.4 (pre-release)
+
+0.5-7
+- Removed obsolete line.doc and line.old.doc files from data directory.
+- New spectrum0.ar provides model-based estimate of spectral density
+ at frequency zero.
+- New effectiveSize diagnostic gives effective sample size.
+- codamenu includes automatic check on effective sample size.
+
+0.5-6
+- Removed further documentation bugs found by "R CMD check coda" using
+ R-1.3 (pre-release)
+
+0.5-5
+- Ironed out last warnings generated by "R CMD check coda"
+
+0.5-4
+
+- Provided documentation for all functions and datasets.
+ (Thanks to Kurt Hornik for the prompting)
+
+0.5-3
+
+- Fixed bug in example for mcmc.list. Row subsetting no longer preserves
+ mcmc objects (Thanks to Kurt Hornik).
+
+0.5-2
+
+- Fixed bug in gelman.transform which did not work for univariate chains
+ (Thanks to Mark A. Beaumont)
+- Fixed confidence limits in geweke.plot
+ (Thanks to Mark A. Beaumont)
+- Allow user to set ylim in densplot (Thanks to Niels Peter Baadsgaard)
+
+0.5-1
+
+- Replaced time series functions with functions from "ts" library
+ (R-base >= 0.65.0)
+- Removed calls to Version() (deprecated). Use is.R() instead.
+- Added new function read.yesno
+- Source files maintained using ESS
+- Allowed restart() in codamenu.options.plot.kernel
+- Simplified print.coda.options
+- Removed "onepage" option in coda.options() (Subsumed in user.layout)
+- Removed "mrows" and "mcols" options in coda.options (use par instead)
+
+Changes to Geweke's diagnostic
+- Uses new function spectrum0() to estimate spectral density at zero
+- Gelman-Rubin-Brooks plot never discards more than half the chain
+ to preserve necessary asymptotic conditions.
+
+Changes to Gelman and Rubin's diagnostic
+- Multivariate psrf added.
+- Documentation for Gelman-Rubin-Brooks plot update to give
+ clearer motivation.
+
+Changes to Heidelberger and Welch's diagnostic
+- Simplified formula for Cramer-von Mises statistic
+- Using new function spectrum0() to estimate spectral density at zero
+- Can set p-value threshold for passing convergence test.
+- p-value is printed in output, using new function pcramer().
+- Prints starting iteration of truncated chain instead of number
+ of iterations to discard.
+
+0.4-7
+
+Fixed bug in read.bugs.interactive() leading to failure when
+user enters both ".ind" and ".out" names (Thanks to John Logsdon).
+
+0.4-6
+
+Archive 0.4-5 was incorrectly compressed with "compress" instead
+of "gzip". Corrected by Friedrich Leisch.
+
+0.4-5
+
+Bug fixes
+- densplot failed with show.obs=TRUE when scale was "positive" or
+ "proportion".
+- as.matrix.mcmc failed to preserve start, end and thin.
+- codamenu did not tidy up on exit.
+
+0.4-4
+
+Started S3 compatibility
+Fixed bug which caused options menus to crash
+Fixed legend bug in gelman.plot
+Confirmed that these bugs are fixed:
+* densplot "missing" and "scale" bugs (Thanks to Greg Warnes)
+* autocorr "improper time series parameters" bug ("acf"
+ function rewritten by Paul Gilbert)
+* integer overflow bug in raftery.diag (Thanks to Morten
+ Frydenberg)
+* read.bugs.interactive will search for the files it needs
+ and print their names.
+codamenu now assigns default variable and chain names to
+data when these are NULL.
+
+0.4-3
+
+Fixed help page errors pointed out by Brian Ripley.
+Fixed coda.credits
+
+0.4-2
+
+Whoops. 0.4-1 was a mistake.
+
+0.4-1
+
+updated manual pages
+new class mcmc.list added to deal with multiple chains. ugrade.mcmc
+ function introduced to deal with old mcmc objects.
+plot functions changed to use the "ask" parameter instead of
+ "pause" functions.
+spec.pgram now handles matrix time-series.
+acf function now calculates cross-correlations. Thanks to Paul Gilbert.
+codamenu functions now use title argument in "menu"
+corrected spelling mistakes in help pages
+changed instances of "T" and "F" to "TRUE" and "FALSE" (R coding standards)
+densplot now recognizes discrete distributions and prints histogram ...
+... also prints histogram if IQR=0 (large mass on one point)
+fixed bug in mcmc which allowed non-integer thinning intervals
+Fixed manual pages with bad use of "alias" command.
+
+0.3-4
+
+Fixed bug in "tspar<-" which breaks much of the code in R-0.62
+
+0.3-3
+
+Package was in obsolete format. Corrected by Fritz.
+
+** Pre-release changes for R version
+
+Created class "mcmc" with associated constructor and extractor
+functions as well as plot, print and summary methods.
+
+Modified diagnostics so they all work on objects of class "mcmc"
+and can be called directly by the user. All diagnostics return
+objects with associated print methods. Renamed some functions
+and arguments for ease of use.
+
+Changed the menu driven interface - now called by the function
+"codamenu" - to avoid recursive calling of menu functions.
+
+Put frequently used code inside utility functions.
+
+Got rid of functions written by Mathsoft
+
+Wrote drop-in replacements for some time series functions which
+are found in S-PLUS but not R.
+
+The logfile facility has been removed. Sorry.
+
+Changed license terms to GPL.
diff --git a/COPYING b/COPYING
new file mode 100644
index 0000000..d60c31a
--- /dev/null
+++ b/COPYING
@@ -0,0 +1,340 @@
+ GNU GENERAL PUBLIC LICENSE
+ Version 2, June 1991
+
+ Copyright (C) 1989, 1991 Free Software Foundation, Inc.
+ 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+ Preamble
+
+ The licenses for most software are designed to take away your
+freedom to share and change it. By contrast, the GNU General Public
+License is intended to guarantee your freedom to share and change free
+software--to make sure the software is free for all its users. This
+General Public License applies to most of the Free Software
+Foundation's software and to any other program whose authors commit to
+using it. (Some other Free Software Foundation software is covered by
+the GNU Library General Public License instead.) You can apply it to
+your programs, too.
+
+ When we speak of free software, we are referring to freedom, not
+price. Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+this service if you wish), that you receive source code or can get it
+if you want it, that you can change the software or use pieces of it
+in new free programs; and that you know you can do these things.
+
+ To protect your rights, we need to make restrictions that forbid
+anyone to deny you these rights or to ask you to surrender the rights.
+These restrictions translate to certain responsibilities for you if you
+distribute copies of the software, or if you modify it.
+
+ For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must give the recipients all the rights that
+you have. You must make sure that they, too, receive or can get the
+source code. And you must show them these terms so they know their
+rights.
+
+ We protect your rights with two steps: (1) copyright the software, and
+(2) offer you this license which gives you legal permission to copy,
+distribute and/or modify the software.
+
+ Also, for each author's protection and ours, we want to make certain
+that everyone understands that there is no warranty for this free
+software. If the software is modified by someone else and passed on, we
+want its recipients to know that what they have is not the original, so
+that any problems introduced by others will not reflect on the original
+authors' reputations.
+
+ Finally, any free program is threatened constantly by software
+patents. We wish to avoid the danger that redistributors of a free
+program will individually obtain patent licenses, in effect making the
+program proprietary. To prevent this, we have made it clear that any
+patent must be licensed for everyone's free use or not licensed at all.
+
+ The precise terms and conditions for copying, distribution and
+modification follow.
+
+ GNU GENERAL PUBLIC LICENSE
+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+ 0. This License applies to any program or other work which contains
+a notice placed by the copyright holder saying it may be distributed
+under the terms of this General Public License. The "Program", below,
+refers to any such program or work, and a "work based on the Program"
+means either the Program or any derivative work under copyright law:
+that is to say, a work containing the Program or a portion of it,
+either verbatim or with modifications and/or translated into another
+language. (Hereinafter, translation is included without limitation in
+the term "modification".) Each licensee is addressed as "you".
+
+Activities other than copying, distribution and modification are not
+covered by this License; they are outside its scope. The act of
+running the Program is not restricted, and the output from the Program
+is covered only if its contents constitute a work based on the
+Program (independent of having been made by running the Program).
+Whether that is true depends on what the Program does.
+
+ 1. You may copy and distribute verbatim copies of the Program's
+source code as you receive it, in any medium, provided that you
+conspicuously and appropriately publish on each copy an appropriate
+copyright notice and disclaimer of warranty; keep intact all the
+notices that refer to this License and to the absence of any warranty;
+and give any other recipients of the Program a copy of this License
+along with the Program.
+
+You may charge a fee for the physical act of transferring a copy, and
+you may at your option offer warranty protection in exchange for a fee.
+
+ 2. You may modify your copy or copies of the Program or any portion
+of it, thus forming a work based on the Program, and copy and
+distribute such modifications or work under the terms of Section 1
+above, provided that you also meet all of these conditions:
+
+ a) You must cause the modified files to carry prominent notices
+ stating that you changed the files and the date of any change.
+
+ b) You must cause any work that you distribute or publish, that in
+ whole or in part contains or is derived from the Program or any
+ part thereof, to be licensed as a whole at no charge to all third
+ parties under the terms of this License.
+
+ c) If the modified program normally reads commands interactively
+ when run, you must cause it, when started running for such
+ interactive use in the most ordinary way, to print or display an
+ announcement including an appropriate copyright notice and a
+ notice that there is no warranty (or else, saying that you provide
+ a warranty) and that users may redistribute the program under
+ these conditions, and telling the user how to view a copy of this
+ License. (Exception: if the Program itself is interactive but
+ does not normally print such an announcement, your work based on
+ the Program is not required to print an announcement.)
+
+These requirements apply to the modified work as a whole. If
+identifiable sections of that work are not derived from the Program,
+and can be reasonably considered independent and separate works in
+themselves, then this License, and its terms, do not apply to those
+sections when you distribute them as separate works. But when you
+distribute the same sections as part of a whole which is a work based
+on the Program, the distribution of the whole must be on the terms of
+this License, whose permissions for other licensees extend to the
+entire whole, and thus to each and every part regardless of who wrote it.
+
+Thus, it is not the intent of this section to claim rights or contest
+your rights to work written entirely by you; rather, the intent is to
+exercise the right to control the distribution of derivative or
+collective works based on the Program.
+
+In addition, mere aggregation of another work not based on the Program
+with the Program (or with a work based on the Program) on a volume of
+a storage or distribution medium does not bring the other work under
+the scope of this License.
+
+ 3. You may copy and distribute the Program (or a work based on it,
+under Section 2) in object code or executable form under the terms of
+Sections 1 and 2 above provided that you also do one of the following:
+
+ a) Accompany it with the complete corresponding machine-readable
+ source code, which must be distributed under the terms of Sections
+ 1 and 2 above on a medium customarily used for software interchange; or,
+
+ b) Accompany it with a written offer, valid for at least three
+ years, to give any third party, for a charge no more than your
+ cost of physically performing source distribution, a complete
+ machine-readable copy of the corresponding source code, to be
+ distributed under the terms of Sections 1 and 2 above on a medium
+ customarily used for software interchange; or,
+
+ c) Accompany it with the information you received as to the offer
+ to distribute corresponding source code. (This alternative is
+ allowed only for noncommercial distribution and only if you
+ received the program in object code or executable form with such
+ an offer, in accord with Subsection b above.)
+
+The source code for a work means the preferred form of the work for
+making modifications to it. For an executable work, complete source
+code means all the source code for all modules it contains, plus any
+associated interface definition files, plus the scripts used to
+control compilation and installation of the executable. However, as a
+special exception, the source code distributed need not include
+anything that is normally distributed (in either source or binary
+form) with the major components (compiler, kernel, and so on) of the
+operating system on which the executable runs, unless that component
+itself accompanies the executable.
+
+If distribution of executable or object code is made by offering
+access to copy from a designated place, then offering equivalent
+access to copy the source code from the same place counts as
+distribution of the source code, even though third parties are not
+compelled to copy the source along with the object code.
+
+ 4. You may not copy, modify, sublicense, or distribute the Program
+except as expressly provided under this License. Any attempt
+otherwise to copy, modify, sublicense or distribute the Program is
+void, and will automatically terminate your rights under this License.
+However, parties who have received copies, or rights, from you under
+this License will not have their licenses terminated so long as such
+parties remain in full compliance.
+
+ 5. You are not required to accept this License, since you have not
+signed it. However, nothing else grants you permission to modify or
+distribute the Program or its derivative works. These actions are
+prohibited by law if you do not accept this License. Therefore, by
+modifying or distributing the Program (or any work based on the
+Program), you indicate your acceptance of this License to do so, and
+all its terms and conditions for copying, distributing or modifying
+the Program or works based on it.
+
+ 6. Each time you redistribute the Program (or any work based on the
+Program), the recipient automatically receives a license from the
+original licensor to copy, distribute or modify the Program subject to
+these terms and conditions. You may not impose any further
+restrictions on the recipients' exercise of the rights granted herein.
+You are not responsible for enforcing compliance by third parties to
+this License.
+
+ 7. If, as a consequence of a court judgment or allegation of patent
+infringement or for any other reason (not limited to patent issues),
+conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License. If you cannot
+distribute so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you
+may not distribute the Program at all. For example, if a patent
+license would not permit royalty-free redistribution of the Program by
+all those who receive copies directly or indirectly through you, then
+the only way you could satisfy both it and this License would be to
+refrain entirely from distribution of the Program.
+
+If any portion of this section is held invalid or unenforceable under
+any particular circumstance, the balance of the section is intended to
+apply and the section as a whole is intended to apply in other
+circumstances.
+
+It is not the purpose of this section to induce you to infringe any
+patents or other property right claims or to contest validity of any
+such claims; this section has the sole purpose of protecting the
+integrity of the free software distribution system, which is
+implemented by public license practices. Many people have made
+generous contributions to the wide range of software distributed
+through that system in reliance on consistent application of that
+system; it is up to the author/donor to decide if he or she is willing
+to distribute software through any other system and a licensee cannot
+impose that choice.
+
+This section is intended to make thoroughly clear what is believed to
+be a consequence of the rest of this License.
+
+ 8. If the distribution and/or use of the Program is restricted in
+certain countries either by patents or by copyrighted interfaces, the
+original copyright holder who places the Program under this License
+may add an explicit geographical distribution limitation excluding
+those countries, so that distribution is permitted only in or among
+countries not thus excluded. In such case, this License incorporates
+the limitation as if written in the body of this License.
+
+ 9. The Free Software Foundation may publish revised and/or new versions
+of the General Public License from time to time. Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+Each version is given a distinguishing version number. If the Program
+specifies a version number of this License which applies to it and "any
+later version", you have the option of following the terms and conditions
+either of that version or of any later version published by the Free
+Software Foundation. If the Program does not specify a version number of
+this License, you may choose any version ever published by the Free Software
+Foundation.
+
+ 10. If you wish to incorporate parts of the Program into other free
+programs whose distribution conditions are different, write to the author
+to ask for permission. For software which is copyrighted by the Free
+Software Foundation, write to the Free Software Foundation; we sometimes
+make exceptions for this. Our decision will be guided by the two goals
+of preserving the free status of all derivatives of our free software and
+of promoting the sharing and reuse of software generally.
+
+ NO WARRANTY
+
+ 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN
+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
+PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
+OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
+MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS
+TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE
+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
+REPAIR OR CORRECTION.
+
+ 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
+INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
+OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
+TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
+YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
+PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGES.
+
+ END OF TERMS AND CONDITIONS
+
+ How to Apply These Terms to Your New Programs
+
+ If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+ To do so, attach the following notices to the program. It is safest
+to attach them to the start of each source file to most effectively
+convey the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+ <one line to give the program's name and a brief idea of what it does.>
+ Copyright (C) <year> <name of author>
+
+ 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+
+Also add information on how to contact you by electronic and paper mail.
+
+If the program is interactive, make it output a short notice like this
+when it starts in an interactive mode:
+
+ Gnomovision version 69, Copyright (C) year name of author
+ Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+ This is free software, and you are welcome to redistribute it
+ under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License. Of course, the commands you use may
+be called something other than `show w' and `show c'; they could even be
+mouse-clicks or menu items--whatever suits your program.
+
+You should also get your employer (if you work as a programmer) or your
+school, if any, to sign a "copyright disclaimer" for the program, if
+necessary. Here is a sample; alter the names:
+
+ Yoyodyne, Inc., hereby disclaims all copyright interest in the program
+ `Gnomovision' (which makes passes at compilers) written by James Hacker.
+
+ <signature of Ty Coon>, 1 April 1989
+ Ty Coon, President of Vice
+
+This General Public License does not permit incorporating your program into
+proprietary programs. If your program is a subroutine library, you may
+consider it more useful to permit linking proprietary applications with the
+library. If this is what you want to do, use the GNU Library General
+Public License instead of this License.
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..b6711ef
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,11 @@
+Package: coda
+Version: 0.9-1
+Date: 2004/12/20
+Title: Output analysis and diagnostics for MCMC
+Author: Martyn Plummer, Nicky Best, Kate Cowles, Karen Vines
+Maintainer: Martyn Plummer <plummer at iarc.fr>
+Depends: R (>= 1.9)
+Description: Output analysis and diagnostics for Markov Chain Monte Carlo simulations.
+License: GPL Version 2 or later.
+URL: http://www-fis.iarc.fr/coda/
+Packaged: Mon Dec 20 12:22:23 2004; martyn
diff --git a/INDEX b/INDEX
new file mode 100644
index 0000000..abf2b96
--- /dev/null
+++ b/INDEX
@@ -0,0 +1,44 @@
+[.mcmc Extract or replace parts of MCMC objects
+as.matrix.mcmc Conversions of MCMC objects
+as.ts.mcmc Coerce mcmc object to time series
+autocorr Autocorrelation function for Markov chains
+autocorr.plot Plot autocorrelations for Markov Chains
+bugs2jags Convert WinBUGS data file to JAGS data file
+coda.options Options settings for the codamenu driver
+codamenu Main menu driver for the coda package
+crosscorr Cross correlations for MCMC output
+crosscorr.plot Plot image of correlation matrix
+cumuplot Cumulative quantile plot
+densplot Probability density function estimate from
+ MCMC output
+effectiveSize Effective sample size for estimating the mean
+gelman.diag Gelman and Rubin's convergence diagnostic
+gelman.plot Gelman-Rubin-Brooks plot
+geweke.diag Geweke's convergence diagnostic
+geweke.plot Geweke-Brooks plot
+heidel.diag Heidelberger and Welch's convergence
+ diagnostic
+line Simple linear regression example
+mcmc Markov Chain Monte Carlo Objects
+mcmc.list Replicated Markov Chain Monte Carlo Objects
+mcmcUpgrade Upgrade mcmc objects in obsolete format
+mcpar Mcpar attribute of MCMC objects
+multi.menu Choose multiple options from a menu
+nchain Dimensions of MCMC objects
+pcramer The Cramer-von Mises Distribution
+plot.mcmc Summary plots of mcmc objects
+raftery.diag Raftery and Lewis's diagnostic
+read.and.check Read data interactively and check that it
+ satisfies conditions
+read.coda Read output files in CODA format
+read.coda.interactive Read CODA output files interactively
+read.openbugs Read CODA output files produced by OpenBUGS
+spectrum0 Estimate spectral density at zero
+spectrum0.ar Estimate spectral density at zero
+summary.mcmc Summary statistics for Markov Chain Monte
+ Carlo chains
+thin Thinning interval
+time.mcmc Time attributes for mcmc objects
+traceplot Trace plot of MCMC output
+varnames Named dimensions of MCMC objects
+window.mcmc Time windows for mcmc objects
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..2216469
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,95 @@
+## export functions
+export(
+ as.array.mcmc.list,
+ as.mcmc,
+ as.ts.mcmc,
+ autocorr,
+ autocorr.plot,
+ bugs2jags,
+ chanames,
+ coda.options,
+ codamenu,
+ crosscorr,
+ crosscorr.plot,
+ cumuplot,
+ densplot,
+ display.coda.options,
+ effectiveSize,
+ gelman.diag,
+ gelman.plot,
+ geweke.diag,
+ geweke.plot,
+ heidel.diag,
+ is.mcmc,
+ is.mcmc.list,
+ mcmc,
+ mcmc.list,
+ mcmcUpgrade,
+ mcpar,
+ multi.menu,
+ nchain,
+ niter,
+ nvar,
+ pcramer,
+ raftery.diag,
+ read.and.check,
+ read.coda.interactive,
+ read.coda,
+ read.jags,
+ read.openbugs,
+ spectrum0,
+ spectrum0.ar,
+ thin,
+ traceplot,
+ varnames,
+ "varnames<-")
+
+## export data
+# export(line)
+# this doesn't need to be exported, not sure why, but it works fine
+# as is.
+
+# [1] [.mcmc as.matrix.mcmc end.mcmc frequency.mcmc plot.mcmc
+# [6] print.mcmc start.mcmc summary.mcmc thin.mcmc time.mcmc
+# [11] window.mcmc
+
+# [1] [.mcmc.list as.matrix.mcmc.list as.mcmc.mcmc.list
+# [4] end.mcmc.list plot.mcmc.list start.mcmc.list
+# [7] summary.mcmc.list thin.mcmc.list time.mcmc.list
+# [10] window.mcmc.list
+
+## register S3 methods
+## mcmc
+S3method("[", mcmc)
+S3method(as.matrix, mcmc)
+S3method(end, mcmc)
+S3method(frequency, mcmc)
+S3method(plot, mcmc)
+S3method(print, mcmc)
+S3method(start, mcmc)
+S3method(summary, mcmc)
+S3method(thin, mcmc)
+S3method(time, mcmc)
+S3method(window, mcmc)
+
+## mcmc.list
+S3method("[", mcmc.list)
+S3method(as.matrix, mcmc.list)
+S3method(as.mcmc, mcmc.list)
+S3method(end, mcmc.list)
+S3method(plot, mcmc.list)
+S3method(start, mcmc.list)
+S3method(summary, mcmc.list)
+S3method(thin, mcmc.list)
+S3method(time, mcmc.list)
+S3method(window, mcmc.list)
+
+## misc
+S3method(print, gelman.diag)
+S3method(print, geweke.diag)
+S3method(print, heidel.diag)
+S3method(print, raftery.diag)
+S3method(print, summary.mcmc)
+
+## Imports
+importFrom(stats, start, end, frequency, time, window)
diff --git a/R/codamenu.R b/R/codamenu.R
new file mode 100644
index 0000000..011a0d4
--- /dev/null
+++ b/R/codamenu.R
@@ -0,0 +1,963 @@
+"codamenu" <- function ()
+{
+ on.exit(tidy.up())
+ coda.options(default=TRUE)
+ file.menu <- c("Read BUGS output files",
+ "Use an mcmc object",
+ "Quit")
+ pick <- menu(file.menu, title = "CODA startup menu")
+ if (pick == 0 || pick == 3)
+ return(invisible())
+ else if (pick == 1) {
+ assign("coda.dat", read.coda.interactive(), pos=1)
+ if (is.null(coda.dat)) {
+ return(invisible())
+ }
+ coda.options(data.saved = FALSE)
+ }
+ else if (pick == 2) {
+ msg <- "\nEnter name of saved object (or type \"exit\" to quit)"
+ repeat {
+ cat(msg, "\n")
+ outname <- read.and.check(what = character())
+ if (outname == "exit" || outname == "\"exit\"") {
+ return(invisible())
+ }
+ else if (!exists(outname))
+ msg <- "Can't find this object"
+ else {
+ work.dat <- eval(parse(text = outname))
+ if (is.mcmc.list(work.dat)) {
+ assign("coda.dat", work.dat, pos=1)
+ break
+ }
+ else if (is.mcmc(work.dat)) {
+ assign("coda.dat", mcmc.list(work.dat), pos=1)
+ break
+ }
+ else
+ msg <- "Not an mcmc or mcmc.list object"
+ }
+ }
+ }
+ else stop("Invalid option")
+ coda.dat <- coda.dat #Create local copy
+ if (is.null(chanames(coda.dat))) {
+ chanames(coda.dat) <- chanames(coda.dat, allow.null = FALSE)
+ }
+ if (is.matrix(coda.dat[[1]]) && is.null(varnames(coda.dat))) {
+ varnames(coda.dat) <- varnames(coda.dat, allow.null = FALSE)
+ }
+ assign("coda.dat", coda.dat, pos=1)
+ rm(coda.dat, inherits=FALSE) #Destroy local copy
+ assign("work.dat", coda.dat, pos=1)
+
+ ## Check for variables that are linear functions of the
+ ## iteration number
+ is.linear <- rep(FALSE, nvar(coda.dat))
+ for (i in 1:nchain(coda.dat)) {
+ for (j in 1:nvar(coda.dat)) {
+ lm.out <- lm(as.matrix(coda.dat[[i]])[,j] ~ time(coda.dat))
+ if (identical(all.equal(var(residuals(lm.out)), 0), TRUE)) {
+ is.linear[j] <- TRUE
+ }
+ }
+ }
+ if (any(is.linear)) {
+ cat("Dropping the following variables, which are linear\n")
+ cat("functions of the iteration number\n")
+ print(varnames(coda.dat)[is.linear])
+ inset <- varnames(coda.dat)[!is.linear]
+ assign("coda.dat", coda.dat[,inset, drop=FALSE], pos=1)
+ assign("work.dat", coda.dat[,inset, drop=FALSE], pos=1)
+ }
+
+ ## Sample size test
+ cat("Checking effective sample size ...")
+ ess <- effectiveSize(gelman.transform(coda.dat))
+ warn.small <- FALSE
+ for (i in 1:nchain(coda.dat))
+ {
+ if (any(ess[[i]] < 200))
+ warn.small <- TRUE
+ }
+ if (warn.small)
+ {
+ cat("\n")
+ cat("*******************************************\n")
+ cat("WARNING !!! \n")
+ cat("Some variables in your chain have an \n")
+ cat("effective sample size of less than 200 \n")
+ cat("This is too small, and may cause errors \n")
+ cat("in the diagnostic tests \n")
+ cat("HINT: \n")
+ cat("Look at plots first to identify variables\n")
+ cat("with slow mixing. (Choose menu Output \n")
+ cat("Analysis then Plots) \n")
+ cat("Re-run your chain with a larger sample \n")
+ cat("size and thinning interval. If possible, \n")
+ cat("reparameterize your model to improve mixing\n")
+ cat("*******************************************\n")
+ }
+ else
+ {
+ cat("OK\n")
+ }
+ current.menu <- "codamenu.main"
+ old.opt <- options(warn=-1, show.error.messages=FALSE)
+ on.exit(options(old.opt))
+ repeat {
+ next.menu <- try(do.call(current.menu, vector("list", 0)))
+ if (!is.null(class(next.menu)) && class(next.menu) == "try-error")
+ {
+ if (current.menu == "codamenu.main")
+ {
+ cat("A crash has occurred in the main menu\nBailing out\n")
+ return();
+ }
+ else
+ {
+ cat("\n\n")
+ cat("**********************\n")
+ cat("An error has occurred\n")
+ cat("Returning to main menu\n")
+ cat("**********************\n")
+ current.menu <- "codamenu.main"
+ }
+ }
+ else
+ {
+ if (next.menu == "quit") {
+ if(read.yesno("Are you sure you want to quit", FALSE))
+ break
+ }
+ else current.menu <- next.menu
+ }
+ }
+ invisible()
+}
+
+"codamenu.anal" <- function ()
+{
+ next.menu <- "codamenu.anal"
+ choices <- c("Plots", "Statistics", "List/Change Options",
+ "Return to Main Menu")
+ next.menu.list <- c("plots", "summary", "codamenu.options",
+ "codamenu.main")
+ cat("\n")
+ pick <- menu(choices, title = "CODA Output Analysis menu")
+ if (pick == 0)
+ next.menu <- "quit"
+ else if (next.menu.list[pick] == "summary") {
+ if (coda.options("combine.stats")) {
+ print(summary(work.dat, quantiles = coda.options("quantiles"),
+ digits = coda.options("digits")))
+ }
+ else for (i in 1:nchain(work.dat)) {
+ cat(chanames(work.dat, allow.null = FALSE)[i], "\n")
+ print(summary(work.dat[[i]], quantiles = coda.options("quantiles"),
+ digits = coda.options("digits")))
+ }
+ }
+ else if (next.menu.list[pick] == "plots") {
+ auto.layout <- !coda.options("user.layout")
+ ask <- TRUE
+ repeat {
+ if (coda.options("combine.plots"))
+ plot(work.dat, trace = coda.options("trace"),
+ density = coda.options("densplot"),
+ smooth = coda.options("lowess"),
+ auto.layout = auto.layout, bwf = coda.options("bandwidth"),
+ combine.chains = !coda.options("combine.plots"),
+ ask = ask)
+ else for (i in 1:nchain(work.dat)) {
+ plot(work.dat[[i]], trace = coda.options("trace"),
+ density = coda.options("densplot"),
+ smooth = coda.options("lowess"),
+ auto.layout = auto.layout, bwf = coda.options("bandwidth"),
+ combine.chains = coda.options("combine.plots"),
+ ask = ask)
+ }
+ codamenu.ps()
+ if (names(dev.cur()) == "postscript")
+ ask <- FALSE
+ else break
+ }
+ }
+ else next.menu <- next.menu.list[pick]
+ return(next.menu)
+}
+
+"codamenu.diags" <- function ()
+{
+ next.menu <- "diags"
+ while (next.menu == "diags") {
+ choices <- c("Geweke", "Gelman and Rubin", "Raftery and Lewis",
+ "Heidelberger and Welch", "Autocorrelations",
+ "Cross-Correlations", "List/Change Options",
+ "Return to Main Menu")
+ next.menu.list <- c("codamenu.diags.geweke", "codamenu.diags.gelman",
+ "codamenu.diags.raftery", "codamenu.diags.heidel",
+ "codamenu.diags.autocorr", "codamenu.diags.crosscorr",
+ "codamenu.options", "codamenu.main")
+ pick <- menu(choices, title = "CODA Diagnostics Menu")
+ if (pick == 0)
+ return("quit")
+ else next.menu <- next.menu.list[pick]
+ }
+ return(next.menu)
+}
+
+"codamenu.diags.autocorr" <- function ()
+{
+ next.menu <- "codamenu.diags"
+ codamenu.output.header("AUTOCORRELATIONS WITHIN EACH CHAIN:")
+ print(autocorr(work.dat), digits = coda.options("digits"))
+ choices <- c("Plot autocorrelations", "Return to Diagnostics Menu")
+ pick <- menu(choices, title = "Autocorrelation Plots Menu")
+ if (pick == 0)
+ next.menu <- "quit"
+ else if (pick == 1) {
+ ask <- TRUE
+ repeat {
+ autocorr.plot(work.dat, auto.layout = !coda.options("user.layout"),
+ ask = ask)
+ codamenu.ps()
+ if (names(dev.cur()) == "postscript")
+ ask <- FALSE
+ else break
+ }
+ }
+ return(next.menu)
+}
+
+"codamenu.diags.crosscorr" <- function ()
+{
+ next.menu <- "codamenu.diags.crosscorr"
+ crosscorr.out <- if (coda.options("combine.corr")) {
+ crosscorr(work.dat)
+ }
+ else lapply(work.dat, crosscorr)
+ if (coda.options("combine.corr") & nchain(work.dat) > 1)
+ cat("Pooling over chains:", chanames(work.dat, allow.null = FALSE),
+ sep = "\n", collapse = "\n")
+ print(crosscorr.out, digits = coda.options("digits"))
+ cat("\n")
+ choices <- c("Change options",
+ "Plot Cross Correlations",
+ "Return to Diagnostics Menu")
+ pick <- menu(choices, title = "Cross correlation plots menu")
+ if (pick == 0)
+ next.menu <- "quit"
+ else
+ switch(pick,
+ change.tfoption("Combine chains", "combine.corr"),
+ {
+ repeat {
+ if (coda.options("combine.corr"))
+ crosscorr.plot(work.dat)
+ else {
+ opar <- par(ask = TRUE)
+ lapply(work.dat, crosscorr.plot)
+ par(opar)
+ }
+ codamenu.ps()
+ if (names(dev.cur()) != "postscript")
+ break
+ }
+ },
+ next.menu <- "codamenu.diags")
+ return(next.menu)
+}
+
+"codamenu.diags.heidel" <- function ()
+{
+ this.menu <- "codamenu.diags.heidel"
+ next.menu <- "codamenu.diags"
+ title <- "HEIDELBERGER AND WELCH STATIONARITY AND INTERVAL HALFWIDTH TESTS"
+ codamenu.output.header(title)
+ cat("Precision of halfwidth test =", coda.options("halfwidth"), "\n\n")
+ heidel.out <- heidel.diag(work.dat, eps = coda.options("halfwidth"))
+ print(heidel.out, digits = coda.options("digits"))
+ choices <- c("Change precision", "Return to diagnostics menu")
+ pick <- menu(choices)
+ if (pick == 0)
+ next.menu <- "quit"
+ else if (pick == 1)
+ next.menu <- codamenu.options.heidel(this.menu)
+ return(next.menu)
+}
+
+"codamenu.diags.raftery" <- function ()
+{
+ next.menu <- this.menu <- "codamenu.diags.raftery"
+ codamenu.output.header("RAFTERY AND LEWIS CONVERGENCE DIAGNOSTIC")
+ print(raftery.diag(work.dat, q = coda.options("q"), r = coda.options("r"),
+ s = coda.options("s")), digits = coda.options("digits"))
+ choices <- c("Change parameters", "Return to diagnostics menu")
+ pick <- menu(choices)
+ next.menu <- if (pick == 0)
+ "quit"
+ else if (pick == 1) {
+ codamenu.options.raftery(this.menu)
+ }
+ else "codamenu.diags"
+ return(next.menu)
+}
+
+"codamenu.main" <- function ()
+{
+ choices <- c("Output Analysis", "Diagnostics", "List/Change Options", "Quit")
+ next.menu.list <- c("codamenu.anal", "codamenu.diags", "codamenu.options",
+ "quit")
+ pick <- menu(choices, title = "CODA Main Menu")
+ if (pick == 0)
+ next.menu <- "quit"
+ else next.menu <- next.menu.list[pick]
+ return(next.menu)
+}
+
+"codamenu.diags.gelman" <- function (tol = 1e-08)
+{
+ next.menu <- this.menu <- "codamenu.diags.gelman"
+ if (nchain(work.dat) == 1) {
+ cat("\nError: you need more than one chain.\n\n")
+ return(next.menu = "codamenu.diags")
+ }
+ else if (niter(work.dat) <= 50) {
+ cat("\nError: you need > 50 iterations in the working data\n")
+ return(next.menu = "codamenu.diags")
+ }
+ z <- window(work.dat, start = niter(work.dat)/2)
+ for (i in 2:nchain(z)) {
+ for (j in 1:(i - 1)) {
+ if (any(apply(as.matrix(z[[i]] - z[[j]]), 2, var)) < tol) {
+ cat("\nError: 2nd halves of",
+ chanames(z, allow.null = FALSE)[c(j, i)],
+ "are identical for at least one variable\n")
+ return(next.menu = "codamenu.diags")
+ }
+ }
+ }
+ codamenu.output.header("GELMAN AND RUBIN DIAGNOSTIC")
+ print(gelman.diag(work.dat, transform = TRUE),
+ digits = coda.options("digits"))
+ choices <- c("Shrink Factor Plots", "Change bin size for shrink plot",
+ "Return to Diagnostics Menu")
+ action.list <- c("ShrinkPlot", "ChangeBin", "Return")
+ while (next.menu == "codamenu.diags.gelman") {
+ pick <- menu(choices, title = "Gelman & Rubin menu")
+ if (pick == 0)
+ next.menu <- "quit"
+ else switch(action.list[pick], ShrinkPlot = {
+ ask <- TRUE
+ repeat {
+ gelman.plot(work.dat, max.bins = coda.options("gr.max"),
+ bin.width = coda.options("gr.bin"),
+ auto.layout = !coda.options("user.layout"),
+ ask = ask)
+ codamenu.ps()
+ if (names(dev.cur()) == "postscript")
+ ask <- FALSE
+ else break
+ }
+ }, ChangeBin = {
+ codamenu.options.gelman(NULL)
+ }, Return = {
+ next.menu <- "codamenu.diags"
+ })
+ }
+ return(next.menu)
+}
+
+"codamenu.diags.geweke" <- function ()
+{
+ next.menu <- "codamenu.diags.geweke"
+ codamenu.output.header("GEWEKE CONVERGENCE DIAGNOSTIC (Z-score)")
+ geweke.out <- geweke.diag(work.dat, frac1 = coda.options("frac1"),
+ frac2 = coda.options("frac2"))
+ print(geweke.out, digits = coda.options("digits"))
+ choices <- c("Change window size", "Plot Z-scores",
+ "Change number of bins for plot",
+ "Return to Diagnostics Menu")
+ action.list <- c("ChangeWindow", "Plot", "ChangeBin", "Return")
+ while (next.menu == "codamenu.diags.geweke") {
+ pick <- menu(choices, title = "Geweke plots menu")
+ if (pick == 0)
+ return("quit")
+ switch(action.list[pick], ChangeWindow = {
+ codamenu.options.geweke.win(NULL)
+ geweke.out <- geweke.diag(work.dat,
+ frac1 = coda.options("frac1"),
+ frac2 = coda.options("frac2"))
+ print(geweke.out, digits = coda.options("digits"))
+ }, Plot = {
+ ask <- TRUE
+ repeat {
+ if(start(work.dat) >= end(work.dat)) {
+ cat("Chain too short: end iteration must be at least twice\n")
+ cat("the start iteration\n")
+ break
+ }
+ geweke.plot(work.dat, frac1 = coda.options("frac1"),
+ frac2 = coda.options("frac2"),
+ nbins = coda.options("geweke.nbin"),
+ auto.layout = !coda.options("user.layout"),
+ ask = ask)
+ codamenu.ps()
+ if (names(dev.cur()) == "postscript")
+ ask <- FALSE
+ else break
+ }
+ }, ChangeBin = {
+ codamenu.options.geweke.bin(NULL)
+ }, Return = {
+ next.menu <- "codamenu.diags"
+ })
+ }
+ return(next.menu)
+}
+
+"codamenu.options" <- function ()
+{
+ next.menu <- "codamenu.options"
+ choices <- c("List current options", "Data Options", "Plot Options",
+ "Summary Statistics Options", "Diagnostics Options",
+ "Output Analysis", "Diagnostics", "Main Menu")
+ action.list <- c("ListOptions", "codamenu.options.data",
+ "codamenu.options.plot", "codamenu.options.stats",
+ "codamenu.options.diag", "codamenu.anal",
+ "codamenu.diags", "codamenu.main")
+ pick <- menu(choices, title = "CODA main options menu")
+ if (pick == 0)
+ return("quit")
+ if (action.list[pick] == "ListOptions") {
+ display.coda.options(data = TRUE, stats = TRUE, plots = TRUE,
+ diags = TRUE)
+ next.menu <- "codamenu.options"
+ }
+ else next.menu <- action.list[pick]
+ return(next.menu)
+}
+
+"codamenu.options.data" <- function ()
+{
+ next.menu <- "codamenu.options.data"
+
+ work.vars <- varnames(work.dat)
+ work.chains <- chanames(work.dat)
+ work.start <- start(work.dat)
+ work.end <- end(work.dat)
+ work.thin <- thin(work.dat)
+
+ choices <- c("List current data options", "Select variables for analysis",
+ "Select chains for analysis", "Select iterations for analysis",
+ "Select thinning interval", "Return to main options menu")
+ action.list <- c("ListDataOptions", "SelectVars", "SelectChains",
+ "SelectIters", "SelectThinInterval", "MainOptionsMenu")
+ pick <- menu(choices, title = "CODA data options menu")
+ if (pick == 0)
+ return("quit")
+ switch(action.list[pick], ListDataOptions = {
+ display.coda.options(data = TRUE)
+ }, SelectVars = {
+ work.vars <- multi.menu(varnames(coda.dat, allow.null = FALSE),
+ "Select variables for analysis",
+ c("VARIABLE NUMBER", "VARIABLE NAME"),
+ allow.zero = FALSE)
+ }, SelectChains = {
+ work.chains <- multi.menu(chanames(coda.dat, allow.null = FALSE),
+ "Select chains for analysis:",
+ c("CHAIN NUMBER", "CHAIN NAME"),
+ allow.zero = FALSE)
+ }, SelectIters = {
+ cat("\nIterations available = ", start(coda.dat), ":",
+ end(coda.dat), "\n", sep = "")
+ work.start <- read.and.check("Enter iteration you wish to start at",
+ lower = start(coda.dat),
+ upper = end(coda.dat),
+ default = start(work.dat))
+ work.end <- read.and.check("Enter iteration you wish to end at",
+ lower = work.start,
+ upper = end(coda.dat),
+ default = end(work.dat))
+ }, SelectThinInterval = {
+ cat("\nThinning interval of full data = ", thin(coda.dat), "\n", sep = "")
+ work.thin <- read.and.check("Enter thinning interval:",
+ lower = thin(coda.dat),
+ default = thin(work.dat))
+ }, MainOptionsMenu = {
+ next.menu <- "codamenu.options"
+ })
+ if (action.list[pick] != "ListDataOptions" && action.list[pick] !=
+ "MainOptionsMenu") {
+ cat("Recreating working data...\n")
+ wd <- window(coda.dat[, work.vars, drop = FALSE], start = work.start,
+ end = work.end, thin = work.thin)
+ assign("work.dat", wd[work.chains, drop=FALSE], pos=1)
+ }
+ return(next.menu)
+}
+
+"codamenu.options.diag" <- function ()
+{
+ next.menu <- this.menu <- "codamenu.options.diag"
+ choices <- c("Display current diagnostic options",
+ "Window sizes for Geweke's diagnostic",
+ "Bin size for plotting Geweke's diagnostic",
+ "Bin size for plotting Gelman & Rubin's diagnostic",
+ "Parameters for Raftery & Lewis' diagnostic",
+ "Halfwidth precision for Heidelberger & Welch's diagnostic",
+ "Combine chains to calculate correlation matrix",
+ "Return to main options menu")
+ pick <- menu(choices, title = "CODA diagnostics options menu")
+ if (pick == 0)
+ return("quit")
+ switch(pick, display.coda.options(diags = TRUE),
+ next.menu <- codamenu.options.geweke.win(this.menu),
+ next.menu <- codamenu.options.geweke.bin(this.menu),
+ next.menu <- codamenu.options.gelman(this.menu),
+ next.menu <- codamenu.options.raftery(this.menu),
+ next.menu <- codamenu.options.heidel(this.menu),
+ {
+ change.tfoption("Do you want to combine all chains to calculate correlation matrix", "combine.corr")
+ }, next.menu <- "codamenu.options")
+ return(next.menu)
+}
+
+"codamenu.options.gelman" <- function (last.menu)
+{
+ choices <- c("Default: bin width = 10; maximum number of bins = 50",
+ "User-specified bin width",
+ "User-specified total number of bins")
+ pick <- menu(choices, title = "Options for defining bin size to plot Gelman-Rubin-Brooks diagnostic")
+ if (pick == 0)
+ return("quit")
+ switch(pick, {
+ coda.options(gr.max = 50)
+ coda.options(gr.bin = 10)
+ }, {
+ coda.options(gr.max = Inf)
+ default <- if (coda.options("gr.bin") == 0)
+ 10
+ else coda.options("gr.bin")
+ msg <- "Enter required bin width:"
+ coda.options(gr.bin = read.and.check(msg, lower = 1,
+ upper = niter(work.dat) - 50, default = default))
+ }, {
+ coda.options(gr.bin = 0)
+ default <- if (is.infinite(coda.options("gr.max")))
+ 50
+ else coda.options("gr.max")
+ msg <- "Enter total number of bins required:"
+ coda.options(gr.max = read.and.check(msg, lower = 1,
+ upper = niter(work.dat) - 50, default = default))
+ })
+ return(last.menu)
+}
+
+"codamenu.options.geweke.bin" <- function (last.menu)
+{
+ msg <- "Enter number of bins for Geweke-Brooks plot"
+ ans <- read.and.check(msg, what=numeric(), lower=1,
+ default=coda.options("geweke.nbin"))
+ coda.options(geweke.nbin = ans)
+ return(last.menu)
+}
+
+"codamenu.options.geweke.win" <-
+function (last.menu)
+{
+ msg1 <- "Enter fraction of chain to include in 1st window:"
+ msg2 <- "Enter fraction of chain to include in 2nd window:"
+ ans1 <- ans2 <- 1
+ while (ans1 + ans2 >= 1) {
+ ans1 <- read.and.check(msg1, lower = 0, upper = 1,
+ default = coda.options("frac1"))
+ ans2 <- read.and.check(msg2, lower = 0, upper = 1,
+ default = coda.options("frac2"))
+ ## Check that sum of fractions doesn't exceed 1.0
+ if (ans1 + ans2 >= 1)
+ cat("Error: Sum of fractions in 1st and 2nd windows must be < 1.0\n")
+ }
+ coda.options(frac1 = ans1, frac2 = ans2)
+ return(last.menu)
+}
+
+"codamenu.options.heidel" <- function (last.menu)
+{
+ coda.options(halfwidth = read.and.check("Enter precision for halfwidth test",
+ lower = 0, default = coda.options("halfwidth")))
+ return(last.menu)
+}
+
+"codamenu.options.plot" <-
+function ()
+{
+ next.menu <- "codamenu.options.plot"
+ choices <- c("Show current plotting options",
+ "Plot trace of samples",
+ "Plot kernel density estimate",
+ "Add smooth line through trace plot",
+ "Combine chains",
+ "Single plot per page",
+ "Specify page layout for plots",
+ "Select bandwidth function for kernel smoothing",
+ "Return to main options menu")
+ pick <- menu(choices, title = "CODA plotting options menu")
+ if (pick == 0)
+ return("quit")
+ switch(pick,
+ display.coda.options(plots = TRUE),
+ change.tfoption(choices[2], "trace"),
+ change.tfoption(choices[3], "densplot"),
+ change.tfoption(choices[4], "lowess"),
+ change.tfoption(choices[5], "combine.plots"),
+ {
+ ans <- read.yesno(choices[6], default=TRUE)
+ if(ans) {
+ coda.options(user.layout = TRUE)
+ par(mfrow = c(1,1))
+ }
+ },
+ {
+ change.tfoption("Do you want to specify your own page layout for the plots", "user.layout")
+ if (coda.options("user.layout")) {
+ mrows <- read.and.check("Enter number of rows per page",
+ lower = 1, upper = 7)
+ mcols <- read.and.check("Enter number of columns per page",
+ lower = 1, upper = 8)
+ par(mfrow = c(mrows, mcols))
+ }
+ }, {
+ next.menu <- "codamenu.options.plot.kernel"
+ }, NULL)
+ if (pick == length(choices))
+ next.menu <- "codamenu.options"
+ return(next.menu)
+}
+
+"codamenu.options.plot.kernel" <-
+function ()
+{
+ if (!coda.options("densplot")) {
+ cat("\nNo density plots requested - this option is irrelevant\n")
+ }
+ else {
+ kernel.menu <- c("Smooth (0.25 * sample range)",
+ "Coarse (Silverman 1986 eqn. 3.28 & 3.30)",
+ "User-defined function",
+ "Return to Plotting Options Menu")
+ pick1 <- menu(kernel.menu, title = "Select kernel bandwidth function")
+ if (pick1 == 0)
+ return("quit")
+ switch(pick1, {
+ bwf <- function(x) {
+ (max(x) - min(x))/4
+ }
+ coda.options(bandwidth = bwf)
+ }, {
+ bwf <- function(x) {
+ 1.06 * min(sd(x), IQR(x)/1.34) * length(x)^-0.2
+ }
+ coda.options(bandwidth = bwf)
+ }, {
+ restart()
+ func.OK <- FALSE
+ while (!func.OK) {
+ cat("Enter bandwidth as an expression in terms of x,\n")
+ cat("the vector of sampled values, e.g. \n")
+ cat("(max(x) - min(x)) / 4\n")
+ ans <- scan(what = character())
+ if (length(ans) > 0) {
+ bwf <- "function(x){"
+ for (i in 1:length(ans)) {
+ bwf <- paste(bwf, ans[i], sep = "")
+ }
+ bwf <- paste(bwf, "}", sep = "")
+ bwf <- eval(parse(text = bwf))
+ ## Carry out simple test to check whether the
+ ## function entered makes sense
+ ##
+ func.OK <- is.numeric(bw <- bwf(1:10)) && (length(bw) == 1)
+ if(!func.OK) {
+ cat("This is not a suitable function: it must return a single\n")
+ cat("numeric value given a numeric vector x.")
+ }
+ }
+ }
+ coda.options(bandwidth = bwf)
+ }, NULL)
+ }
+ return("codamenu.options.plot")
+}
+
+"codamenu.options.plot.ps" <- function ()
+{
+ choices <- c("Portrait", "Landscape")
+ pick <- menu(choices, "Select options for saving plots to PostScript files")
+ if (pick == 0)
+ return("quit")
+ else coda.options(ps.orientation = c("portrait", "landscape")[pick])
+ if (.Device == "X11")
+ x11(orientation = coda.options("ps.orientation"))
+ else if (.Device == "Win32")
+ windows(orientation = coda.options("ps.orientation"))
+ return("codamenu.options.plot")
+}
+
+"codamenu.options.raftery" <- function (last.menu)
+{
+ coda.options(q = read.and.check("Enter quantile to be estimated:",
+ lower = 0, upper = 1, default = coda.options("q")))
+ coda.options(r = read.and.check("Enter required precision:",
+ upper = coda.options("q"), default = coda.options("r")))
+ coda.options(s = read.and.check("Enter required probability:",
+ lower = 0, upper = 1, default = coda.options("s")))
+ return(last.menu)
+}
+
+"codamenu.options.stats" <- function ()
+{
+ next.menu <- "codamenu.options.stats"
+ choices <- c("Display current statistics options", "Combine chains for summary statistics",
+ "Quantiles for summary statistics", "Number of significant digits for printing",
+ "Return to main options menu")
+ pick <- menu(choices, title = "CODA options for summary statistics")
+ if (pick == 0)
+ return("quit")
+ switch(pick, display.coda.options(stats = TRUE), {
+ mssg <- "Do you want to combine all chains when calculating summary statistics"
+ change.tfoption(mssg, "combine.stats")
+ }, {
+ mssg <- paste("Enter quantiles required, separated by commas\n(Default =",
+ paste(coda.options("quantiles"), collapse = ", "))
+ repeat {
+ cat("\n", mssg, "\n")
+ ans <- as.numeric(scan(what = character(), sep = ",",
+ quiet = TRUE, nlines = 1))
+ if (length(ans) == 0)
+ ans <- coda.options("quantiles")
+ if (any(is.na(ans)))
+ mssg <- "You must enter numeric values"
+ else if (any(ans >= 1) || any(ans <= 0))
+ mssg <- "You must enter values between 0 and 1"
+ else break
+ }
+ if (length(ans) > 0)
+ coda.options(quantiles = sort(ans))
+ }, {
+ mssg <- "Enter number of significant digits to be printed"
+ ans <- read.and.check(mssg, what = integer(), lower = 0,
+ default = coda.options("digits"))
+ coda.options(digits = ans)
+ }, {
+ next.menu <- "codamenu.options"
+ })
+ return(next.menu)
+}
+
+"display.coda.options" <-
+ function (data = FALSE, stats = FALSE, plots = FALSE, diags = FALSE)
+{
+ cat("\nCurrent option settings:")
+ cat("\n=======================\n\n")
+ if (data) {
+ cat("WORKING DATA\n")
+ cat("============\n")
+ cat("Variables selected : ",
+ paste(varnames(work.dat, allow.null = FALSE), collapse=", ")
+ ,"\n", sep="")
+ cat("Chains selected : ",
+ paste(chanames(work.dat, allow.null = FALSE), collapse=", ")
+ , "\n", sep="")
+ cat("Iterations - start : ", start(work.dat), "\n", sep="")
+ cat(" end : ", end(work.dat), "\n", sep="")
+ cat("Thinning interval : ", thin(work.dat), "\n", sep="")
+ cat("\n")
+ }
+ if (stats) {
+ cat("SUMMARY STATISTICS OPTIONS\n")
+ cat("==========================\n\n")
+ cat("Combine chains : ", coda.options("combine.stats"), "\n", sep="")
+ cat("Quantiles : ",
+ paste(coda.options("quantiles") * 100, "%", sep="", collapse = ", "),
+ "\n", sep="")
+ cat("Significant digits : ", coda.options("digits"), "\n", sep="")
+ cat("\n")
+ }
+ if (plots) {
+ cat("PLOTTING OPTIONS\n")
+ cat("================\n\n")
+ cat("Trace : ", coda.options("trace"), "\n", sep="")
+ cat("Density : ", coda.options("densplot"), "\n", sep="")
+ cat("Smooth lines : ", coda.options("lowess"), "\n", sep="")
+ cat("Combine chains : ", coda.options("combine.plots"), "\n", sep="")
+ cat("User-defined layout : ", coda.options("user.layout"), "\n", sep="")
+ if(coda.options("user.layout")) {
+ cat(" : ", paste(par("mfrow"), collapse=" X "), "\n", sep="")
+ }
+ cat("Bandwidth function :\n")
+ print(coda.options("bandwidth"))
+ cat("\n")
+ }
+ if (diags) {
+ cat("DIAGNOSTICS OPTIONS\n")
+ cat("===================\n\n")
+ cat("Geweke\n")
+ cat("------\n")
+ cat("Window 1 fraction : ", coda.options("frac1"), "\n", sep="")
+ cat("Window 2 fraction : ", coda.options("frac2"), "\n", sep="")
+ cat("Number of bins : ", coda.options("geweke.nbin"), "\n", sep="")
+ cat("\n")
+
+ cat("Gelman & Rubin\n")
+ cat("--------------\n")
+ cat("Bin width : ", coda.options("gr.bin"), "\n", sep="")
+ cat("Max number of bins : ", coda.options("gr.max"), "\n", sep="")
+ cat("\n")
+
+ cat("Raftery & Lewis\n")
+ cat("---------------\n")
+ cat("Quantile (q) : ", coda.options("q"), "\n", sep="")
+ cat("Precision (+/- r) : ", coda.options("r"), "\n", sep="")
+ cat("Probability (s) : ", coda.options("s"), "\n", sep="")
+ cat("\n")
+
+ cat("Cross-correlations\n")
+ cat("------------------\n")
+ cat("Combine chains : ", coda.options("combine.corr"), "\n", sep="")
+ cat("\n")
+ }
+ invisible()
+}
+
+"read.coda.interactive" <- function ()
+{
+ repeat {
+ cat("Enter CODA index file name\n")
+ cat("(or a blank line to exit)\n")
+ index.file <- scan(what = character(), sep = "\n", strip.white = TRUE,
+ quiet = TRUE, nlines=1)
+ if (length(index.file) == 0)
+ return()
+ cat("Enter CODA output file names, separated by return key\n")
+ cat("(leave a blank line when you have finished)\n")
+ output.files <- scan(what = character(), sep = "\n", strip.white = TRUE,
+ quiet = TRUE)
+ all.files <- c(index.file, output.files)
+ if (any(!file.exists(all.files))) {
+ cat("The following files were not found:\n")
+ cat(paste(all.files[!file.exists(all.files)],
+ collapse = "\n"), "\n\n")
+ }
+ else break
+ }
+ nfiles <- length(output.files)
+ chains <- vector("list", nfiles)
+ names(chains) <- output.files
+ for (i in 1:nfiles) chains[[i]] <- read.coda(output.files[i], index.file)
+ return(mcmc.list(chains))
+}
+
+"tidy.up" <- function ()
+{
+ cat("\nQuitting codamenu ...\n")
+ if (!coda.options("data.saved")) {
+ ans <- read.yesno("Do you want to save the mcmc output", default=FALSE)
+ if (ans == TRUE) {
+ cat("Enter name you want to call this object file:\n")
+ fname <- scan(what = character(), nmax = 1, strip.white = TRUE)
+ assign(fname, coda.dat, pos=1)
+ }
+ }
+ coda.objects <- c("coda.dat", "work.dat")
+ for (i in coda.objects) {
+ if (exists(i))
+ rm(list=i, pos = 1)
+ }
+}
+
+"codamenu.ps" <- function ()
+{
+ if (names(dev.cur()) == "postscript") {
+ dev.off()
+ }
+ else {
+ cat("\nSave plots as a postscript file (y/N) ?\n")
+ ans <- readline()
+ if (length(ans) == 0)
+ ans <- "n"
+ if (ans == "Y" | ans == "y") {
+ repeat {
+ mssg <- "Enter name you want to call this postscript file"
+ ps.name <- read.and.check(mssg, what = character(),
+ default = "Rplots.ps")
+ if (file.exists(ps.name)) {
+ pick <- menu(title = "File exists", choices = c("overwrite",
+ "choose another file name"))
+ if (pick == 1)
+ break
+ }
+ }
+ postscript(file = ps.name)
+ }
+ }
+ return(dev.cur())
+}
+
+"codamenu.output.header" <- function (title)
+{
+ ##
+ ## A short header: common to most codamenu output
+ ##
+ cat("\n", title, sep = "")
+ cat("\n", paste(rep("=", nchar(title)), collapse = ""), "\n\n",
+ sep = "")
+ cat("Iterations used = ", start(work.dat), ":", end(work.dat),
+ "\n", sep = "")
+ cat("Thinning interval =", thin(work.dat), "\n")
+ cat("Sample size per chain =", niter(work.dat), "\n\n")
+ invisible()
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/R/cumuplot.R b/R/cumuplot.R
new file mode 100644
index 0000000..4b16ff3
--- /dev/null
+++ b/R/cumuplot.R
@@ -0,0 +1,40 @@
+cumuplot <- function(x, probs=c(0.025,0.5,0.975), ylab="", lty=c(2,1),
+ lwd=c(1,2), type="l", ask=TRUE, auto.layout=TRUE,
+ col=1, ...)
+{
+ cquantile <- function(z, probs)
+ {
+ ## Calculates cumulative quantile of a vector
+ cquant <- matrix(0, nrow=length(z), length(probs))
+ for(i in seq(along=z)) # for loop proved faster than apply here
+ cquant[i,] <- quantile(z[1:i], probs=probs, names=FALSE)
+ cquant <- as.data.frame(cquant)
+ names(cquant) <- paste(formatC(100*probs,format="fg",wid=1,digits=7),
+ "%", sep="") # just like quantile.default
+ return(cquant)
+ }
+
+ oldpar <- par(ask = ask)
+ on.exit(par(oldpar))
+ if (auto.layout) {
+ oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x),
+ Nparms = nvar(x)))
+ oldpar <- c(oldpar, par(ask = ask))
+ }
+
+ if (!is.mcmc.list(x))
+ x <- mcmc.list(as.mcmc(x))
+
+ Iterations <- time(x)
+ for (i in 1:nchain(x)) {
+ for (j in 1:nvar(x)) {
+ Y <- cquantile(as.matrix(x[[i]])[,j], probs=probs)
+ matplot(Iterations, Y, ylab=ylab, lty=lty, lwd=lwd, type=type,
+ col=col, ...)
+ title(paste(varnames(x)[j], ifelse(is.null(chanames(x)),
+ "", ":"), chanames(x)[i], sep = ""))
+ }
+ }
+}
+
+
diff --git a/R/gelman.R b/R/gelman.R
new file mode 100644
index 0000000..fead51e
--- /dev/null
+++ b/R/gelman.R
@@ -0,0 +1,236 @@
+"gelman.diag" <- function (x, confidence = 0.95, transform = FALSE)
+ ## Gelman and Rubin's diagnostic
+ ## Gelman, A. and Rubin, D (1992). Inference from iterative simulation
+ ## using multiple sequences. Statistical Science, 7, 457-551.
+ ##
+ ## Correction and Multivariate generalization:
+ ## Brooks, S.P. and Gelman, A. (1997) General methods for monitoring
+ ## convergence of iterative simulations. Journal of Computational and
+ ## Graphical Statistics, 7, 434-455.
+
+{
+ x <- as.mcmc.list(x)
+ if (nchain(x) < 2)
+ stop("You need at least two chains")
+ if (start(x) < end(x)/2)
+ x <- window(x, start = end(x)/2 + 1)
+ Niter <- niter(x)
+ Nchain <- nchain(x)
+ Nvar <- nvar(x)
+ xnames <- varnames(x)
+
+ if(transform)
+ x <- gelman.transform(x)
+ ##
+ ## Estimate mean within-chain variance (W) and between-chain variance
+ ## (B/Niter), and calculate sampling variances and covariance of the
+ ## estimates (varW, varB, covWB)
+ ##
+ ## Multivariate (upper case)
+ x <- lapply(x, as.matrix)
+ S2 <- array(sapply(x, var, simplify=TRUE), dim=c(Nvar,Nvar,Nchain))
+ W <- apply(S2, c(1,2), mean)
+ xbar <- matrix(sapply(x, apply, 2, mean, simplify=TRUE), nrow=Nvar,
+ ncol=Nchain)
+ B <- Niter * var(t(xbar))
+
+ if(Nvar > 1) {
+ emax <- eigen(qr.solve(W,B), symmetric=FALSE, only.values=TRUE)$values[1]
+ mpsrf <- sqrt( (1 - 1/Niter) + (1 + 1/Nvar) * emax/Niter )
+ }
+ else
+ mpsrf <- NULL
+ ## Univariate (lower case)
+ w <- diag(W)
+ b <- diag(B)
+
+
+ s2 <- matrix(apply(S2, 3, diag), nrow=Nvar, ncol=Nchain)
+ muhat <- apply(xbar,1,mean)
+ var.w <- apply(s2, 1, var)/Nchain
+ var.b <- (2 * b^2)/(Nchain - 1)
+ cov.wb <- (Niter/Nchain) * diag(var(t(s2), t(xbar^2)) -
+ 2 * muhat * var(t(s2), t(xbar)))
+ ##
+ ## Posterior interval combines all uncertainties in a t interval with
+ ## center muhat, scale sqrt(V), and df.V degrees of freedom.
+ ##
+ V <- (Niter - 1) * w / Niter + (1 + 1/Nchain) * b/ Niter
+ var.V <- ((Niter - 1)^2 * var.w + (1 + 1/Nchain)^2 *
+ var.b + 2 * (Niter - 1) * (1 + 1/Nchain) * cov.wb)/Niter^2
+ df.V <- (2 * V^2)/var.V
+ ##
+ ## Potential scale reduction factor (that would be achieved by
+ ## continuing simulations forever) is estimated by
+ ## R = sqrt(V/W) * df.adj
+ ## where df.adj is a degrees of freedom adjustment for the width
+ ## of the t-interval.
+ ##
+ ## To calculate upper confidence interval we divide R2 = R^2 into two
+ ## parts, fixed and random. The upper limit of the random part is
+ ## calculated assuming that B/W has an F distribution.
+ ##
+ df.adj <- (df.V + 3)/(df.V + 1)
+ B.df <- Nchain - 1
+ W.df <- (2 * w^2)/var.w
+ R2.fixed <- (Niter - 1)/Niter
+ R2.random <- (1 + 1/Nchain) * (1/Niter) * (b/w)
+ R2.estimate <- R2.fixed + R2.random
+ R2.upper <- R2.fixed + qf((1 + confidence)/2, B.df, W.df) * R2.random
+ psrf <- cbind(sqrt(df.adj * R2.estimate), sqrt(df.adj * R2.upper))
+ dimnames(psrf) <- list(xnames,
+ c("Point est.",
+ paste(50 * (1+confidence), "% quantile", sep = ""))
+ )
+
+
+ out <- list(psrf = psrf, mpsrf=mpsrf)
+ class(out) <- "gelman.diag"
+ out
+}
+
+"gelman.transform" <- function(x)
+ ## Gelman and Rubin diagnostic assumes a normal distribution. To
+ ## improve the normal approximation, variables on [0, Inf) are log
+ ## transformed, and variables on [0,1] are logit-transformed.
+{
+ if (nvar(x) == 1) {
+ z <- data.frame(lapply(x, unclass))
+ if (min(z) > 0) {
+ y <- if(max(z) < 1)
+ log(z/(1-z))
+ else log(z)
+ for (j in 1:nchain(x)) x[[j]] <- y[,j]
+ }
+ }
+ else for (i in 1:nvar(x)) {
+ z <- data.frame(lapply(x[, i], unclass))
+ if (min(z) > 0) {
+ y <- if (max(z) < 1)
+ log(z/(1 - z))
+ else log(z)
+ for (j in 1:nchain(x)) x[[j]][, i] <- y[, j]
+ }
+ }
+ return(x)
+}
+
+"gelman.mv.diag" <- function (x, confidence = 0.95, transform = FALSE)
+{
+ s2 <- sapply(x, var, simplify=TRUE)
+ W <- matrix(apply(s2, 1, mean), nvar(x), nvar(x))
+ xbar <- sapply(x, apply, 2, mean, simplify=TRUE)
+ B <- niter(x) * var(t(xbar))
+ emax <- eigen(qr.solve(W,B), symmetric=FALSE, only.values=TRUE)$values[1]
+ mpsrf <- sqrt( (1 - 1/niter(x)) + (1 + 1/nvar(x)) * emax )
+ return(mpsrf)
+}
+
+
+"print.gelman.diag" <-
+ function (x, digits = 3, ...)
+{
+ cat("Potential scale reduction factors:\n\n")
+ print.default(x$psrf, digits = digits, ...)
+ if(!is.null(x$mpsrf)) {
+ cat("\nMultivariate psrf\n\n")
+ cat(format(x$mpsrf,digits = digits))
+ }
+ cat("\n")
+}
+
+"gelman.plot" <-
+ function (x, bin.width = 10, max.bins = 50, confidence = 0.95,
+ transform = FALSE, auto.layout = TRUE, ask = TRUE,
+ col = 1:2, lty = 1:2, xlab = "last iteration in chain",
+ ylab = "shrink factor", type = "l", ...)
+{
+ x <- as.mcmc.list(x)
+ oldpar <- NULL
+ on.exit(par(oldpar))
+ if (auto.layout)
+ oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x), Nparms = nvar(x)))
+ y <- gelman.preplot(x, bin.width = bin.width, max.bins = max.bins,
+ confidence = confidence)
+ all.na <- apply(is.na(y$shrink[, , 1, drop = FALSE]), 2, all)
+ if (!any(all.na))
+ for (j in 1:nvar(x)) {
+ matplot(y$last.iter, y$shrink[, j, ], col = col,
+ lty = lty, xlab = xlab, ylab = ylab, type = type,
+ ...)
+ abline(h = 1)
+ ymax <- max(c(1, y$shrink[, j, ]), na.rm = TRUE)
+ leg <- dimnames(y$shrink)[[3]]
+ xmax <- max(y$last.iter)
+ legend(xmax, ymax, legend = leg, lty = lty, bty = "n",
+ col = col, xjust = 1, yjust = 1)
+ title(main = varnames(x)[j])
+ if (j==1)
+ oldpar <- c(oldpar, par(ask = ask))
+ }
+ return(invisible(y))
+}
+
+"gelman.preplot" <-
+ function (x, confidence = 0.95, transform = FALSE, bin.width = 10,
+ max.bins = 50)
+{
+ x <- as.mcmc.list(x)
+ if (niter(x) <= 50)
+ stop("Less than 50 iterations in chain")
+ nbin <- min(floor((niter(x) - 50)/thin(x)), max.bins)
+ binw <- floor((niter(x) - 50)/nbin)
+ last.iter <- c(seq(from = start(x) + 50 * thin(x), by = binw *
+ thin(x), length = nbin), end(x))
+ shrink <- array(dim = c(nbin + 1, nvar(x), 2))
+ dimnames(shrink) <- list(last.iter, varnames(x),
+ c("median", paste(50 * (confidence + 1), "%",
+ sep = ""))
+ )
+ for (i in 1:(nbin + 1)) {
+ shrink[i, , ] <- gelman.diag(window(x, end = last.iter[i]),
+ confidence = confidence,
+ transform = transform)$psrf
+ }
+ all.na <- apply(is.na(shrink[, , 1, drop = FALSE]), 2, all)
+ if (any(all.na)) {
+ cat("\n******* Error: *******\n")
+ cat("Cannot compute Gelman & Rubin's diagnostic for any chain \n")
+ cat("segments for variables", varnames(x)[all.na], "\n")
+ cat("This indicates convergence failure\n")
+ }
+ return(list(shrink = shrink, last.iter = last.iter))
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/R/geweke.R b/R/geweke.R
new file mode 100644
index 0000000..d5f52fb
--- /dev/null
+++ b/R/geweke.R
@@ -0,0 +1,81 @@
+"geweke.diag" <-
+ function (x, frac1 = 0.1, frac2 = 0.5)
+ ##
+{
+ if (is.mcmc.list(x))
+ return(lapply(x, geweke.diag, frac1, frac2))
+ x <- as.mcmc(x)
+ xstart <- c(start(x), end(x) - frac2 * (end(x) - start(x)))
+ xend <- c(start(x) + frac1 * (end(x) - start(x)), end(x))
+ y.variance <- y.mean <- vector("list", 2)
+ for (i in 1:2) {
+ y <- window(x, start = xstart[i], end = xend[i])
+ y.mean[[i]] <- apply(as.matrix(y), 2, mean)
+ y.variance[[i]] <- spectrum0(y)$spec/niter(y)
+ }
+ z <- (y.mean[[1]] - y.mean[[2]])/sqrt(y.variance[[1]] + y.variance[[2]])
+ out <- list(z = z, frac = c(frac1, frac2))
+ class(out) <- "geweke.diag"
+ return(out)
+}
+
+"geweke.plot" <-
+ function (x, frac1 = 0.1, frac2 = 0.5, nbins = 20,
+ pvalue = 0.05, auto.layout = TRUE, ask = TRUE, ...)
+{
+ x <- as.mcmc.list(x)
+ oldpar <- NULL
+ on.exit(par(oldpar))
+ if (auto.layout)
+ oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x),
+ Nparms = nvar(x)))
+ ystart <- seq(from = start(x), to = (start(x) + end(x))/2, length = nbins)
+ gcd <- array(dim = c(length(ystart), nvar(x), nchain(x)),
+ dimnames = c(ystart, varnames(x), chanames(x)))
+ for (n in 1:length(ystart)) {
+ geweke.out <- geweke.diag(window(x, start = ystart[n]),
+ frac1 = frac1, frac2 = frac2)
+ for (k in 1:nchain(x)) gcd[n, , k] <- geweke.out[[k]]$z
+ }
+ climit <- qnorm(1 - pvalue/2)
+ for (k in 1:nchain(x)) for (j in 1:nvar(x)) {
+ ylimit <- max(c(climit, abs(gcd[, j, k])))
+ plot(ystart, gcd[, j, k], type = "p", xlab = "First iteration in segment",
+ ylab = "Z-score", pch = 4, ylim = c(-ylimit, ylimit),
+ ...)
+ abline(h = c(climit, -climit), lty = 2)
+ if (nchain(x) > 1) {
+ title(main = paste(varnames(x, allow.null = FALSE)[j],
+ " (", chanames(x, allow.null = FALSE)[k], ")",
+ sep = ""))
+ }
+ else {
+ title(main = paste(varnames(x, allow.null = FALSE)[j],
+ sep = ""))
+ }
+ if (k==1 && j==1)
+ oldpar <- c(oldpar, par(ask = ask))
+ }
+ invisible(list(start.iter = ystart, z = gcd))
+}
+
+"print.geweke.diag" <- function (x, digits = min(4, .Options$digits), ...)
+ ## Print method for output from geweke.diag
+{
+ cat("\nFraction in 1st window =", x$frac[1])
+ cat("\nFraction in 2nd window =", x$frac[2], "\n\n")
+ print.default(x$z, digits = digits, ...)
+ cat("\n")
+ invisible(x)
+}
+
+
+
+
+
+
+
+
+
+
+
diff --git a/R/heidel.R b/R/heidel.R
new file mode 100644
index 0000000..4282f2d
--- /dev/null
+++ b/R/heidel.R
@@ -0,0 +1,200 @@
+"heidel.diag" <- function (x, eps = 0.1, pvalue=0.05)
+{
+ if (is.mcmc.list(x))
+ return(lapply(x, heidel.diag, eps))
+ x <- as.mcmc(as.matrix(x))
+ HW.mat0 <- matrix(0, ncol = 6, nrow = nvar(x))
+ dimnames(HW.mat0) <- list(varnames(x),
+ c("stest", "start", "pvalue", "htest",
+ "mean", "halfwidth"))
+ HW.mat <- HW.mat0
+ for (j in 1:nvar(x)) {
+ start.vec <- seq(from=start(x), to = end(x)/2, by=niter(x)/10)
+ Y <- x[, j, drop = TRUE]
+ n1 <- length(Y)
+ ## Schruben's test for convergence, applied sequentially
+ ##
+ S0 <- spectrum0(window(Y, start=end(Y)/2))$spec
+ converged <- FALSE
+ for (i in seq(along = start.vec)) {
+ Y <- window(Y, start = start.vec[i])
+ n <- niter(Y)
+ ybar <- mean(Y)
+ B <- cumsum(Y) - ybar * (1:n)
+ Bsq <- (B * B)/(n * S0)
+ I <- sum(Bsq)/n
+ if(converged <- !is.na(I) && pcramer(I) < 1 - pvalue)
+ break
+ }
+ ## Recalculate S0 using section of chain that passed convergence test
+ S0ci <- spectrum0(Y)$spec
+ halfwidth <- 1.96 * sqrt(S0ci/n)
+ passed.hw <- !is.na(halfwidth) & (abs(halfwidth/ybar) <= eps)
+ if (!converged || is.na(I) || is.na(halfwidth)) {
+ nstart <- NA
+ passed.hw <- NA
+ halfwidth <- NA
+ ybar <- NA
+ }
+ else {
+ nstart <- start(Y)
+ }
+ HW.mat[j, ] <- c(converged, nstart, 1 - pcramer(I),
+ passed.hw, ybar, halfwidth)
+ }
+ class(HW.mat) <- "heidel.diag"
+ return(HW.mat)
+}
+
+"print.heidel.diag" <-
+ function (x, digits = 3, ...)
+{
+ HW.title <- matrix(c("Stationarity", "test", "start", "iteration",
+ "p-value", "",
+ "Halfwidth", "test", "Mean", "", "Halfwidth", ""),
+ nrow = 2)
+ y <- matrix("", nrow = nrow(x), ncol = 6)
+ for (j in 1:ncol(y)) {
+ y[, j] <- format(x[, j], digits = digits)
+ }
+ y[, c(1, 4)] <- ifelse(x[, c(1, 4)], "passed", "failed")
+ y <- rbind(HW.title, y)
+ vnames <- if (is.null(rownames(x)))
+ paste("[,", 1:nrow(x), "]", sep = "")
+ else rownames(x)
+ dimnames(y) <- list(c("", "", vnames), rep("", 6))
+ print.default(y[, 1:3], quote = FALSE, ...)
+ print.default(y[, 4:6], quote = FALSE, ...)
+ invisible(x)
+}
+
+"spectrum0.ar" <- function(x)
+{
+ x <- as.matrix(x)
+ v0 <- order <- numeric(ncol(x))
+ names(v0) <- names(order) <- colnames(x)
+ z <- 1:nrow(x)
+ for (i in 1:ncol(x))
+ {
+ lm.out <- lm(x[,i] ~ z)
+ if (identical(all.equal(var(residuals(lm.out)), 0), TRUE)) {
+ v0[i] <- 0
+ order[i] <- 0
+ }
+ else {
+ ar.out <- ar(x[,i], aic=TRUE)
+ v0[i] <- ar.out$var.pred/(1 - sum(ar.out$ar))^2
+ order[i] <- ar.out$order
+ }
+ }
+ return(list(spec=v0, order=order))
+}
+
+effectiveSize <- function(x)
+{
+ if (is.mcmc.list(x))
+ {
+ ans <- vector("list", nchain(x))
+ for (i in 1:nchain(x))
+ ans[[i]] <- effectiveSize(x[[i]])
+ }
+ else
+ {
+ x <- as.mcmc(x)
+ x <- as.matrix(x)
+ spec <- spectrum0.ar(x)$spec
+ ans <- ifelse(spec==0, 0, nrow(x) * apply(x, 2, var)/spec)
+ }
+ return(ans)
+}
+
+"spectrum0" <- function(x, max.freq=0.5, order=1, max.length=200)
+{
+ x <- as.matrix(x)
+ if (!is.null(max.length) && nrow(x) > max.length) {
+ batch.size <- ceiling(nrow(x)/max.length)
+ x <- aggregate(ts(x, frequency=batch.size), nfreq = 1, FUN=mean)
+ }
+ else {
+ batch.size <- 1
+ }
+
+ out <- do.spectrum0(x, max.freq=max.freq, order=order)
+ out$spec <- out$spec * batch.size
+ return(out)
+}
+
+"do.spectrum0" <- function(x, max.freq=0.5, order=1)
+{
+ ## Estimate spectral density of time series x at frequency 0.
+ ## spectrum0(x)/length(x) estimates the variance of mean(x)
+ ##
+ ## NB We do NOT use the same definition of spectral density
+ ## as in spec.pgram.
+ ##
+ fmla <- switch(order+1,
+ spec ~ one,
+ spec ~ f1,
+ spec ~ f1 + f2)
+ if(is.null(fmla))
+ stop("invalid order")
+
+ N <- nrow(x)
+ Nfreq <- floor(N/2)
+ freq <- seq(from = 1/N, by = 1/N, length = Nfreq)
+ f1 <- sqrt(3) * (4 * freq - 1)
+ f2 <- sqrt(5) * (24 * freq^2 - 12 * freq + 1)
+ v0 <- numeric(ncol(x))
+ for(i in 1:ncol(x)) {
+ y <- x[,i]
+ if (var(y) == 0) {
+ v0[i] <- 0
+ }
+ else {
+ yfft <- fft(y)
+ spec <- Re(yfft * Conj(yfft))/ N
+ spec.data <- data.frame(one = rep(1, Nfreq), f1=f1, f2=f2,
+ spec = spec[1 + (1:Nfreq)],
+ inset = I(freq<=max.freq))
+
+ glm.out <- glm(fmla, family=Gamma(link="log"), data=spec.data,
+ subset=inset)
+ v0[i] <- predict(glm.out, type="response",
+ newdata=data.frame(spec=0,one=1,f1=-sqrt(3),f2=sqrt(5)))
+ }
+ }
+ return(list(spec=v0))
+}
+
+"pcramer" <- function (q, eps=1.0e-5)
+{
+ ## Distribution function of the Cramer-von Mises statistic
+ ##
+ log.eps <- log(eps)
+ y <- matrix(0, nrow=4, ncol=length(q))
+ for(k in 0:3) {
+ z <- gamma(k + 0.5) * sqrt(4*k + 1)/(gamma(k+1) * pi^(3/2) * sqrt(q))
+ u <- (4*k + 1)^2/(16*q)
+ y[k+1,] <- ifelse(u > -log.eps, 0, z * exp(-u) * besselK(x = u, nu=1/4))
+ }
+ return(apply(y,2,sum))
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/R/jags.R b/R/jags.R
new file mode 100644
index 0000000..412df9d
--- /dev/null
+++ b/R/jags.R
@@ -0,0 +1,22 @@
+"read.jags" <- function (file = "jags.out", start, end, thin, quiet=FALSE)
+{
+ read.coda(file, start, end, thin, quiet)
+}
+
+bugs2jags <- function(infile, outfile)
+{
+ ## Convert S-style data for WinBUGS into the R dump format
+ ## used by JAGS.
+ bugs.dat <- dget(infile)
+ for (bugs.variable.name in names(bugs.dat)) {
+ if(!is.null(dim(bugs.dat[[bugs.variable.name]]))) {
+ ## Manually reverse order of dimensions of arrays
+ dim(bugs.dat[[bugs.variable.name]]) <-
+ rev(dim(bugs.dat[[bugs.variable.name]]))
+ ## Then transpose
+ bugs.dat[[bugs.variable.name]] <- aperm(bugs.dat[[bugs.variable.name]])
+ }
+ assign(bugs.variable.name, bugs.dat[[bugs.variable.name]])
+ }
+ dump(names(bugs.dat), file=outfile)
+}
diff --git a/R/mcextractor.R b/R/mcextractor.R
new file mode 100644
index 0000000..47c52c8
--- /dev/null
+++ b/R/mcextractor.R
@@ -0,0 +1,94 @@
+"chanames" <-
+function (x, allow.null = TRUE)
+{
+ if (is.mcmc.list(x)) {
+ if (is.null(names(x)))
+ if (allow.null)
+ NULL
+ else paste("chain", 1:length(x), sep = "")
+ else names(x)
+ }
+ else NULL
+}
+
+"chanames<-" <-
+function (x, value)
+{
+ if (is.mcmc.list(x))
+ names(x) <- value
+ else stop("Not an mcmc.list object")
+ x
+}
+
+"varnames" <-
+function (x, allow.null = TRUE)
+{
+ if (!is.mcmc(x) && !is.mcmc.list(x))
+ return(NULL)
+ y <- if (is.mcmc(x))
+ dimnames(x)[[2]]
+ else if (is.mcmc.list(x))
+ dimnames(x[[1]])[[2]]
+ if (is.null(y) && !allow.null)
+ y <- paste("var", 1:nvar(x), sep = "")
+ return(y)
+}
+
+"varnames<-" <-
+function (x, value)
+{
+ if (is.mcmc(x)) {
+ dimnames(x)[[2]] <- value
+ }
+ else if (is.mcmc.list(x))
+ for (i in 1:nchain(x)) varnames(x[[i]]) <- value
+ else stop("Not an mcmc or mcmc.list object")
+ x
+}
+
+"nchain" <-
+function (x)
+{
+ if (is.mcmc(x))
+ 1
+ else if (is.mcmc.list(x))
+ length(x)
+ else NULL
+}
+
+"nvar" <-
+function (x)
+{
+
+ if (is.mcmc(x)) {
+ if (is.matrix(x)) ncol(x) else 1
+ }
+ else if (is.mcmc.list(x)) {
+ if (is.matrix(x[[1]])) ncol(x[[1]]) else 1
+ }
+ else NULL
+}
+
+"niter" <-
+function (x)
+{
+ if (is.mcmc(x)) {
+ if (is.matrix(x)) nrow(x) else length(x)
+ }
+ else if (is.mcmc.list(x)) {
+ if (is.matrix(x[[1]])) nrow(x[[1]]) else length(x[[1]])
+ }
+ else NULL
+}
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/R/mcmc.R b/R/mcmc.R
new file mode 100644
index 0000000..026c1de
--- /dev/null
+++ b/R/mcmc.R
@@ -0,0 +1,310 @@
+"[.mcmc" <- function (x, i, j, drop = missing(i))
+{
+ xstart <- start(x)
+ xthin <- thin(x)
+ y <- NextMethod("[")
+ if (length(y) == 0 || is.null(y))
+ return(y)
+ if (missing(i))
+ return(mcmc(y, start = xstart, thin = xthin))
+ else
+ return(y)
+}
+
+"as.mcmc" <- function (x)
+ UseMethod("as.mcmc")
+
+"as.mcmc.default" <- function (x)
+ if (is.mcmc(x)) x else mcmc(x)
+
+"as.ts.mcmc" <- function (x)
+{
+ x <- as.mcmc(x)
+ y <- ts(x, start = start(x), end = end(x), deltat = thin(x))
+ attr(y, "mcpar") <- NULL
+ return(y)
+}
+
+"start.mcmc" <- function (x, ...)
+{
+ mcpar(as.mcmc(x))[1]
+}
+
+"end.mcmc" <- function (x, ...)
+{
+ mcpar(as.mcmc(x))[2]
+}
+
+"frequency.mcmc" <- function (x, ...)
+{
+ 1/thin.mcmc(x)
+}
+
+"thin.mcmc" <- function (x, ...)
+{
+ mcpar(as.mcmc(x))[3]
+}
+
+"is.mcmc" <- function (x)
+{
+ if (inherits(x, "mcmc"))
+ if (length(dim(x)) == 3)
+ stop("Obsolete mcmc object\nUpdate with a command like\nx <- mcmcUpgrade(x)")
+ else TRUE
+ else FALSE
+}
+
+"mcmc" <- function (data = NA, start = 1, end = numeric(0), thin = 1)
+{
+ if (is.matrix(data)) {
+ niter <- nrow(data)
+ nvar <- ncol(data)
+ }
+ else if (is.data.frame(data)) {
+ if (!all(sapply(data, is.numeric))) {
+ stop ("Data frame contains non-numeric values")
+ }
+ data <- as.matrix(data)
+ niter <- nrow(data)
+ nvar <- ncol(data)
+ }
+ else {
+ niter <- length(data)
+ nvar <- 1
+ }
+ thin <- round(thin)
+ if (length(start) > 1)
+ stop("Invalid start")
+ if (length(end) > 1)
+ stop("Invalid end")
+ if (length(thin) != 1)
+ stop("Invalid thin")
+ if (missing(end))
+ end <- start + (niter - 1) * thin
+ else if (missing(start))
+ start <- end - (niter - 1) * thin
+ nobs <- floor((end - start)/thin + 1.0) ### patch
+ if (niter < nobs)
+ stop("Start, end and thin incompatible with data")
+ else {
+ end <- start + thin * (nobs - 1)
+ if (nobs < niter)
+ data <- data[1:nobs, , drop = FALSE]
+ }
+ attr(data, "mcpar") <- c(start, end, thin)
+ attr(data, "class") <- "mcmc"
+ data
+}
+
+"print.mcmc" <- function (x, ...)
+{
+ x.orig <- x
+ cat("Markov Chain Monte Carlo (MCMC) output:\nStart =", start(x),
+ "\nEnd =", end(x), "\nThinning interval =", thin(x), "\n")
+ attr(x, "mcpar") <- NULL
+ attr(x, "class") <- NULL
+ NextMethod("print", ...)
+ invisible(x.orig)
+}
+
+
+"as.matrix.mcmc" <- function (x, iters = FALSE)
+{
+ y <- matrix(nrow = niter(x), ncol = nvar(x) + iters)
+ var.cols <- iters + 1:nvar(x)
+ if (iters)
+ y[, 1] <- as.vector(time(x))
+ y[, var.cols] <- x
+ rownames <- character(ncol(y))
+ if (iters)
+ rownames[1] <- "ITER"
+ rownames[var.cols] <- varnames(x, allow.null = FALSE)
+ dimnames(y) <- list(NULL, rownames)
+ return(y)
+}
+
+"time.mcmc" <- function (x, ...)
+{
+ x <- as.mcmc(x)
+ ts(seq(from = start(x), to = end(x), by = thin(x)), start = start(x),
+ end = end(x), deltat = thin(x))
+}
+
+"window.mcmc" <- function (x, start, end, thin, ...)
+{
+ ts.eps <- getOption("ts.eps")
+ xmcpar <- mcpar(x)
+ xstart <- xmcpar[1]
+ xend <- xmcpar[2]
+ xthin <- xmcpar[3]
+ if (missing(thin))
+ thin <- xthin
+ else if (thin%%xthin != 0) {
+ thin <- xthin
+ warning("Thin value not changed")
+ }
+ xtime <- as.vector(time(x))
+ if (missing(start))
+ start <- xstart
+ else if (length(start) != 1)
+ stop("bad value for start")
+ else if (start < xstart) {
+ start <- xstart
+ warning("start value not changed")
+ }
+ if (missing(end))
+ end <- xend
+ else if (length(end) != 1)
+ stop("bad value for end")
+ else if (end > xend) {
+ end <- xend
+ warning("end value not changed")
+ }
+ if (start > end)
+ stop("start cannot be after end")
+ if (all(abs(xtime - start) > abs(start) * ts.eps)) {
+ start <- xtime[(xtime > start) & ((start + xthin) > xtime)]
+ }
+ if (all(abs(end - xtime) > abs(end) * ts.eps)) {
+ end <- xtime[(xtime < end) & ((end - xthin) < xtime)]
+ }
+ use <- 1:niter(x)
+ use <- use[use >= trunc((start - xstart)/xthin + 1.5) &
+ use <= trunc((end - xstart)/xthin + 1.5) &
+ (use - trunc((start- xstart)/xthin + 1.5))%%(thin%/%xthin) == 0]
+ y <- if (is.matrix(x))
+ x[use, , drop = FALSE]
+ else x[use]
+ return(mcmc(y, start=start, end=end, thin=thin))
+}
+
+"mcpar" <- function (x)
+{
+ attr(x, "mcpar")
+}
+
+"mcmcUpgrade" <- function (x)
+{
+ if (inherits(x, "mcmc")) {
+ if (length(dim(x)) == 3) {
+ nchain <- dim(x)[3]
+ xtspar <- attr(x, "tspar")
+ xstart <- xtspar[1]
+ xend <- xtspar[2]
+ xthin <- xtspar[3]
+ out <- vector("list", nchain)
+ for (i in 1:nchain) {
+ y <- unclass(x)[, , 1, drop = TRUE]
+ attr(y, "title") <- NULL
+ attr(y, "tspar") <- NULL
+ out[[i]] <- mcmc(y, start = xstart, end = xend,
+ thin = xthin)
+ }
+ if (nchain == 1)
+ return(out[[1]])
+ else return(mcmc.list(out))
+ }
+ else return(x)
+ }
+ else stop("Can't upgrade")
+}
+
+"thin" <-
+function (x, ...)
+ UseMethod("thin")
+
+"set.mfrow" <-
+function (Nchains = 1, Nparms = 1, nplots = 1, sepplot = FALSE)
+{
+ ## Set up dimensions of graphics window:
+ ## If only density plots OR trace plots are requested, dimensions are:
+ ## 1 x 1 if Nparms = 1
+ ## 1 X 2 if Nparms = 2
+ ## 2 X 2 if Nparms = 3 or 4
+ ## 3 X 2 if Nparms = 5 or 6 or 10 - 12
+ ## 3 X 3 if Nparms = 7 - 9 or >= 13
+ ## If both density plots AND trace plots are requested, dimensions are:
+ ## 1 x 2 if Nparms = 1
+ ## 2 X 2 if Nparms = 2
+ ## 3 X 2 if Nparms = 3, 5, 6, 10, 11, or 12
+ ## 4 x 2 if Nparms otherwise
+ ## If separate plots are requested for each chain, dimensions are:
+ ## 1 x 2 if Nparms = 1 & Nchains = 2
+ ## 2 X 2 if Nparms = 2 & Nchains = 2 OR Nparms = 1 & Nchains = 3 or 4
+ ## 3 x 2 if Nparms = 3 or >= 5 & Nchains = 2
+ ## OR Nchains = 5 or 6 or 10 - 12 (and any Nparms)
+ ## 2 x 3 if Nparms = 2 or 4 & Nchains = 3
+ ## 4 x 2 if Nparms = 4 & Nchains = 2
+ ## OR Nchains = 4 & Nparms > 1
+ ## 3 x 3 if Nparms = 3 or >= 5 & Nchains = 3
+ ## OR Nchains = 7 - 9 or >= 13 (and any Nparms)
+ mfrow <- if (sepplot && Nchains > 1 && nplots == 1) {
+ ## Separate plots per chain
+ ## Only one plot per variable
+ if (Nchains == 2) {
+ switch(min(Nparms, 5),
+ c(1,2),
+ c(2,2),
+ c(3,2),
+ c(4,2),
+ c(3,2))
+ }
+ else if (Nchains == 3) {
+ switch(min(Nparms, 5),
+ c(2,2),
+ c(2,3),
+ c(3,3),
+ c(2,3),
+ c(3,3))
+ }
+ else if (Nchains == 4) {
+ if (Nparms == 1)
+ c(2,2)
+ else
+ c(4,2)
+ }
+ else if (any(Nchains == c(5,6,10,11,12)))
+ c(3,2)
+ else if (any(Nchains == c(7,8,9)) || Nchains >=13)
+ c(3,3)
+
+ }
+ else {
+ if (nplots==1) {
+ ## One plot per variable
+ mfrow <- switch(min(Nparms,13),
+ c(1,1),
+ c(1,2),
+ c(2,2),
+ c(2,2),
+ c(3,2),
+ c(3,2),
+ c(3,3),
+ c(3,3),
+ c(3,3),
+ c(3,2),
+ c(3,2),
+ c(3,2),
+ c(3,3))
+ }
+ else {
+ ## Two plot per variable
+ ##
+ mfrow <- switch(min(Nparms, 13),
+ c(1,2),
+ c(2,2),
+ c(3,2),
+ c(4,2),
+ c(3,2),
+ c(3,2),
+ c(4,2),
+ c(4,2),
+ c(4,2),
+ c(3,2),
+ c(3,2),
+ c(3,2),
+ c(4,2))
+ }
+ }
+ return(mfrow)
+}
diff --git a/R/mcmclist.R b/R/mcmclist.R
new file mode 100644
index 0000000..b4d1fcd
--- /dev/null
+++ b/R/mcmclist.R
@@ -0,0 +1,191 @@
+"[.mcmc.list" <- function (x, i, j, drop = TRUE)
+{
+ ## Trying to squeeze too much functionality in here
+ ## x[p:q] will subset the list
+ ## x[p,], x[,q], x[p,q] will be recursively applied to
+ ## the elements of the list, even if they are vectors
+ if (nargs() < 3 + !missing(drop)) {
+ ## Subset the list
+ y <- NextMethod("[")
+ }
+ else {
+ ## Subset the elements of the list
+ y <- vector("list", length(x))
+ names(y) <- names(x)
+ for (k in 1:length(y)) {
+ drop1 <- drop | !is.matrix(x[[k]])
+ y[[k]] <- if (missing(i) && missing(j))
+ x[[k]]
+ else if (missing(i))
+ mcmc(as.matrix(x[[k]])[, j, drop = drop1], start(x), end(x), thin(x))
+ else if (missing(j))
+ as.matrix(x[[k]])[i, , drop = drop1]
+ else as.matrix(x[[k]])[i, j, drop = drop1]
+ }
+ }
+ if (is.list(y) && all(sapply(y, is.mcmc, simplify = TRUE)))
+ y <- mcmc.list(y)
+ return(y)
+}
+
+
+"mcmc.list" <- function (...)
+{
+ x <- list(...)
+ if (length(x) == 1 && is.list(x[[1]]))
+ x <- x[[1]]
+ if (!all(unlist(lapply(x, is.mcmc))))
+ stop("Arguments must be mcmc objects")
+ nargs <- length(x)
+ if (nargs >= 2) {
+ xmcpar <- lapply(x, mcpar)
+ if (!all(unlist(lapply(xmcpar, "==", xmcpar[[1]]))))
+ stop("Different start, end or thin values in each chain")
+ xnvar <- lapply(x, nvar)
+ if (!all(unlist(lapply(xnvar, "==", xnvar[[1]]))))
+ stop("Different number of variables in each chain")
+ xvarnames <- lapply(x, varnames, allow.null = FALSE)
+ if (!all(unlist(lapply(xvarnames, "==", xvarnames[[1]]))))
+ stop("Different variable names in each chain")
+ }
+ class(x) <- "mcmc.list"
+ return(x)
+}
+
+"start.mcmc.list" <- function (x, ...)
+{
+ start(x[[1]])
+}
+
+"end.mcmc.list" <- function (x, ...)
+{
+ end(x[[1]])
+}
+
+"thin.mcmc.list" <- function (x, ...)
+{
+ thin(x[[1]])
+}
+
+"is.mcmc.list" <- function (x)
+ inherits(x, "mcmc.list")
+
+"plot.mcmc.list" <-
+ function (x, trace = TRUE, density = TRUE, smooth = TRUE, bwf,
+ auto.layout = TRUE, ask = TRUE, ...)
+{
+ oldpar <- NULL
+ on.exit(par(oldpar))
+ if (auto.layout) {
+ mfrow <- set.mfrow(Nchains = nchain(x), Nparms = nvar(x),
+ nplots = trace + density)
+ oldpar <- par(mfrow = mfrow)
+ }
+ for (i in 1:nvar(x)) {
+ if (trace)
+ traceplot(x[, i, drop = FALSE], smooth = smooth)
+ if (density) {
+ if (missing(bwf))
+ densplot(x[, i, drop = FALSE])
+ else densplot(x[, i, drop = FALSE], bwf = bwf)
+ }
+ if (i==1)
+ oldpar <- c(oldpar, par(ask = ask))
+ }
+}
+
+"summary.mcmc.list" <-
+ function (object, quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975), ...)
+{
+ x <- mcmc.list(object)
+ statnames <- c("Mean", "SD", "Naive SE", "Time-series SE")
+ varstats <- matrix(nrow = nvar(x), ncol = length(statnames),
+ dimnames = list(varnames(x), statnames))
+ xtsvar <- matrix(nrow = nchain(x), ncol = nvar(x))
+ if (is.matrix(x[[1]])) {
+ for (i in 1:nchain(x))
+ for(j in 1:nvar(x))
+ xtsvar[i, j] <- spectrum0(x[[i]][,j])$spec
+ xlong <- do.call("rbind", x)
+ }
+ else {
+ for (i in 1:nchain(x))
+ xtsvar[i, ] <- spectrum0(x[[i]])$spec
+ xlong <- as.matrix(x)
+ }
+
+ xmean <- apply(xlong, 2, mean)
+ xvar <- apply(xlong, 2, var)
+ xtsvar <- apply(xtsvar, 2, mean)
+ varquant <- t(apply(xlong, 2, quantile, quantiles))
+ varstats[, 1] <- xmean
+ varstats[, 2] <- sqrt(xvar)
+ varstats[, 3] <- sqrt(xvar/niter(x))
+ varstats[, 4] <- sqrt(xtsvar/niter(x))
+ varquant <- drop(varquant)
+ varstats <- drop(varstats)
+ out <- list(statistics = varstats, quantiles = varquant,
+ start = start(x), end = end(x), thin = thin(x),
+ nchain = nchain(x))
+ class(out) <- "summary.mcmc"
+ return(out)
+}
+
+"as.matrix.mcmc.list" <-
+ function (x, iters = FALSE, chains = FALSE)
+{
+ x <- mcmc.list(x)
+ y <- matrix(nrow = niter(x) * nchain(x), ncol = nvar(x) +
+ chains + iters)
+ var.cols <- chains + iters + 1:nvar(x)
+ for (i in 1:nchain(x)) {
+ use.rows <- niter(x) * (i - 1) + 1:niter(x)
+ if (chains)
+ y[use.rows, 1] <- i
+ if (iters)
+ y[use.rows, chains + 1] <- as.vector(time(x))
+ y[use.rows, var.cols] <- x[[i]]
+ }
+ rownames <- character(ncol(y))
+ if (chains)
+ rownames[1] <- "CHAIN"
+ if (iters)
+ rownames[1 + chains] <- "ITER"
+ rownames[var.cols] <- varnames(x, allow.null = FALSE)
+ dimnames(y) <- list(NULL, rownames)
+ return(y)
+}
+
+"as.mcmc.mcmc.list" <- function (x)
+{
+ if (nchain(x) == 1)
+ return(x[[1]])
+ else stop("Can't coerce mcmc.list to mcmc object:\n more than 1 chain")
+}
+
+"time.mcmc.list" <- function (x, ...)
+ time(x[[1]])
+
+"window.mcmc.list" <- function (x, ...)
+{
+ structure(lapply(x, window.mcmc, ...), class = "mcmc.list")
+}
+
+"as.mcmc.list" <- function (x, ...)
+ UseMethod("as.mcmc.list")
+
+"as.mcmc.list.default" <- function (x, ...)
+ if (is.mcmc.list(x)) x else mcmc.list(x)
+
+"as.array.mcmc.list" <- function(x, drop=TRUE, ...)
+{
+ y <- array(dim=c(niter(x), nvar(x), nchain(x)),
+ dimnames = list(iter=time(x), var=varnames(x), chain=chanames(x)))
+ for(i in 1:nchain(x))
+ y[,,i] <- x[[i]]
+ if(drop)
+ return(drop(y))
+ else
+ return(y)
+}
+
diff --git a/R/output.R b/R/output.R
new file mode 100644
index 0000000..dea78e6
--- /dev/null
+++ b/R/output.R
@@ -0,0 +1,386 @@
+"autocorr" <-
+function (x, lags = c(0, 1, 5, 10, 50), relative = TRUE)
+{
+ if (relative)
+ lags <- lags * thin(x)
+ else if (any(lags%%thin(x) != 0))
+ stop("Lags do not conform to thinning interval")
+ lags <- lags[lags < niter(x) * thin(x)]
+ if (is.mcmc.list(x))
+ return(lapply(x, autocorr, lags, relative))
+ x <- as.mcmc(x)
+ y <- array(dim = c(length(lags), nvar(x), nvar(x)))
+ dimnames(y) <- list(paste("Lag", lags), varnames(x), varnames(x))
+ acf.out <- acf(as.ts.mcmc(x), lag.max = max(lags), plot = FALSE)$acf
+ y[, , ] <- if (is.array(acf.out))
+ acf.out[lags%/%thin(x) + 1, , ]
+ else acf.out[lags%/%thin(x) + 1]
+ return(y)
+}
+
+"autocorr.plot" <-
+function (x, lag.max, auto.layout = TRUE, ask = TRUE, ...)
+{
+ oldpar <- NULL
+ on.exit(par(oldpar))
+ if (auto.layout)
+ oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x),
+ Nparms = nvar(x)))
+ if (!is.mcmc.list(x))
+ x <- mcmc.list(as.mcmc(x))
+ for (i in 1:nchain(x)) {
+ xacf <- if (missing(lag.max))
+ acf(as.ts.mcmc(x[[i]]), plot = FALSE)
+ else acf(as.ts.mcmc(x[[i]]), lag.max = lag.max, plot = FALSE)
+ for (j in 1:nvar(x)) {
+ plot(xacf$lag[, j, j], xacf$acf[, j, j], type = "h",
+ ylab = "Autocorrelation", xlab = "Lag", ylim = c(-1,
+ 1), ...)
+ title(paste(varnames(x)[j], ifelse(is.null(chanames(x)),
+ "", ":"), chanames(x)[i], sep = ""))
+ if (i==1 && j==1)
+ oldpar <- c(oldpar, par(ask = ask))
+ }
+ }
+ invisible(x)
+}
+
+"crosscorr" <-
+function (x)
+{
+ cor(as.matrix(x))
+}
+
+"crosscorr.plot" <-
+function (x, col = topo.colors(10), ...)
+{
+ Nvar <- nvar(x)
+ pcorr <- crosscorr(x)
+ dens <- ((pcorr + 1) * length(col))%/%2 + (pcorr < 1) + (pcorr <
+ -1)
+ cutoffs <- format(seq(from = 1, to = -1, length = length(col) +
+ 1), digits = 2)
+ leg <- paste("(", cutoffs[-1], ",", cutoffs[-length(cutoffs)],
+ "]", sep = "")
+ oldpar <- NULL
+ on.exit(par(oldpar))
+ oldpar <- c(par(pty = "s", adj = 0.5), oldpar)
+ plot(0, 0, type = "n", xlim = c(0, Nvar), ylim = c(0, Nvar),
+ xlab = "", ylab = "", xaxt = "n", yaxt = "n", ...)
+ axis(1, at = 1:Nvar - 0.5, labels = abbreviate(varnames(x,
+ allow.null = FALSE), minlength = 7))
+ axis(2, at = 1:Nvar - 0.5, labels = abbreviate(varnames(x,
+ allow.null = FALSE), minlength = 7)[Nvar:1])
+ for (cl in 1:Nvar) {
+ for (rw in 1:(Nvar - cl + 1)) polygon(y = c(cl - 1, cl -
+ 1, cl, cl, cl - 1), x = c(rw - 1, rw, rw, rw - 1,
+ rw - 1), col = col[dens[nrow(dens) - cl + 1, rw]])
+ }
+ yval <- seq(from = Nvar/2, to = Nvar, length = length(col) +
+ 1)
+ ydelta <- Nvar/(2 * (length(col) + 1))
+ for (i in 1:length(col)) {
+ polygon(y = c(yval[i], yval[i + 1], yval[i + 1], yval[i],
+ yval[i]), col = col[i], x = c(Nvar - ydelta, Nvar -
+ ydelta, Nvar, Nvar, Nvar - ydelta))
+ }
+ text(Nvar - ydelta, Nvar, "1", adj = c(1, 1))
+ text(Nvar - ydelta, 0.5 * Nvar, "-1", adj = c(1, 0))
+ text(Nvar - ydelta, 0.75 * Nvar, "0", adj = c(1, 0.5))
+ return()
+}
+
+"densplot" <-
+function (x, show.obs = TRUE, bwf, main = "", ylim, ...)
+{
+ xx <- as.matrix(x)
+ for (i in 1:nvar(x)) {
+ y <- xx[, i, drop = TRUE]
+ if (missing(bwf))
+ bwf <- function(x) {
+ x <- x[!is.na(as.vector(x))]
+ return(1.06 * min(sd(x), IQR(x)/1.34) * length(x)^-0.2)
+ }
+ bw <- bwf(y)
+ width <- 4 * bw
+ if (max(abs(y - floor(y))) == 0 || bw == 0)
+ hist(y, prob = TRUE, main = main, ...)
+ else {
+ scale <- "open"
+ if (max(y) <= 1 && 1 - max(y) < 2 * bw) {
+ if (min(y) >= 0 && min(y) < 2 * bw) {
+ scale <- "proportion"
+ y <- c(y, -y, 2 - y)
+ }
+ }
+ else if (min(y) >= 0 && min(y) < 2 * bw) {
+ scale <- "positive"
+ y <- c(y, -y)
+ }
+ else scale <- "open"
+ dens <- density(y, width = width)
+ if (scale == "proportion") {
+ dens$y <- 3 * dens$y[dens$x >= 0 & dens$x <=
+ 1]
+ dens$x <- dens$x[dens$x >= 0 & dens$x <= 1]
+ }
+ else if (scale == "positive") {
+ dens$y <- 2 * dens$y[dens$x >= 0]
+ dens$x <- dens$x[dens$x >= 0]
+ }
+ if(missing(ylim))
+ ylim <- c(0, max(dens$y))
+ plot(dens, ylab = "", main = main, type = "l",
+ xlab = paste("N =", niter(x), " Bandwidth =", formatC(dens$bw)),
+ ylim = ylim, ...)
+ if (show.obs)
+ lines(y[1:niter(x)], rep(max(dens$y)/100, niter(x)),
+ type = "h")
+ }
+ if (!is.null(varnames(x)) && is.null(list(...)$main))
+ title(paste("Density of", varnames(x)[i]))
+ }
+ return(invisible(x))
+}
+
+"read.jags" <- function (file = "jags.out", start, end, thin, quiet=FALSE)
+{
+ nc <- nchar(file)
+ if (nc > 3 && substring(file, nc - 3, nc) == ".out")
+ root <- substring(file, 1, nc - 4)
+ else root <- file
+ index.file = paste(root, ".ind", sep="")
+
+ read.coda(file, index.file, start, end, thin, quiet)
+}
+
+"read.openbugs" <- function(stem="", start, end, thin, quiet=FALSE)
+{
+ nchain <- 0
+ while (TRUE) {
+ output.file <- paste(stem,"CODAchain",nchain+1,".txt",sep="")
+ if (file.exists(output.file))
+ nchain <- nchain + 1
+ else
+ break
+ }
+
+ if (nchain==0)
+ stop("No output files found")
+
+ index.file <- paste(stem,"CODAindex.txt",sep="")
+ ans <- vector("list",nchain)
+ for (i in 1:nchain) {
+ output.file <- paste(stem,"CODAchain",i,".txt",sep="")
+ ans[[i]] <- read.coda(output.file, index.file, start, end, thin, quiet)
+ }
+ return(mcmc.list(ans))
+}
+
+"read.coda" <- function (output.file, index.file, start, end, thin,quiet=FALSE) {
+ index <- read.table(index.file,
+ row.names = 1, col.names = c("", "begin", "end"))
+ vnames <- row.names(index)
+ temp <- scan(output.file, what = list(iter = 0, val = 0), quiet = TRUE)
+ ## Do one pass through the data to see if we can construct
+ ## a regular time series easily
+ ##
+ start.vec <- end.vec <- thin.vec <- numeric(nrow(index))
+ for (i in 1:length(vnames)) {
+ iter.i <- temp$iter[index[i, "begin"]:index[i, "end"]]
+ thin.i <- unique(diff(iter.i))
+ thin.vec[i] <- if (length(thin.i) == 1)
+ thin.i
+ else NA
+ start.vec[i] <- iter.i[1]
+ end.vec[i] <- iter.i[length(iter.i)]
+ }
+ if (any(is.na(start.vec)) || any(thin.vec != thin.vec[1]) ||
+ any((start.vec - start.vec[1])%%thin.vec[1] != 0)) {
+ ##
+ ## Do it the brute force way
+ ##
+ iter <- sort(unique(temp$iter))
+ old.thin <- unique(diff(iter))
+ if (length(old.thin) == 1)
+ is.regular <- TRUE
+ else {
+ if (all(old.thin%%min(old.thin) == 0))
+ old.thin <- min(old.thin)
+ else old.thin <- 1
+ is.regular <- FALSE
+ }
+ }
+ else {
+ iter <- seq(from = min(start.vec), to = max(end.vec),
+ by = thin.vec[1])
+ old.thin <- thin.vec[1]
+ is.regular <- TRUE
+ }
+ if (missing(start))
+ start <- min(start.vec)
+ else if (start < min(start.vec)) {
+ warning("start not changed")
+ start <- min(start.vec)
+ }
+ else if (start > max(end.vec))
+ stop("Start after end of data")
+ else iter <- iter[iter >= start]
+ if (missing(end))
+ end <- max(end.vec)
+ else if (end > max(end.vec)) {
+ warning("end not changed")
+ end <- max(end.vec)
+ }
+ else if (end < min(start.vec))
+ stop("End before start of data")
+ else iter <- iter[iter <= end]
+ if (missing(thin))
+ thin <- old.thin
+ else if (thin%%old.thin != 0) {
+ thin <- old.thin
+ warning("thin not changed")
+ }
+ else {
+ new.iter <- iter[(iter - start)%%thin == 0]
+ new.thin <- unique(diff(new.iter))
+ if (length(new.thin) != 1 || new.thin != thin)
+ warning("thin not changed")
+ else {
+ iter <- new.iter
+ end <- max(iter)
+ is.regular <- TRUE
+ }
+ }
+ out <- matrix(NA, nrow = length(iter), ncol = nrow(index))
+ dimnames(out) <- list(iter, vnames)
+ for (v in vnames) {
+ if(!quiet)
+ cat("Abstracting", v, "... ")
+ inset <- index[v, "begin"]:index[v, "end"]
+ iter.v <- temp$iter[inset]
+ if (!is.regular) {
+ use.v <- duplicated(c(iter, iter.v))[-(1:length(iter))]
+ use <- duplicated(c(iter.v, iter))[-(1:length(iter.v))]
+ }
+ else {
+ use.v <- (iter.v - start)%%thin == 0 & iter.v >=
+ start & iter.v <= end
+ use <- (iter.v[use.v] - start)%/%thin + 1
+ }
+ if (any(use) & any(use.v))
+ out[use, v] <- temp$val[inset[use.v]]
+ if(!quiet)
+ cat(length(use), "valid values\n")
+ }
+ if (is.regular)
+ out <- mcmc(out, start = start, end = end, thin = thin)
+ else warning("not returning an mcmc object")
+ return(out)
+}
+
+"traceplot" <-
+function (x, smooth = TRUE, col = 1:6, type = "l", ylab = "",
+ ...)
+{
+ x <- mcmc.list(x)
+ args <- list(...)
+ for (j in 1:nvar(x)) {
+ xp <- as.vector(time(x))
+ yp <- if (nvar(x) > 1)
+ x[, j, drop = TRUE]
+ else x
+ yp <- do.call("cbind", yp)
+ matplot(xp, yp, xlab = "Iterations", ylab = ylab, type = type,
+ col = col, ...)
+ if (!is.null(varnames(x)) && is.null(list(...)$main))
+ title(paste("Trace of", varnames(x)[j]))
+ if (smooth) {
+ scol <- rep(col, length = nchain(x))
+ for (k in 1:nchain(x)) lines(lowess(xp, yp[, k]),
+ col = scol[k])
+ }
+ }
+}
+
+"plot.mcmc" <- function (x, trace = TRUE, density = TRUE, smooth = TRUE, bwf,
+ auto.layout = TRUE, ask = TRUE, ...)
+{
+ oldpar <- NULL
+ on.exit(par(oldpar))
+ if (auto.layout) {
+ mfrow <- set.mfrow(Nchains = nchain(x), Nparms = nvar(x),
+ nplots = trace + density)
+ oldpar <- par(mfrow = mfrow)
+ }
+ for (i in 1:nvar(x)) {
+ y <- mcmc(as.matrix(x)[, i, drop=FALSE], start(x), end(x), thin(x))
+ if (trace)
+ traceplot(y, smooth = smooth)
+ if (density) {
+ if (missing(bwf))
+ densplot(y)
+ else
+ densplot(y, bwf = bwf)
+ }
+ if (i==1)
+ oldpar <- c(oldpar, par(ask=ask))
+ }
+}
+
+"summary.mcmc" <-
+ function (object, quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975), ...)
+{
+ x <- as.mcmc(object)
+ statnames <- c("Mean", "SD", "Naive SE", "Time-series SE")
+ varstats <- matrix(nrow = nvar(x), ncol = length(statnames),
+ dimnames = list(varnames(x), statnames))
+ sp0 <- function(x) spectrum0(x)$spec
+ if (is.matrix(x)) {
+ xmean <- apply(x, 2, mean)
+ xvar <- apply(x, 2, var)
+ xtsvar <- apply(x, 2, sp0)
+ varquant <- t(apply(x, 2, quantile, quantiles))
+ }
+ else {
+ xmean <- mean(x, na.rm = TRUE)
+ xvar <- var(x, na.rm = TRUE)
+ xtsvar <- sp0(x)
+ varquant <- quantile(x, quantiles)
+ }
+ varstats[, 1] <- xmean
+ varstats[, 2] <- sqrt(xvar)
+ varstats[, 3] <- sqrt(xvar/niter(x))
+ varstats[, 4] <- sqrt(xtsvar/niter(x))
+ varstats <- drop(varstats)
+ varquant <- drop(varquant)
+ out <- list(statistics = varstats, quantiles = varquant,
+ start = start(x), end = end(x), thin = thin(x), nchain = 1)
+ class(out) <- "summary.mcmc"
+ return(out)
+}
+
+"print.summary.mcmc" <-
+ function (x, digits = max(3, .Options$digits - 3), ...)
+{
+ cat("\n", "Iterations = ", x$start, ":", x$end, "\n", sep = "")
+ cat("Thinning interval =", x$thin, "\n")
+ cat("Number of chains =", x$nchain, "\n")
+ cat("Sample size per chain =", (x$end - x$start)/x$thin +
+ 1, "\n")
+ cat("\n1. Empirical mean and standard deviation for each variable,")
+ cat("\n plus standard error of the mean:\n\n")
+ print(x$statistics, digits = digits, ...)
+ cat("\n2. Quantiles for each variable:\n\n")
+ print(x$quantiles, digits = digits, ...)
+ cat("\n")
+ invisible(x)
+}
+
+
+
+
+
+
+
+
diff --git a/R/raftery.R b/R/raftery.R
new file mode 100644
index 0000000..d0f9088
--- /dev/null
+++ b/R/raftery.R
@@ -0,0 +1,97 @@
+"raftery.diag" <-
+function (data, q = 0.025, r = 0.005, s = 0.95, converge.eps = 0.001)
+{
+ if (is.mcmc.list(data))
+ return(lapply(data, raftery.diag, q, r, s, converge.eps))
+ data <- as.mcmc(data)
+ resmatrix <- matrix(nrow = nvar(data), ncol = 4, dimnames = list(varnames(data,
+ allow.null = TRUE), c("M", "N", "Nmin", "I")))
+ phi <- qnorm(0.5 * (1 + s))
+ nmin <- as.integer(ceiling((q * (1 - q) * phi^2)/r^2))
+ if (nmin > niter(data))
+ resmatrix <- c("Error", nmin)
+ else for (i in 1:nvar(data)) {
+ # First need to find the thinning parameter kthin
+ #
+ if (is.matrix(data)) {
+ quant <- quantile(data[, i, drop = TRUE], probs = q)
+ dichot <- mcmc(data[, i, drop = TRUE] <= quant, start = start(data),
+ end = end(data), thin = thin(data))
+ }
+ else {
+ quant <- quantile(data, probs = q)
+ dichot <- mcmc(data <= quant, start = start(data),
+ end = end(data), thin = thin(data))
+ }
+ kthin <- 0
+ bic <- 1
+ while (bic >= 0) {
+ kthin <- kthin + thin(data)
+ testres <- as.vector(window.mcmc(dichot, thin = kthin))
+ newdim <- length(testres)
+ testtran <- table(testres[1:(newdim - 2)], testres[2:(newdim -
+ 1)], testres[3:newdim])
+ testtran <- array(as.double(testtran), dim = dim(testtran))
+ g2 <- 0
+ for (i1 in 1:2) {
+ for (i2 in 1:2) {
+ for (i3 in 1:2) {
+ if (testtran[i1, i2, i3] != 0) {
+ fitted <- (sum(testtran[i1, i2, 1:2]) *
+ sum(testtran[1:2, i2, i3]))/(sum(testtran[1:2,
+ i2, 1:2]))
+ g2 <- g2 + testtran[i1, i2, i3] * log(testtran[i1,
+ i2, i3]/fitted) * 2
+ }
+ }
+ }
+ }
+ bic <- g2 - log(newdim - 2) * 2
+ }
+ #
+ # then need to find length of burn-in and No of iterations for required precision
+ #
+ finaltran <- table(testres[1:(newdim - 1)], testres[2:newdim])
+ alpha <- finaltran[1, 2]/(finaltran[1, 1] + finaltran[1,
+ 2])
+ beta <- finaltran[2, 1]/(finaltran[2, 1] + finaltran[2,
+ 2])
+ tempburn <- log((converge.eps * (alpha + beta))/max(alpha,
+ beta))/(log(abs(1 - alpha - beta)))
+ nburn <- as.integer(ceiling(tempburn) * kthin)
+ tempprec <- ((2 - alpha - beta) * alpha * beta * phi^2)/(((alpha +
+ beta)^3) * r^2)
+ nkeep <- as.integer(ceiling(tempprec) * kthin)
+ iratio <- (nburn + nkeep)/nmin
+ resmatrix[i, 1] <- nburn
+ resmatrix[i, 2] <- nkeep + nburn
+ resmatrix[i, 3] <- nmin
+ resmatrix[i, 4] <- signif(iratio, digits = 3)
+ }
+ y <- list(params = c(r = r, s = s, q = q), resmatrix = resmatrix)
+ class(y) <- "raftery.diag"
+ return(y)
+}
+"print.raftery.diag" <-
+function (x, digits = 3, ...)
+{
+ cat("\nQuantile (q) =", x$params["q"])
+ cat("\nAccuracy (r) = +/-", x$params["r"])
+ cat("\nProbability (s) =", x$params["s"], "\n")
+ if (x$resmatrix[1] == "Error")
+ cat("\nYou need a sample size of at least", x$resmatrix[2],
+ "with these values of q, r and s\n")
+ else {
+ out <- x$resmatrix
+ for (i in ncol(out)) out[, i] <- format(out[, i], digits = digits)
+ out <- rbind(matrix(c("Burn-in ", "Total", "Lower bound ",
+ "Dependence", "(M)", "(N)", "(Nmin)", "factor (I)"),
+ byrow = TRUE, nrow = 2), out)
+ if (!is.null(rownames(x$resmatrix)))
+ out <- cbind(c("", "", rownames(x$resmatrix)), out)
+ dimnames(out) <- list(rep("", nrow(out)), rep("", ncol(out)))
+ print.default(out, quote = FALSE, ...)
+ cat("\n")
+ }
+ invisible(x)
+}
diff --git a/R/thin.R b/R/thin.R
new file mode 100644
index 0000000..e69de29
diff --git a/R/util.R b/R/util.R
new file mode 100644
index 0000000..32ee612
--- /dev/null
+++ b/R/util.R
@@ -0,0 +1,255 @@
+"read.yesno" <-
+function (string, default=TRUE)
+{
+ wrd <- ifelse(default, " (Y/n)?\n:", " (y/N)?\n:")
+ cat("\n", string, wrd, sep = "")
+ ans <- readline()
+ val <- if (default)
+ pmatch(ans, c("no","NO"), nomatch=0) == 0
+ else
+ pmatch(ans, c("yes","YES"), nomatch=0) != 0
+ return(val)
+}
+
+"change.tfoption" <-
+function (string, option)
+{
+ current.value <- coda.options(option)
+ if (!is.logical(current.value))
+ stop("Invalid option: must take logical values")
+ new.value <- read.yesno(string, current.value)
+ if (new.value != current.value) {
+ arg <- list(new.value)
+ names(arg) <- option
+ coda.options(arg)
+ }
+ return()
+}
+
+"coda.options" <-
+function (...)
+{
+ ## Set and display coda options
+ single <- FALSE
+ if (!exists(".Coda.Options", frame = 1))
+ .Coda.Options <<- .Coda.Options.Default
+ if (nargs() == 0) {
+ return(.Coda.Options)
+ }
+ else {
+ args <- list(...)
+ if (length(args) == 1) {
+ if (is.list(args[[1]]))
+ args <- args[[1]]
+ else if (is.null(names(args)))
+ single <- TRUE
+ }
+ }
+ if (is.null(names(args))) {
+ ## Display options
+ args <- unlist(args)
+ value <- vector("list", length(args))
+ names(value) <- args
+ for (v in args) if (any(v == names(.Coda.Options)))
+ value[v] <- .Coda.Options[v]
+ if (single)
+ return(value[[1]])
+ else return(value)
+ }
+ else {
+ ## Set options
+ oldvalue <- vector("list", length(args))
+ names(oldvalue) <- names(args)
+ if (any(names(args) == "default") && args$default ==
+ TRUE)
+ .Coda.Options <<- .Coda.Options.Default
+ for (v in names(args)) if (any(v == names(.Coda.Options))) {
+ oldvalue[v] <- .Coda.Options[v]
+ if (is.null(args[[v]]))
+ .Coda.Options[v] <<- list(NULL)
+ else if (mode(.Coda.Options[[v]]) == mode(args[[v]]))
+ .Coda.Options[v] <<- args[v]
+ }
+ invisible(oldvalue)
+ }
+}
+
+"multi.menu" <- function (choices, title, header, allow.zero = TRUE)
+{
+ ## Select more than one value from a menu
+ ##
+ if (!missing(title))
+ cat(title, "\n\n")
+ mat <- matrix(c(1:length(choices), choices), ncol = 2)
+ if (!missing(header)) {
+ if (length(header) == 2)
+ mat <- rbind(header, mat)
+ else stop("header is wrong length")
+ }
+ cat(paste(format(mat[, 1]), format(mat[, 2])), sep = "\n")
+ repeat {
+ cat("\nEnter relevant number(s), separated by commas",
+ "Ranges such as 3:7 may be specified)", sep = "\n")
+ if (allow.zero)
+ cat("(Enter 0 for none)\n")
+ ans <- scan(what = character(), sep = ",", strip.white = TRUE,
+ nlines = 1, quiet = TRUE)
+ if (length(ans) > 0) {
+ out <- numeric(0)
+ for (i in 1:length(ans)) {
+ nc <- nchar(ans[i])
+ wrd <- substring(ans[i], 1:nc, 1:nc)
+ colons <- wrd == ":"
+ err <- any(is.na(as.numeric(wrd[!colons]))) |
+ sum(colons) > 1 | colons[1] | colons[nc]
+ if (err) {
+ cat("Error: you have specified a non-numeric value!\n")
+ break
+ }
+ else {
+ out <- c(out, eval(parse(text = ans[i])))
+ if (min(out) < ifelse(allow.zero, 0, 1) | max(out) >
+ length(choices) | (any(out == 0) & length(out) >
+ 1)) {
+ err <- TRUE
+ cat("Error: you have specified variable number(s) out of range!\n")
+ break
+ }
+ }
+ }
+ if (!err)
+ break
+ }
+ }
+ return(out)
+}
+
+"read.and.check" <-
+ function (message = "", what = numeric(), lower, upper, answer.in, default)
+{
+ ## Read data from the command line and check that it satisfies
+ ## certain conditions. The function will loop until it gets
+ ## and answer satisfying the conditions. This entails extensive
+ ## checking of the conditions to make sure they are consistent
+ ## so we don't end up in an infinite loop.
+ have.lower <- !missing(lower)
+ have.upper <- !missing(upper)
+ have.ans.in <- !missing(answer.in)
+ have.default <- !missing(default)
+ if (have.lower | have.upper) {
+ if (!is.numeric(what))
+ stop("Can't have upper or lower limits with non numeric input")
+ if (have.lower && !is.numeric(lower))
+ stop("lower limit not numeric")
+ if (have.upper && !is.numeric(upper))
+ stop("upper limit not numeric")
+ if ((have.upper & have.lower) && upper < lower)
+ stop("lower limit greater than upper limit")
+ }
+ if (have.ans.in) {
+ if (mode(answer.in) != mode(what))
+ stop("inconsistent values of what and answer.in")
+ if (have.lower)
+ answer.in <- answer.in[answer.in >= lower]
+ if (have.upper)
+ answer.in <- answer.in[answer.in <= upper]
+ if (length(answer.in) == 0)
+ stop("No possible response matches conditions")
+ }
+ if (have.default) {
+ if (mode(default) != mode(what))
+ stop("inconsistent values of what and default")
+ if (have.lower && default < lower)
+ stop("default value below lower limit")
+ if (have.upper && default > upper)
+ stop("default value above upper limit")
+ if (have.ans.in && !any(answer.in == default))
+ stop("default value does not satisfy conditions")
+ }
+ err <- TRUE
+ while (err) {
+ if (nchar(message) > 0) {
+ cat("\n", message, "\n", sep = "")
+ if (have.default)
+ cat("(Default = ", default, ")\n", sep = "")
+ }
+ repeat {
+ cat("1:")
+ ans <- readline()
+ if (length(ans) == 1 && nchar(ans) > 0)
+ break
+ else if (have.default) {
+ ans <- default
+ break
+ }
+ }
+ if (is.numeric(what)) {
+ err1 <- TRUE
+ ans <- as.numeric(ans)
+ message <- "You must enter a number"
+ if (is.na(ans))
+ NULL
+ else if ((have.lower & have.upper) && (ans < lower |
+ ans > upper))
+ message <- paste(message, "between", lower, "and",
+ upper)
+ else if (have.lower && ans < lower)
+ message <- paste(message, ">=", lower)
+ else if (have.upper && ans > upper)
+ message <- paste(message, "<=", upper)
+ else err1 <- FALSE
+ }
+ else err1 <- FALSE
+ if (have.ans.in) {
+ if (!is.na(ans) && !any(ans == answer.in)) {
+ message <- paste("You must enter one of the following:",
+ paste(answer.in, collapse = ","))
+ err2 <- TRUE
+ }
+ else err2 <- FALSE
+ }
+ else err2 <- FALSE
+ err <- err1 | err2
+ }
+ return(ans)
+}
+
+".Coda.Options.Default" <-
+ list(trace = TRUE,
+ densplot = TRUE,
+ lowess = FALSE,
+ combine.plots = TRUE,
+ bandwidth = function (x)
+ {
+ x <- x[!is.na(x)]
+ 1.06 * min(sd(x), IQR(x)/1.34) * length(x)^-0.2
+ },
+ digits = 3,
+ quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975),
+ frac1 = 0.1,
+ frac2 = 0.5,
+ q = 0.025,
+ r = 0.005,
+ s = 0.95,
+ combine.stats = FALSE,
+ combine.corr = FALSE,
+ halfwidth = 0.1,
+ user.layout = FALSE,
+ gr.bin = 10,
+ geweke.nbin = 20,
+ gr.max = 50,
+ data.saved = TRUE
+ )
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/README b/README
new file mode 100644
index 0000000..60afc5c
--- /dev/null
+++ b/README
@@ -0,0 +1,37 @@
+CODA is a set of tools for analyzing the output of Markov Chain
+Monte Carlo (MCMC) simulations and diagnosing lack of convergence.
+
+The S original can be downloaded from
+http://www.mrc-bsu.cam.ac.uk
+and is Copyright (C) MRC Biostatistics Unit 1995.
+
+The CODA S manual contained these acknowledgements
+"The support of the Economic and Social Research Council (UK) is
+gratefully acknowledged. The work was funded in part by ESRC (UK)
+award number H519 25 5023. We are also grateful to Brad Carlin
+for many helpful comments and ideas concerning the CODA software
+and manual, and to Steve Brooks for suggesting the graphical
+implementations of the Geweke and Gelman & Rubin convergence
+diagnostics."
+
+See "CHANGELOG" for information on the changes in the R version.
+
+Martyn Plummer <plummer at iarc.fr> 20/5/1998
+
+######################################################################
+Copying CODA for R
+
+`CODA' 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, or (at your option) any later
+version.
+
+`CODA' 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.
+
+A copy of the GNU General Public License is available on the World Wide
+Web at http://www.gnu.org/copyleft/gpl.html. You can also obtain it by
+writing to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge,
+MA 02139, USA.
diff --git a/data/line.R b/data/line.R
new file mode 100644
index 0000000..d9c397a
--- /dev/null
+++ b/data/line.R
@@ -0,0 +1,230 @@
+"line" <-
+structure(list(line1 = structure(c(7.17313, 2.95253, 3.66989,
+3.31522, 3.70544, 3.5791, 2.70206, 2.96136, 3.53406, 2.09471,
+3.06524, 2.73901, 3.08929, 2.70413, 2.52174, 3.73216, 2.67339,
+3.01779, 2.70123, 2.97837, 2.87646, 3.19134, 2.79708, 3.25934,
+2.73704, 2.71863, 3.01749, 2.984, 2.73648, 2.55982, 2.94331,
+3.27086, 3.55436, 2.80668, 3.18001, 3.07977, 2.77286, 2.74566,
+3.67006, 2.24932, 3.08757, 2.56762, 2.80828, 3.76523, 3.10116,
+3.01902, 2.65918, 3.52618, 2.43358, 2.57103, 3.55572, 3.08102,
+3.06265, 2.57951, 2.69831, 2.6604, 3.1768, 2.92196, 1.73505,
+3.46348, 2.99482, 3.18707, 3.2433, 3.3528, 3.12026, 3.27727,
+3.30783, 2.68048, 3.31009, 2.8779, 3.21578, 3.27217, 3.21431,
+3.4163, 2.72993, 3.83811, 2.7769, 3.10277, 0.858013, 1.72668,
+2.2129, 3.13023, 3.11833, 2.52852, 3.19425, 2.33602, 3.86727,
+2.45181, 2.53356, 2.33611, 2.71644, 3.36774, 2.49752, 2.98218,
+2.49989, 2.92943, 2.86987, 4.08664, 2.92904, 2.87937, 3.5247,
+2.6994, 3.2374, 2.79161, 3.0065, 2.84989, 2.92812, 3.09409, 2.98947,
+3.05823, 3.16941, 3.00681, 2.87709, 2.76276, 3.18145, 2.36106,
+2.81812, 3.08823, 2.65743, 2.8572, 3.33833, 2.80197, 2.95402,
+2.95727, 3.19123, 3.43451, 2.59671, 3.60602, 2.68363, 3.56141,
+3.16793, 3.2417, 3.3542, 2.98939, 2.34584, 3.03147, 3.16633,
+2.56506, 3.01896, 3.17286, 2.80791, 2.92504, 3.02911, 2.85911,
+3.28071, 3.43015, 2.3898, 1.67055, 3.81164, 2.43832, 2.92065,
+1.37547, 3.23811, 2.96203, 2.83998, 2.36596, 3.51832, 2.59301,
+3.08633, 2.87675, 3.16977, 2.81822, 2.48206, 3.28322, 2.51132,
+3.01153, 3.25603, 2.98599, 2.93238, 2.12555, 3.5398, 2.96786,
+2.73922, 2.76325, 2.89286, 3.41665, 2.84412, 2.69208, 3.87634,
+2.70917, 3.23347, 3.25577, 2.70044, 2.83369, 3.2501, 3.12599,
+3.16206, 2.90532, 3.35174, 3.36722, 2.73076, 3.26922, 2.81658,
+2.95479, 2.71722, 3.02135, 3.13198, 4.07261, 2.84968, 2.69468,
+-1.5662, 1.50337, 0.628157, 1.18272, 0.490437, 0.20697, 0.882553,
+1.08515, 1.06926, 1.48077, 0.378349, 0.444712, -0.01306, 1.15585,
+0.856368, 0.657461, 0.594411, 0.88777, 0.745918, 0.993619, 1.02806,
+0.553611, 1.01886, 0.470345, 0.520167, 0.89238, 0.68696, 0.798954,
+0.688647, 0.579934, 0.844434, 0.810385, 0.872753, 0.895408, 0.827359,
+0.919973, 0.926044, 0.522444, 1.05851, 1.32707, 0.674701, 0.746217,
+1.14114, 1.24147, 0.791581, 0.938256, 0.522189, 0.718395, 0.928954,
+0.57796, 1.07794, 0.812967, 1.04835, 0.993271, 0.781727, 0.306433,
+0.58444, 0.55981, 0.88081, 1.01651, 0.876193, 0.650802, 0.75265,
+0.948834, 0.557148, 0.99243, 0.6986, 0.78952, 0.570903, 0.725568,
+1.16533, 0.453463, 0.871724, 0.889113, 0.507765, 0.414795, 0.87789,
+0.633498, 0.736897, -0.568545, 0.770004, 0.869198, 0.665951,
+0.758507, 0.742198, 1.02311, 0.969699, 0.580553, 1.14088, 0.582853,
+0.677865, 0.652573, 0.952581, 1.1587, 0.999131, 1.09971, 1.02759,
+0.571576, 0.756901, 0.927734, 0.753859, 0.494313, 0.510782, 0.793101,
+1.10834, 0.658604, 0.788326, 0.866193, 0.884201, 0.507964, 0.179421,
+0.718163, 0.874446, 1.01166, 0.626119, 0.808673, 0.52271, 1.06545,
+0.902426, 0.380372, 0.659257, 0.68071, 0.794731, 0.680552, 0.508912,
+0.494933, 1.17763, 1.22197, 0.545816, 0.965099, 0.677669, 0.688658,
+1.06926, 0.344008, 1.09068, 1.02111, 0.691322, 0.364369, 0.915177,
+0.566777, 1.15452, 0.625684, 0.708469, 0.984351, 0.576846, 0.228625,
+1.40307, 1.21012, 0.143615, 1.4046, 1.32799, 0.0254974, 1.17005,
+1.28567, 1.51562, 0.842245, 0.852613, 0.780749, 0.967637, 0.810536,
+0.858167, 0.960625, 0.613509, 1.38389, 1.22199, 0.93934, 0.479022,
+0.684411, 0.907639, 0.694256, 0.762828, 1.08536, 0.737203, 0.79847,
+0.605098, 1.01106, 1.29919, -0.116987, 0.645465, 0.868414, 0.934646,
+0.841893, 0.864243, 1.03498, 0.784894, 1.28125, 0.594917, 1.23255,
+0.788097, 0.929231, 0.823524, 0.748782, 0.625626, 0.742992, 0.733012,
+0.598408, 0.654669, 1.19037, 0.472082, 0.669647, 11.2331, 4.88649,
+1.39734, 0.662879, 1.36213, 1.0435, 1.29043, 0.459322, 0.634257,
+0.912919, 1.1885, 0.576963, 2.40033, 1.79423, 0.943078, 0.903465,
+0.731041, 0.622143, 0.790993, 0.740969, 0.58258, 0.710707, 0.870071,
+1.00433, 1.03022, 1.08587, 0.572809, 0.716056, 0.636699, 0.563391,
+0.848744, 0.516244, 0.642656, 0.554613, 0.695308, 0.617855, 0.425071,
+0.603549, 0.764704, 0.907059, 0.840521, 0.812011, 0.612573, 1.70232,
+0.644559, 0.481628, 0.641244, 0.805918, 0.957325, 0.676466, 0.743322,
+0.951599, 0.595749, 0.462307, 0.746122, 0.958645, 0.764432, 0.527422,
+1.74938, 1.09919, 0.493587, 0.473889, 0.442078, 0.590437, 0.913574,
+0.667202, 0.772119, 0.511031, 0.47078, 0.551933, 0.785661, 0.902284,
+0.420702, 0.872218, 1.51717, 1.27235, 0.70261, 0.504571, 2.0531,
+2.2977, 1.38751, 0.818279, 0.326157, 0.79268, 0.783918, 0.714941,
+0.756979, 0.972415, 0.917166, 0.822013, 1.10835, 1.14671, 0.949399,
+1.01348, 1.07893, 0.564832, 0.704571, 1.1337, 1.58693, 0.463972,
+0.734966, 0.897794, 1.30018, 0.909768, 0.640371, 0.541355, 0.685204,
+1.09255, 0.531485, 0.693595, 0.814454, 1.08119, 0.497489, 0.497872,
+0.822327, 0.4964, 0.566723, 0.606257, 0.711031, 0.683499, 1.26004,
+1.13398, 0.529169, 0.538057, 0.564507, 0.493755, 1.00047, 1.24981,
+0.973904, 0.859757, 0.462606, 0.588109, 0.691391, 0.666279, 0.927818,
+1.05845, 0.45543, 0.79399, 0.885761, 0.513252, 0.637226, 0.698318,
+0.404093, 0.979977, 1.42097, 0.761441, 1.0802, 2.09035, 2.85473,
+2.44305, 1.81988, 1.78082, 1.69385, 0.643862, 1.5311, 0.960032,
+0.69531, 0.821253, 0.49187, 0.394229, 0.806075, 0.557793, 0.694918,
+1.46079, 1.56847, 0.823876, 0.73726, 0.698766, 0.800476, 1.45114,
+1.15917, 0.692247, 0.610882, 0.709356, 0.630025, 0.81496, 2.35916,
+1.30581, 1.1819, 1.11921, 0.618407, 0.821751, 0.76678, 0.845531,
+0.519626, 0.636336, 0.82726, 0.896255, 0.786207, 0.959944, 0.528469,
+0.607334, 0.434402, 0.527194, 1.04028, 0.40121, 0.693556, 1.07389,
+1.40933, 0.702048), .Dim = c(200, 3), title = "Line example from BUGS manual", .Dimnames = list(
+ c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11",
+ "12", "13", "14", "15", "16", "17", "18", "19", "20", "21",
+ "22", "23", "24", "25", "26", "27", "28", "29", "30", "31",
+ "32", "33", "34", "35", "36", "37", "38", "39", "40", "41",
+ "42", "43", "44", "45", "46", "47", "48", "49", "50", "51",
+ "52", "53", "54", "55", "56", "57", "58", "59", "60", "61",
+ "62", "63", "64", "65", "66", "67", "68", "69", "70", "71",
+ "72", "73", "74", "75", "76", "77", "78", "79", "80", "81",
+ "82", "83", "84", "85", "86", "87", "88", "89", "90", "91",
+ "92", "93", "94", "95", "96", "97", "98", "99", "100", "101",
+ "102", "103", "104", "105", "106", "107", "108", "109", "110",
+ "111", "112", "113", "114", "115", "116", "117", "118", "119",
+ "120", "121", "122", "123", "124", "125", "126", "127", "128",
+ "129", "130", "131", "132", "133", "134", "135", "136", "137",
+ "138", "139", "140", "141", "142", "143", "144", "145", "146",
+ "147", "148", "149", "150", "151", "152", "153", "154", "155",
+ "156", "157", "158", "159", "160", "161", "162", "163", "164",
+ "165", "166", "167", "168", "169", "170", "171", "172", "173",
+ "174", "175", "176", "177", "178", "179", "180", "181", "182",
+ "183", "184", "185", "186", "187", "188", "189", "190", "191",
+ "192", "193", "194", "195", "196", "197", "198", "199", "200"
+ ), c("alpha", "beta", "sigma")), mcpar = c(1, 200, 1), class = "mcmc"),
+ line2 = structure(c(2.0665, 2.55073, 3.72809, 2.90303, 3.06816,
+ 2.98528, 2.3895, 3.34781, 2.72255, 3.03422, 2.63764, 3.06537,
+ 2.81315, 2.81342, 3.29965, 3.34629, 2.49348, 3.10066, 3.36426,
+ 3.12697, 3.34863, 3.10297, 3.74189, 2.83072, 2.9638, 2.79184,
+ 3.32553, 2.68574, 3.27625, 3.34029, 3.02742, 2.65711, 3.84285,
+ 1.96594, 3.97947, 2.29779, 2.84843, 3.13026, 2.90137, 2.99816,
+ 3.11093, 2.90318, 3.45998, 3.13388, 3.07965, 3.18346, 2.24565,
+ 3.16838, 3.20878, 2.97041, 2.03604, 1.92021, 3.88952, 2.55067,
+ 3.21949, 3.14698, 3.05184, 2.87204, 3.32694, 3.26103, 3.57177,
+ 3.91897, 4.0938, 2.88649, 2.22315, 3.23495, 3.14733, 3.32536,
+ 2.39064, 2.84845, 2.92939, 2.7519, 2.95218, 3.2797, 3.18452,
+ 2.25467, 3.0187, 3.139, 3.52716, 2.874, 3.15899, 3.32123,
+ 2.88329, 2.81646, 2.65933, 3.35541, 2.77517, 2.72268, 3.23746,
+ 2.95345, 3.31671, 3.03783, 3.17906, 2.68964, 2.90325, 3.1247,
+ 3.29658, 2.17373, 3.20434, 2.73006, 3.07409, 3.16792, 3.15965,
+ 2.92753, 2.12887, 2.3563, 3.34344, 3.06189, 3.18122, 2.88713,
+ 3.28045, 3.10015, 3.00738, 3.02902, 2.77719, 3.27406, 2.93017,
+ 1.78484, 2.30908, 3.26235, 1.69871, 2.17906, 4.32621, 1.92995,
+ 4.06669, 3.20001, 2.57108, 3.02103, 3.15325, 3.27674, 2.0916,
+ 2.09049, 3.31506, 3.65649, 2.70318, 2.81422, 3.24301, 2.97542,
+ 3.30149, 2.28222, 3.0605, 2.92056, 3.04161, 3.08682, 2.67589,
+ 3.05074, 2.65665, 2.76893, 3.19402, 2.64889, 3.52223, 3.4263,
+ 2.66415, 3.14702, 3.22297, 2.72178, 3.26504, 3.1093, 3.09429,
+ 3.0888, 2.73249, 2.4253, 3.11828, 3.42255, 3.3503, 3.19407,
+ 3.29954, 3.2595, 2.88135, 2.64103, 2.73082, 2.79848, 3.84443,
+ 3.88629, 3.05079, 2.82624, 3.18078, 2.32302, 2.64219, 2.68628,
+ 3.10314, 0.995679, 3.86631, 3.17386, 3.4696, 2.47629, 2.91303,
+ 3.31464, 3.17376, 3.48446, 2.39874, 2.9073, 3.51113, 3.25472,
+ 3.2314, 3.40695, 3.07729, 3.1929, 3.16783, 3.04118, 0.94983,
+ 0.533868, 0.520681, 0.976253, 0.936261, 0.663807, 0.520043,
+ 0.944512, 1.05836, 0.829194, 0.427276, 1.00063, 1.0042, 0.676282,
+ 0.724315, 0.378694, 0.854646, 0.816283, 0.680898, 0.729862,
+ 0.466017, 0.474517, 0.394259, 0.643551, 1.41783, 0.83389,
+ 0.770501, 0.60024, 0.931097, 1.03179, 0.520711, 0.576376,
+ 0.757683, -0.0119585, 0.150777, 1.24154, 0.586833, 0.687874,
+ 0.972612, 0.616836, 0.951993, 0.654136, 0.835886, 0.760564,
+ 1.31897, 0.689678, 0.906503, 1.17308, 0.784572, 1.08129,
+ 1.45739, 1.78349, 0.508111, 0.935584, 1.13076, 0.642275,
+ 0.774581, 0.598707, 0.62828, 1.04412, 0.41416, 1.19133, 1.51109,
+ 1.76453, 1.83201, 1.22612, 1.06988, 0.671452, 0.554991, 0.975127,
+ 0.790377, 0.902487, 0.731783, 1.0601, 1.34457, 0.285294,
+ 0.391361, 0.46474, 0.975072, 0.698322, 0.526451, 0.643863,
+ 0.671277, 0.793108, 0.691103, 1.02774, 0.396104, 0.521036,
+ 0.802256, 0.695039, 1.14746, 0.583661, 0.502895, 1.073, 0.48887,
+ 0.622952, 0.652116, 0.925127, 0.563, 0.94581, 0.863503, 0.818531,
+ 0.467898, 1.20696, 0.989445, 0.251698, 1.30748, 1.20417,
+ 0.971032, 1.02162, 0.853551, 0.52759, 0.655414, 0.944856,
+ 0.851747, 0.795153, 1.00725, 0.393811, 0.492267, 0.513625,
+ 0.429794, 0.46022, 0.0291329, 1.29284, 0.965173, 0.986593,
+ 0.884349, 0.450296, 0.979947, 1.03558, 1.19383, 0.00606791,
+ 1.5062, 0.536381, 0.815264, 0.408098, 0.729014, 0.544239,
+ 1.06494, 0.918886, 1.08346, 0.807998, 0.721469, 0.791212,
+ 0.797347, 1.12917, 0.948997, 0.610172, 1.00315, 0.936662,
+ 0.99261, 0.714404, 0.87753, 0.801131, 1.04282, 0.856488,
+ 0.974811, 0.980089, 0.808844, 0.648608, 1.46944, 0.480841,
+ 0.637054, 0.858352, 0.58538, 0.498672, 0.549712, 0.950073,
+ 1.0674, 0.788781, 1.17857, 1.38812, 0.822348, 1.22831, 1.12554,
+ 0.978673, 0.777699, 0.725505, 0.490426, 1.8001, -0.423911,
+ 0.401579, 0.462218, 1.06891, 0.799904, 0.242961, 0.935728,
+ 0.785555, 0.598434, 0.794666, 1.04866, 0.121866, 1.4973,
+ 1.05516, 0.890826, 1.08431, 0.894609, 0.650802, 0.73541,
+ 0.711719, 2.85379, 1.71367, 0.785154, 0.948167, 0.900588,
+ 0.353863, 1.08933, 0.970812, 0.687744, 0.638359, 0.685602,
+ 1.08756, 0.945919, 1.0119, 0.396287, 0.703951, 1.09031, 0.520708,
+ 0.508404, 0.653224, 0.697648, 1.38661, 1.70559, 0.911467,
+ 0.791402, 2.10078, 0.615926, 0.712817, 0.507729, 0.607383,
+ 0.683574, 0.607109, 0.861146, 0.956142, 3.11952, 1.72381,
+ 0.775019, 0.549051, 0.378787, 0.606155, 1.15792, 0.96476,
+ 0.860189, 0.543162, 0.581212, 0.597179, 0.883207, 0.612349,
+ 0.453471, 0.560563, 1.83515, 2.71558, 1.86329, 0.95504, 0.819032,
+ 0.856319, 0.513872, 0.615228, 0.458746, 0.812665, 0.880437,
+ 1.57162, 1.35811, 1.77924, 3.2142, 1.47968, 0.747181, 0.824633,
+ 0.722867, 0.985975, 0.44006, 0.857335, 0.8669, 0.883645,
+ 0.685724, 0.607935, 1.16049, 0.705714, 1.20252, 1.07224,
+ 0.661277, 0.577825, 0.684492, 0.472472, 0.590746, 0.919054,
+ 1.0548, 0.931485, 1.10547, 0.587061, 0.815852, 1.00739, 0.98054,
+ 0.76707, 0.677419, 0.432805, 0.686665, 0.836332, 0.829481,
+ 0.943772, 0.642839, 0.589167, 0.785612, 0.631841, 1.6026,
+ 1.15419, 1.48332, 0.819761, 0.99212, 0.679848, 0.706879,
+ 0.611507, 0.703271, 0.575749, 0.483907, 0.459734, 0.774127,
+ 2.67997, 1.33679, 1.14551, 1.28857, 1.58951, 2.05222, 3.26272,
+ 1.34444, 0.701195, 0.93947, 0.657195, 0.759601, 0.554493,
+ 1.08315, 2.5578, 1.17271, 1.28146, 0.575432, 0.915557, 1.02713,
+ 0.680568, 0.732328, 0.759689, 1.41759, 0.656574, 0.626601,
+ 0.542805, 0.762587, 0.581579, 1.42557, 0.575034, 0.782904,
+ 0.737252, 0.764629, 0.528021, 1.02315, 0.402441, 0.505607,
+ 0.885029, 0.428758, 0.503349, 0.531532, 0.502562, 1.31589,
+ 2.06585, 0.699819, 0.747237, 0.646591, 0.681749, 0.893148,
+ 0.529854, 0.353713, 0.919825, 1.03411, 1.47011, 1.20727,
+ 0.836974, 0.897857, 1.0101, 0.66153, 0.772771, 0.708564,
+ 2.06246, 2.36283, 2.6266, 2.18514, 1.4795, 0.959802, 1.67397,
+ 1.14405, 0.46466, 0.943977, 0.853832, 1.00929, 1.14416, 2.24906,
+ 1.0981, 0.606293, 1.30693, 0.846828, 0.465129, 0.672417,
+ 0.639787), .Dim = c(200, 3), title = "Line example from BUGS manual", .Dimnames = list(
+ c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10",
+ "11", "12", "13", "14", "15", "16", "17", "18", "19",
+ "20", "21", "22", "23", "24", "25", "26", "27", "28",
+ "29", "30", "31", "32", "33", "34", "35", "36", "37",
+ "38", "39", "40", "41", "42", "43", "44", "45", "46",
+ "47", "48", "49", "50", "51", "52", "53", "54", "55",
+ "56", "57", "58", "59", "60", "61", "62", "63", "64",
+ "65", "66", "67", "68", "69", "70", "71", "72", "73",
+ "74", "75", "76", "77", "78", "79", "80", "81", "82",
+ "83", "84", "85", "86", "87", "88", "89", "90", "91",
+ "92", "93", "94", "95", "96", "97", "98", "99", "100",
+ "101", "102", "103", "104", "105", "106", "107", "108",
+ "109", "110", "111", "112", "113", "114", "115", "116",
+ "117", "118", "119", "120", "121", "122", "123", "124",
+ "125", "126", "127", "128", "129", "130", "131", "132",
+ "133", "134", "135", "136", "137", "138", "139", "140",
+ "141", "142", "143", "144", "145", "146", "147", "148",
+ "149", "150", "151", "152", "153", "154", "155", "156",
+ "157", "158", "159", "160", "161", "162", "163", "164",
+ "165", "166", "167", "168", "169", "170", "171", "172",
+ "173", "174", "175", "176", "177", "178", "179", "180",
+ "181", "182", "183", "184", "185", "186", "187", "188",
+ "189", "190", "191", "192", "193", "194", "195", "196",
+ "197", "198", "199", "200"), c("alpha", "beta", "sigma"
+ )), mcpar = c(1, 200, 1), class = "mcmc")), .Names = c("line1",
+"line2"), class = "mcmc.list")
diff --git a/man/Cramer.Rd b/man/Cramer.Rd
new file mode 100644
index 0000000..ad23ea5
--- /dev/null
+++ b/man/Cramer.Rd
@@ -0,0 +1,28 @@
+\name{Cramer}
+\alias{pcramer}
+\title{The Cramer-von Mises Distribution}
+\description{
+ Distribution function of the Cramer-von Mises distribution.
+}
+\usage{
+ pcramer(q, eps)
+}
+\arguments{
+ \item{q}{vector of quantiles.}
+ \item{eps}{accuracy required}
+}
+\value{
+ \code{pcramer} gives the distribution function,
+}
+\references{
+Anderson TW. and Darling DA. Asymptotic theory of certain `goodness
+of fit' criteria based on stochastic processes. \emph{Ann. Math.
+Statist.}, \bold{23}, 192-212 (1952).
+
+Csorgo S. and Faraway, JJ. The exact and asymptotic distributions of
+the Cramer-von Mises statistic. J. Roy. Stat. Soc. (B), \bold{58},
+221-234 (1996).
+}
+
+\keyword{distribution}
+
diff --git a/man/as.ts.mcmc.Rd b/man/as.ts.mcmc.Rd
new file mode 100644
index 0000000..d5c60b2
--- /dev/null
+++ b/man/as.ts.mcmc.Rd
@@ -0,0 +1,24 @@
+\name{as.ts.mcmc}
+\alias{as.ts.mcmc}
+\title{Coerce mcmc object to time series}
+
+\usage{as.ts.mcmc(x)}
+
+\arguments{
+ \item{x}{an mcmc object}
+}
+
+\description{
+\code{as.ts.mcmc} will coerce an mcmc object to a time series.}
+
+\note{
+This is pretending to be an \code{as.ts} method for mcmc objects. But
+\code{as.ts} is not generic.}
+
+\author{Martyn Plummer}
+
+\seealso{
+\code{\link{as.ts}}
+}
+\keyword{ts}
+\keyword{utilities}
diff --git a/man/autocorr.Rd b/man/autocorr.Rd
new file mode 100644
index 0000000..3f111ac
--- /dev/null
+++ b/man/autocorr.Rd
@@ -0,0 +1,30 @@
+\name{autocorr}
+\alias{autocorr}
+\title{Autocorrelation function for Markov chains}
+
+\usage{autocorr(mcmc.obj, lags = c(0, 1, 5, 10, 50), relative=TRUE }
+
+\description{
+\code{autocorr} calculates the autocorrelation function for the
+Markov chain \code{mcmc.obj} at the lags given by \code{lags}.
+The lag values are taken to be relative to the thinning interval
+if \code{relative=TRUE}.
+
+High autocorrelations within chains indicate slow mixing and, usually,
+slow convergence. It may be useful to thin out a chain with high
+autocorrelations before calculating summary statistics: a thinned
+chain may contain most of the information, but take up less space in
+memory. Re-running the MCMC sampler with a different parameterization
+may help to reduce autocorrelation.
+}
+
+\value{
+A vector or array containing the autocorrelations.
+}
+
+\author{Martyn Plummer}
+
+\seealso{
+\code{\link{acf}}, \code{\link{autocorr.plot}}.}
+}
+\keyword{ts}
diff --git a/man/autocorr.plot.Rd b/man/autocorr.plot.Rd
new file mode 100644
index 0000000..a4cd5b0
--- /dev/null
+++ b/man/autocorr.plot.Rd
@@ -0,0 +1,24 @@
+\name{autocorr.plot}
+\alias{autocorr.plot}
+\title{Plot autocorrelations for Markov Chains}
+
+\usage{autocorr.plot(x, lag.max, auto.layout = TRUE, ask=TRUE, ...)}
+
+\arguments{
+\item{x}{A Markov Chain}
+\item{lag.max}{Maximum value at which to calculate acf}
+\item{auto.layout}{If \code{TRUE} then, set up own layout for
+plots, otherwise use existing one.}
+\item{ask}{If \code{TRUE} then prompt user before displaying
+ each page of plots.}
+\item{...}{graphical parameters}
+}
+
+\description{
+Plots the autocorrelation function for each variable in each chain in x.
+}
+
+\seealso{
+ \code{\link{autocorr}}.
+}
+\keyword{hplot}
diff --git a/man/bugs2jags.Rd b/man/bugs2jags.Rd
new file mode 100644
index 0000000..4dffc22
--- /dev/null
+++ b/man/bugs2jags.Rd
@@ -0,0 +1,33 @@
+\name{bugs2jags}
+\alias{bugs2jags}
+\title{Convert WinBUGS data file to JAGS data file}
+\usage{
+bugs2jags(infile, outfile)
+}
+\arguments{
+ \item{infile}{name of the input file}
+ \item{outfile}{name of the output file}
+}
+\description{
+ \code{bugs2jags} converts a WinBUGS data in the format called "S-Plus"
+ (i.e. the format created by the \code{dput} function) and writes it in
+ \code{dump} format used by JAGS.
+
+ NB WinBUGS stores its arrays in row order. This is different
+ from R and JAGS which both store arrays in column order. This
+ difference is taken into account by \code{bugs2jags} which will
+ automatically reorder the data in arrays, without changing the
+ dimension.
+}
+\references{
+ Spiegelhalter DJ, Thomas A, Best NG and Lunn D (2003).
+ \emph{WinBUGS version 1.4 user manual}
+ MRC Biostatistics Unit, Cambridge, UK.
+}
+\author{Martyn Plummer}
+\seealso{
+ \code{\link{dput}},
+ \code{\link{dump}}.
+}
+\keyword{file}
+
diff --git a/man/coda.options.Rd b/man/coda.options.Rd
new file mode 100644
index 0000000..49342a9
--- /dev/null
+++ b/man/coda.options.Rd
@@ -0,0 +1,77 @@
+\name{coda.options}
+\alias{coda.options}
+\alias{display.coda.options}
+\alias{.Coda.Options}
+\alias{.Coda.Options.Default}
+\title{Options settings for the codamenu driver}
+\usage{
+coda.options(...)
+display.coda.options(data = FALSE, stats = FALSE, plots = FALSE, diags = FALSE)
+.Coda.Options
+.Coda.Options.Default
+}
+\arguments{
+ \item{data}{logical flag: show data options?}
+ \item{stats}{logical flag: show summary statistic options?}
+ \item{plots}{logical flag: show plotting options?}
+ \item{diags}{logical flag: show plotting options?}
+ \item{...}{list of options}
+}
+\description{
+ \code{coda.options} is a utility function that queries and sets
+ options for the \code{codamenu()} function. These settings affect
+ the behaviour of the functions in the coda library only when they
+ are called via the \code{codamenu()} interface.
+
+ The \code{coda.options()} function behaves just like the
+ \code{options()} function in the base library, with the additional
+ feature that \code{coda.options(default=TRUE)} will reset all options
+ to the default values.
+
+ Options can be pretty-printed using the \code{display.coda.options()}
+ function, which groups the options into sections.
+
+ Available options are
+ \describe{
+ \item{bandwidth}{Bandwidth function used when smoothing samples to
+ produce density estimates. Defaults to Silverman's ``Rule of thumb''.}
+ \item{combine.corr}{Logical option that determines whether to
+ combine multiple chains when calculating cross-correlations.}
+ \item{combine.plots}{Logical option that determines whether to
+ combine multiple chains when plotting.}
+ \item{combine.plots}{Logical option that determines whether to
+ combine multiple chains when calculating summary statistics.}
+ \item{data.saved}{For internal use only.}
+ \item{densplot}{Logical option that determines whether to plot
+ a density plot when plot methods are called for mcmc objects.}
+ \item{digits}{Number of significant digits to use when printing.}
+ \item{frac1}{For Geweke diagnostic, fraction to use from
+ start of chain. Defaults to 0.1}
+ \item{frac2}{For Geweke diagnostic, fraction to use from
+ end of chain. Default to 0.5.}
+ \item{gr.bin}{For Geweke-Brooks plot, number of iterations to use
+ per bin.}
+ \item{gr.max}{For Geweke-Brooks plot, maximum number of bins to
+ use. This option overrides \code{gr.bin}.}
+ \item{halfwidth}{For Heidelberger and Welch diagnostic, the target
+ value for the ratio of half width to sample mean.}
+ \item{lowess}{Logical option that controls whether to plot a
+ smooth line through a trace plot when plotting MCMC output.}
+ \item{q}{For Raftery and Lewis diagnostic, the target quantile to
+ be estimated}
+ \item{r}{For Raftery and Lewis diagnostic, the required precision.}
+ \item{s}{For Raftery and Lewis diagnostic, the probability of
+ obtaining an estimate in the interval (q-r, q+r).}
+ \item{quantiles}{Vector of quantiles to print when calculating
+ summary statistics for MCMC output.}
+ \item{trace}{Logical option that determines whether to plot a trace
+ of the sampled output when plotting MCMC output.}
+ \item{user.layout}{Logical option that determines whether current
+ value of par("mfrow") should be used for plots (TRUE) or whether
+ the optimal layout should be calculated (FALSE).}
+ }
+}
+\seealso{
+ \code{\link{options}}
+}
+\keyword{utilities}
diff --git a/man/codamenu.Rd b/man/codamenu.Rd
new file mode 100644
index 0000000..c16f8f7
--- /dev/null
+++ b/man/codamenu.Rd
@@ -0,0 +1,15 @@
+\name{codamenu}
+\alias{codamenu}
+\title{Main menu driver for the coda package}
+
+\usage{
+codamenu()
+}
+
+\description{
+\code{codamenu} presents a simple menu-based interface to the functions
+in the coda package. It is designed for users who know nothing about
+the R language.}
+
+\author{Kate Cowles, Nicky Best, Karen Vines, Martyn Plummer}
+\keyword{utilities}
diff --git a/man/crosscorr.Rd b/man/crosscorr.Rd
new file mode 100644
index 0000000..fdaaa86
--- /dev/null
+++ b/man/crosscorr.Rd
@@ -0,0 +1,25 @@
+\name{crosscorr}
+\alias{crosscorr}
+\title{Cross correlations for MCMC output}
+
+\usage{crosscorr(x)}
+
+\arguments{
+ \item{x}{an \code{mcmc} or \code{mcmc.list} object.}
+}
+\description{
+ \code{crosscorr} calculates cross-correlations between
+ variables in Markov Chain Monte Carlo output. If \code{x}
+ is an mcmc.list then all chains in \code{x} are combined
+ before calculating the correlation.
+}
+
+\value{
+A matrix or 3-d array containing the correlations.
+}
+
+\seealso{
+ \code{\link{crosscorr.plot}}, \code{\link{autocorr}}.
+}
+\keyword{multivariate}
+\keyword{array}
diff --git a/man/crosscorr.plot.Rd b/man/crosscorr.plot.Rd
new file mode 100644
index 0000000..18df54a
--- /dev/null
+++ b/man/crosscorr.plot.Rd
@@ -0,0 +1,27 @@
+\name{crosscorr.plot}
+\alias{crosscorr.plot}
+\title{Plot image of correlation matrix}
+
+\usage{crosscorr.plot (x, col = topo.colors(10), ...)}
+
+\arguments{
+\item{x}{an \code{mcmc} or \code{mcmc.list} object}
+\item{col}{color palette to use}
+\item{...}{graphical parameters}
+}
+
+\description{
+\code{crosscorr.plot} provides an image of the correlation matrix for
+\code{x}. If \code{x} is an \code{mcmc.list} object, then all chains
+are combined.
+
+The range [-1,1] is divided into a number of equal-length categories
+given by the length of \code{col} and assigned the corresponding color.
+By default, topographic colours are used as this makes it easier to
+distinguish positive and negative correlations.
+}
+
+\seealso{
+\code{\link{crosscorr}}, \code{\link{image}}, \code{\link{topo.colors}}.
+}
+\keyword{hplot}
diff --git a/man/cumuplot.Rd b/man/cumuplot.Rd
new file mode 100644
index 0000000..0e60a4a
--- /dev/null
+++ b/man/cumuplot.Rd
@@ -0,0 +1,30 @@
+\name{cumuplot}
+\alias{cumuplot}
+\title{Cumulative quantile plot}
+
+\usage{
+ cumuplot(x, probs=c(0.025,0.5,0.975), ylab="",
+ lty=c(2,1), lwd=c(1,2), type="l", ask=TRUE, auto.layout=TRUE,
+ col=1, ...)
+}
+
+\arguments{
+ \item{x}{an \code{mcmc} object}
+ \item{probs}{vector of desired quantiles}
+ \item{ylab, lty, lwd, type, col}{graphical parameters}
+ \item{auto.layout}{If \code{TRUE}, then set up own layout for
+ plots, otherwise use existing one.}
+ \item{ask}{If \code{TRUE} then prompt user before displaying
+ each page of plots.}
+ \item{...}{further graphical parameters}
+}
+
+\description{
+ Plots the evolution of the sample quantiles as a function of the
+ number of iterations.
+}
+
+\author{
+ Arni Magnusson <arnima at u.washington.edu>
+}
+\keyword{hplot}
diff --git a/man/densplot.Rd b/man/densplot.Rd
new file mode 100644
index 0000000..9baf342
--- /dev/null
+++ b/man/densplot.Rd
@@ -0,0 +1,37 @@
+\name{densplot}
+\alias{densplot}
+\title{Probability density function estimate from MCMC output}
+
+\usage{densplot(x, show.obs=TRUE, bwf, main, ylim, ...)}
+
+\arguments{
+ \item{x}{An \code{mcmc} or \code{mcmc.list} object}
+ \item{show.obs}{Show observations along the x-axis}
+ \item{bwf}{Function for calculating the bandwidth. If omitted,
+ the bandwidth is calculate by 1.06 times the minimum of the standard
+ deviation and the interquartile range divided by 1.34 times the sample
+ size to the negative one fifth power}
+ \item{main}{Title. See \code{par()}}
+ \item{ylim}{Limits on y axis. See \code{par()}}
+ \item{...}{Further graphical parameters}
+}
+
+\description{
+ Displays a plot of the density estimate for each variable in \code{x},
+ calculated by the \code{density} function.
+}
+
+\note{
+ You can call this function directly, but it is more usually
+ called by the \code{plot.mcmc} function.
+
+ If a variable is bounded below at 0, or bounded in the interval [0,1],
+ then the data are reflected at the boundary before being passed to the
+ density() function. This allows correct estimation of a non-zero
+ density at the boundary.
+}
+
+\seealso{
+ \code{\link{density}}, \code{\link{plot.mcmc}}.
+}
+\keyword{hplot}
diff --git a/man/effectiveSize.Rd b/man/effectiveSize.Rd
new file mode 100644
index 0000000..56359c2
--- /dev/null
+++ b/man/effectiveSize.Rd
@@ -0,0 +1,28 @@
+\name{effectiveSize}
+\alias{effectiveSize}
+\title{Effective sample size for estimating the mean}
+\description{
+Sample size adjusted for autocorrelation.
+}
+\usage{
+effectiveSize(x)
+}
+\arguments{
+\item{x}{An mcmc or mcmc.list object.}
+}
+\details{
+For a time series \code{x} of length \code{N}, the standard error of
+the mean is \code{var(x)/n} where \code{n} is the effective sample
+size. \code{n = N} only when there is no autocorrelation.
+
+Estimation of the effective sample size requires estimating the
+spectral density at frequency zero. This is done by the function
+\code{spectrum0.ar}
+}
+\value{
+A vector giving the effective sample size for each column of \code{x}.
+}
+\seealso{
+ \code{\link{spectrum0.ar}}.
+}
+\keyword{ts}
diff --git a/man/gelman.diag.Rd b/man/gelman.diag.Rd
new file mode 100644
index 0000000..e58492f
--- /dev/null
+++ b/man/gelman.diag.Rd
@@ -0,0 +1,96 @@
+\name{gelman.diag}
+\alias{gelman.diag}
+%\alias{gelman.mv.diag}
+%\alias{gelman.transform}
+%\alias{print.gelman.diag}
+\title{Gelman and Rubin's convergence diagnostic}
+
+\usage{gelman.diag(x, confidence = 0.95, transform=FALSE)}
+
+\arguments{
+\item{x}{An \code{mcmc.list} object with more than one chain,
+and with starting values that are overdispersed with respect to the
+posterior distribution.}
+\item{confidence}{the coverage probability of the confidence interval for the
+potential scale reduction factor}
+\item{transform}{a logical flag indicating whether variables in
+\code{x} should be transformed to improve the normality of the
+distribution. If set to TRUE, a log transform or logit transform, as
+appropriate, will be applied.}
+}
+
+\description{
+The `potential scale reduction factor' is calculated for each variable in
+\code{x}, together with upper and lower confidence limits. Approximate
+convergence is diagnosed when the upper limit is close to 1. For
+multivariate chains, a multivariate value is calculated that bounds
+above the potential scale reduction factor for any linear combination
+of the (possibly transformed) variables.
+
+The confidence limits are based on the assumption that the stationary
+distribution of the variable under examination is normal. Hence the
+`transform' parameter may be used to improve the normal approximation.
+}
+
+\section{Theory}{
+Gelman and Rubin (1992) propose a general approach to monitoring
+convergence of MCMC output in which two or more parallel chains are run
+with starting values that are overdispersed relative to the posterior
+distribution. Convergence is diagnosed when the chains have `forgotten'
+their initial values, and the output from all chains is indistinguishable.
+The \code{gelman.diag} diagnostic is applied to a single variable from
+the chain. It is based a comparison of within-chain and between-chain
+variances, and is similar to a classical analysis of variance.
+
+There are two ways to estimate the variance of the stationary distribution:
+the mean of the empirical variance within each chain, \eqn{W}, and
+the empirical variance from all chains combined, which can be expressed as
+\deqn{ \widehat{\sigma}^2 =
+ \frac{(n-1) B }{n} + \frac{W}{n} }{ sigma.hat^2 = (n-1)B/n + W/n }
+where \eqn{B} is the empirical between-chain variance.
+
+If the chains have converged, then both estimates are unbiased. Otherwise
+the first method will \emph{underestimate} the variance, since the
+individual chains have not had time to range all over the stationary
+distribution, and the second method will \emph{overestimate} the variance,
+since the starting points were chosen to be overdispersed.
+
+The convergence diagnostic is based on the assumption that the
+target distribution is normal. A Bayesian credible interval can
+be constructed using a t-distribution with mean
+\deqn{\widehat{\mu}=\mbox{Sample mean of all chains
+combined}}{mu.hat = Sample mean of all chains combined}
+and variance
+\deqn{\widehat{V}=\widehat{\sigma}^2 + \frac{B}{mn}}{V.hat=sigma.hat2 + B/(mn)}
+where \eqn{m} is the number of chains, and degrees of freedom
+estimated by the method of moments
+\deqn{d = \frac{2*\widehat{V}}{\mbox{Var}(\widehat{V})}}{d = 2*V.hat/Var(V.hat)}
+Use of the t-distribution accounts for the fact that the mean
+and variance of the posterior distribution are estimated.
+
+The convergence diagnostic itself is
+\deqn{R=\sqrt{\frac{(d+3) \widehat{V}}{(d+1)W}}}{R=sqrt((d+3) V.hat /((d+1)W)}
+Values substantially above 1 indicate lack of convergence. If the
+chains have not converged, Bayesian credible intervals based on the
+t-distribution are too wide, and have the potential to shrink by this
+factor if the MCMC run is continued.
+}
+
+\note{
+The multivariate a version of Gelman and Rubin's diagnostic was
+proposed by Brooks and Gelman (1997).
+}
+
+\references{
+Gelman, A and Rubin, DB (1992) Inference from iterative simulation
+using multiple sequences, \emph{Statistical Science}, \bold{7}, 457-511.
+
+Brooks, SP. and Gelman, A. (1997) General methods for monitoring
+convergence of iterative simulations. \emph{Journal of Computational and
+Graphical Statistics}, \bold{7}, 434-455.
+}
+
+\seealso{
+ \code{\link{gelman.plot}}.
+}
+\keyword{htest}
diff --git a/man/gelman.plot.Rd b/man/gelman.plot.Rd
new file mode 100644
index 0000000..2c5d40d
--- /dev/null
+++ b/man/gelman.plot.Rd
@@ -0,0 +1,60 @@
+\name{gelman.plot}
+\alias{gelman.plot}
+%\alias{gelman.preplot}
+\title{Gelman-Rubin-Brooks plot}
+
+\usage{
+gelman.plot(x, bin.width = 10, max.bins = 50,
+confidence = 0.95, transform = FALSE, auto.layout = TRUE, ask = TRUE,
+col, lty, xlab, ylab, type, ...)
+}
+
+\arguments{
+ \item{x}{an mcmc object}
+ \item{bin.width}{Number of observations per segment, excluding the
+ first segment which always has at least 50 iterations.}
+ \item{max.bins}{Maximum number of bins, excluding the last one.}
+ \item{confidence}{Coverage probability of confidence interval.}
+ \item{transform}{Automatic variable transformation (see \code{gelman.diag})}
+ \item{auto.layout}{If \code{TRUE} then, set up own layout for
+ plots, otherwise use existing one.}
+ \item{ask}{Prompt user before displaying each page of plots.}
+ \item{col}{graphical parameter (see \code{par})}
+ \item{lty}{graphical parameter (see \code{par})}
+ \item{xlab}{graphical parameter (see \code{par})}
+ \item{ylab}{graphical parameter (see \code{par})}
+ \item{type}{graphical parameter (see \code{par})}
+ \item{...}{further graphical parameters.}
+}
+
+\description{
+This plot shows the evolution of Gelman and Rubin's shrink factor as
+the number of iterations increases.
+}
+
+\details{
+The Markov chain is divided into bins according to the arguments
+\code{bin.width} and \code{max.bins}. Then the Gelman-Rubin shrink factor
+is repeatedly calculated. The first shrink factor is calculated with
+observations 1:50, the second with observations \eqn{1:(50+n)} where n is
+the bin width, the third contains samples \eqn{1:(50+2n)} and so on.
+}
+
+\references{
+Brooks, S P. and Gelman, A. (1998) General Methods for Monitoring
+Convergence of Iterative Simulations. Journal of Computational and
+Graphical Statistics. 7. p434-455.
+}
+
+\section{Theory}{
+A potential problem with \code{gelman.diag} is that it may mis-diagnose
+convergence if the shrink factor happens to be close to 1 by chance.
+By calculating the shrink factor at several points in time,
+\code{gelman.plot} shows if the shrink factor has really converged, or
+whether it is still fluctuating.
+}
+
+\seealso{
+ \code{\link{gelman.diag}}.
+}
+\keyword{hplot}
diff --git a/man/geweke.diag.Rd b/man/geweke.diag.Rd
new file mode 100644
index 0000000..b365c8f
--- /dev/null
+++ b/man/geweke.diag.Rd
@@ -0,0 +1,48 @@
+\name{geweke.diag}
+\alias{geweke.diag}
+%\alias{print.geweke.diag}
+\title{Geweke's convergence diagnostic}
+
+\usage{geweke.diag(x, frac1=0.1, frac2=0.5)}
+
+\arguments{
+ \item{x}{an mcmc object}
+ \item{frac1}{fraction to use from beginning of chain}
+ \item{frac2}{fraction to use from end of chain}
+}
+
+\value{
+Z-scores for a test of equality of means
+between the first and last parts of the chain. A separate
+statistic is calculated for each variable in each chain.
+}
+
+\description{
+Geweke (1992) proposed a convergence diagnostic for Markov chains
+based on a test for equality of the means of the first and last part
+of a Markov chain (by default the first 10\% and the last 50\%). If the
+samples are drawn from the stationary distribution of the chain, the two
+means are equal and Geweke's statistic has an asymptotically standard
+normal distribution.
+
+The test statistic is a standard Z-score: the difference between the
+two sample means divided by its estimated standard error. The standard
+error is estimated from the spectral density at zero and so takes into
+account any autocorrelation.
+
+The Z-score is calculated under the assumption that the two parts of
+the chain are asymptotically independent, which requires that the sum
+of \code{frac1} and \code{frac2} be strictly less than 1.
+}
+
+\seealso{
+ \code{\link{geweke.plot}}.
+}
+
+\references{
+Geweke, J. Evaluating the accuracy of sampling-based approaches
+to calculating posterior moments. In \emph{Bayesian Statistics 4}
+(ed JM Bernado, JO Berger, AP Dawid and AFM Smith). Clarendon Press,
+Oxford, UK.}
+
+\keyword{htest}
diff --git a/man/geweke.plot.Rd b/man/geweke.plot.Rd
new file mode 100644
index 0000000..bd54ee2
--- /dev/null
+++ b/man/geweke.plot.Rd
@@ -0,0 +1,47 @@
+\name{geweke.plot}
+\alias{geweke.plot}
+\title{Geweke-Brooks plot}
+
+\usage{
+geweke.plot(x, frac1 = 0.1, frac2 = 0.5, nbins = 20,
+ pvalue = 0.05, auto.layout = TRUE, ask = TRUE, ...)
+}
+
+\arguments{
+ \item{x}{an mcmc object}
+ \item{frac1}{fraction to use from beginning of chain.}
+ \item{frac2}{fraction to use from end of chain.}
+ \item{nbins}{Number of segments.}
+ \item{pvalue}{p-value used to plot confidence limits for the null hypothesis.}
+ \item{auto.layout}{If \code{TRUE} then, set up own layout for
+ plots, otherwise use existing one.}
+ \item{ask}{Prompt user before displaying each page of plots.}
+ \item{...}{Graphical parameters.}
+}
+
+\description{
+If \code{geweke.diag} indicates that the first and last part of a sample
+from a Markov chain are not drawn from the same distribution, it may
+be useful to discard the first few iterations to see if the rest of the
+chain has "converged". This plot shows what happens to Geweke's Z-score
+when successively larger numbers of iterations are discarded from the
+beginning of the chain. To preserve the asymptotic conditions required
+for Geweke's diagnostic, the plot never discards more than half the chain.
+
+The first half of the Markov chain is divided into \code{nbins - 1}
+segments, then Geweke's Z-score is repeatedly calculated. The first
+Z-score is calculated with all iterations in the chain, the second
+after discarding the first segment, the third after discarding the first
+two segments, and so on. The last Z-score is calculated using only the
+samples in the second half of the chain.
+}
+
+\note{
+The graphical implementation of Geweke's diagnostic was suggested
+by Steve Brooks.
+}
+
+\seealso{
+ \code{\link{geweke.diag}}.
+}
+\keyword{hplot}
diff --git a/man/heidel.diag.Rd b/man/heidel.diag.Rd
new file mode 100644
index 0000000..aa3222e
--- /dev/null
+++ b/man/heidel.diag.Rd
@@ -0,0 +1,70 @@
+\name{heidel.diag}
+\alias{heidel.diag}
+%\alias{print.heidel.diag}
+\title{Heidelberger and Welch's convergence diagnostic}
+
+\usage{heidel.diag(x, eps=0.1, pvalue=0.05)}
+
+\arguments{
+\item{x} {An \code{mcmc} object}
+\item{eps}{Target value for ratio of halfwidth to sample mean}
+\item{pvalue}{significance level to use}
+}
+
+\description{
+\code{heidel.diag} is a run length control diagnostic based on a criterion
+of relative accuracy for the estimate of the mean. The default setting
+corresponds to a relative accuracy of two significant digits.
+
+\code{heidel.diag} also implements a convergence diagnostic, and removes
+up to half the chain in order to ensure that the means are estimated
+from a chain that has converged.
+}
+
+\details{
+The convergence test uses the Cramer-von-Mises statistic to test the null
+hypothesis that the sampled values come from a stationary distribution.
+The test is successively applied, firstly to the whole chain, then after
+discarding the first 10\%, 20\%, ... of the chain until either the null
+hypothesis is accepted, or 50\% of the chain has been discarded. The latter
+outcome constitutes `failure' of the stationarity test and indicates
+that a longer MCMC run is needed. If the stationarity test is passed,
+the number of iterations to keep and the number to discard are reported.
+
+The half-width test calculates a 95\% confidence interval for the
+mean, using the portion of the chain which passed the stationarity test.
+Half the width of this interval is compared with the estimate of the mean.
+If the ratio between the half-width and the mean is lower than \code{eps},
+the halfwidth test is passed. Otherwise the length of the sample is
+deemed not long enough to estimate the mean with sufficient accuracy.
+}
+
+\section{Theory}{
+The \code{heidel.diag} diagnostic is based on the work of Heidelberger
+and Welch (1983), who combined their earlier work on simulation run
+length control (Heidelberger and Welch, 1981) with the work of Schruben
+(1982) on detecting initial transients using Brownian bridge theory.
+}
+
+\note{
+If the half-width test fails then the run should be extended. In order
+to avoid problems caused by sequential testing, the test should not
+be repeated too frequently. Heidelberger and Welch (1981) suggest
+increasing the run length by a factor I > 1.5, each time, so that
+estimate has the same, reasonably large, proportion of new data.
+}
+
+\references{
+Heidelberger P and Welch PD. A spectral method for confidence interval
+generation and run length control in simulations. Comm. ACM. \bold{24},
+233-245 (1981)
+
+Heidelberger P and Welch PD. Simulation run length control in the
+presence of an initial transient. \emph{Opns Res.}, \bold{31},
+1109-44 (1983)
+
+Schruben LW. Detecting initialization bias in simulation experiments.
+\emph{Opns. Res.}, \bold{30}, 569-590 (1982).
+}
+
+\keyword{htest}
diff --git a/man/linepost.Rd b/man/linepost.Rd
new file mode 100644
index 0000000..28944c4
--- /dev/null
+++ b/man/linepost.Rd
@@ -0,0 +1,16 @@
+\name{line}
+\docType{data}
+\alias{line}
+\title{Simple linear regression example}
+\description{
+ Sample MCMC output from a simple linear regression model given in
+ the BUGS manual.
+}
+\usage{data(line)}
+\format{An \code{mcmc} object}
+\source{
+ Spiegelhalter, D.J., Thomas, A., Best, N.G. and Gilks, W.R. (1995)
+ BUGS: Bayesian inference using Gibbs Sampling, Version 0.5,
+ MRC Biostatistics Unit, Cambridge.
+}
+\keyword{datasets}
diff --git a/man/mcmc.Rd b/man/mcmc.Rd
new file mode 100644
index 0000000..9e7a982
--- /dev/null
+++ b/man/mcmc.Rd
@@ -0,0 +1,56 @@
+\name{mcmc}
+\alias{mcmc}
+\alias{as.mcmc}
+\alias{as.mcmc.default}
+\alias{is.mcmc}
+\alias{print.mcmc}
+\title{Markov Chain Monte Carlo Objects}
+
+\usage{
+mcmc(data= NA, start = 1, end = numeric(0), thin = 1)
+as.mcmc(x)
+is.mcmc(x)
+}
+
+\arguments{
+ \item{data}{a vector or matrix of MCMC output}
+ \item{start}{the iteration number of the first observation}
+ \item{end}{the iteration number of the last observation}
+ \item{thin}{the thinning interval between consecutive observations}
+ \item{x}{An object that may be coerced to an mcmc object}
+}
+
+\description{
+
+ The function `mcmc' is used to create a Markov Chain Monte Carlo object.
+ The data are taken to be a vector, or a matrix with one column per
+ variable.
+
+ An mcmc object may be summarized by the \code{summary} function
+ and visualized with the \code{plot} function.
+
+ MCMC objects resemble time series (\code{ts}) objects and have
+ methods for the generic functions \code{time}, \code{start},
+ \code{end}, \code{frequency} and \code{window}.
+}
+
+\author{Martyn Plummer}
+
+\note{
+ The format of the mcmc class has changed between coda version 0.3
+ and 0.4. Older mcmc objects will now cause \code{is.mcmc} to
+ fail with an appropriate warning message. Obsolete mcmc objects can
+ be upgraded with the \code{mcmcUpgrade} function.
+}
+
+
+\seealso{
+ \code{\link{mcmc.list}},
+ \code{\link{mcmcUpgrade}},
+ \code{\link{thin}},
+ \code{\link{window.mcmc}},
+ \code{\link{summary.mcmc}},
+ \code{\link{plot.mcmc}}.
+}
+
+\keyword{ts}
diff --git a/man/mcmc.convert.Rd b/man/mcmc.convert.Rd
new file mode 100644
index 0000000..b671318
--- /dev/null
+++ b/man/mcmc.convert.Rd
@@ -0,0 +1,46 @@
+\name{mcmc.convert}
+\alias{as.matrix.mcmc}
+\alias{as.matrix.mcmc.list}
+\alias{as.array.mcmc.list}
+\alias{as.mcmc.mcmc.list}
+\title{Conversions of MCMC objects}
+
+\usage{
+as.matrix(x, iters = FALSE, chains = FALSE))
+as.array.mcmc.list(x, drop, ...)
+}
+
+\arguments{
+ \item{x}{An \code{mcmc} or \code{mcmc.list} object}
+ \item{iters}{logical flag: add column for iteration number?}
+ \item{chains}{logical flag: add column for chain number? (if mcmc.list)}
+ \item{drop}{logical flag: if \code{TRUE} the result is coerced to the
+ lowest possible dimension}
+ \item{...}{further arguments for future methods}
+}
+
+\description{
+ These are methods for the generic functions \code{as.matrix()},
+ \code{as.array()} and \code{as.mcmc()}.
+
+ \code{as.matrix()} strips the MCMC attributes from an \code{mcmc}
+ object and returns a matrix. If \code{iters = TRUE} then a
+ column is added with the iteration number. For \code{mcmc.list}
+ objects, the rows of multiple chains are concatenated and, if
+ \code{chains = TRUE} a column is added with the chain number.
+
+ \code{mcmc.list} objects can be coerced to 3-dimensional arrays
+ with the \code{as.array()} function.
+
+ An \code{mcmc.list} object with a single chain can be coerced
+ to an \code{mcmc} object with \code{as.mcmc()}. If the argument
+ has multiple chains, this causes an error.
+}
+
+\seealso{
+ \code{\link{as.matrix}},
+ \code{\link{as.array}},
+ \code{\link{as.mcmc}},
+}
+
+\keyword{array}
diff --git a/man/mcmc.list.Rd b/man/mcmc.list.Rd
new file mode 100644
index 0000000..865b792
--- /dev/null
+++ b/man/mcmc.list.Rd
@@ -0,0 +1,56 @@
+\name{mcmc.list}
+\alias{mcmc.list}
+\alias{as.mcmc.list}
+\alias{as.mcmc.list.default}
+\alias{is.mcmc.list}
+\alias{plot.mcmc.list}
+\title{Replicated Markov Chain Monte Carlo Objects}
+
+\usage{
+mcmc.list(...)
+as.mcmc.list(x, ...)
+is.mcmc.list(x)
+}
+
+\arguments{
+ \item{...}{a list of mcmc objects}
+ \item{x}{an object that may be coerced to mcmc.list}
+}
+
+
+\description{
+
+The function `mcmc.list' is used to represent parallel runs of the same
+chain, with different starting values and random seeds. The list must
+be balanced: each chain in the list must have the same iterations and
+the same variables.
+
+Diagnostic functions which act on \code{mcmc} objects may also be applied
+to \code{mcmc.list} objects. In general, the chains will be combined,
+if this makes sense, otherwise the diagnostic function will be applied
+separately to each chain in the list.
+
+Since all the chains in the list have the same iterations, a single time
+dimension can be ascribed to the list. Hence there are time series methods
+\code{time}, \code{window}, \code{start}, \code{end}, \code{frequency}
+and \code{thin} for \code{mcmc.list} objects.
+
+An \code{mcmc.list} can be indexed as if it were a single mcmc object
+using the \code{[} operator (see examples below). The \code{[[} operator
+selects a single \code{mcmc} object from the list.
+}
+
+\author{Martyn Plummer}
+
+\seealso{
+ \code{\link{mcmc}}.
+}
+
+\examples{
+data(line)
+x1 <- line[[1]] #Select first chain
+x2 <- line[,1, drop=FALSE] #Select first var from all chains
+varnames(x2) == varnames(line)[1] #TRUE
+}
+
+\keyword{ts}
diff --git a/man/mcmc.subset.Rd b/man/mcmc.subset.Rd
new file mode 100644
index 0000000..100c55b
--- /dev/null
+++ b/man/mcmc.subset.Rd
@@ -0,0 +1,33 @@
+\name{mcmc.subset}
+\alias{[.mcmc}
+\alias{[.mcmc.list}
+\title{Extract or replace parts of MCMC objects}
+
+\usage{
+x[i,j]
+}
+
+\arguments{
+ \item{x}{An \code{mcmc} object}
+ \item{i}{Row to extract}
+ \item{j}{Column to extract}
+}
+
+\description{
+
+ These are methods for subsetting \code{mcmc} objects. You can select
+ iterations using the first dimension and variables using the second
+ dimension. Selecting iterations will return a vector or matrix, not an
+ \code{mcmc} object. If you want to do row-subsetting of an \code{mcmc}
+ object and preserve its dimensions, use the \code{window} function.
+
+ Subsetting applied to an \code{mcmc.list} object will simultaneously
+ affect all the parallel chains in the object.
+}
+
+\seealso{
+ \code{\link{[}},
+ \code{\link{window.mcmc}}
+}
+
+\keyword{ts}
diff --git a/man/mcmcUpgrade.Rd b/man/mcmcUpgrade.Rd
new file mode 100644
index 0000000..ce55fcf
--- /dev/null
+++ b/man/mcmcUpgrade.Rd
@@ -0,0 +1,29 @@
+\name{mcmcUpgrade}
+\alias{mcmcUpgrade}
+\title{Upgrade mcmc objects in obsolete format}
+
+\usage{
+ mcmcUpgrade(x)
+}
+
+\arguments{
+ \item{x}{an obsolete \code{mcmc} object.}
+}
+
+\description{
+ In previous releases of CODA, an \code{mcmc} object could
+ be a single or multiple chains. A new class \code{mcmc.list}
+ has now been introduced to deal with multiple chains and
+ \code{mcmc} objects can only have data from a single chain.
+
+ Objects stored in the old format are now obsolete
+ and must be upgraded.
+}
+
+\author{Martyn Plummer}
+
+\seealso{
+ \code{\link{mcmc}}.
+}
+
+\keyword{ts}
diff --git a/man/mcpar.Rd b/man/mcpar.Rd
new file mode 100644
index 0000000..2261087
--- /dev/null
+++ b/man/mcpar.Rd
@@ -0,0 +1,26 @@
+\name{mcpar}
+\alias{mcpar}
+\title{Mcpar attribute of MCMC objects}
+
+\usage{
+mcpar(x)
+}
+
+\arguments{
+\item{x}{An \code{mcmcm} or \code{mcmc.list} object}
+}
+
+\description{
+ The `mcpar' attribute of an MCMC object gives the start iteration
+ the end iteration and the thinning interval of the chain.
+
+ It resembles the `tsp' attribute of time series (\code{ts}) objects.
+}
+
+\seealso{
+ \code{\link{ts}},
+ \code{\link{mcmc}},
+ \code{\link{mcmc.list}},
+}
+
+\keyword{ts}
diff --git a/man/multi.menu.Rd b/man/multi.menu.Rd
new file mode 100644
index 0000000..39f48ce
--- /dev/null
+++ b/man/multi.menu.Rd
@@ -0,0 +1,37 @@
+\name{multi.menu}
+\alias{multi.menu}
+\title{Choose multiple options from a menu}
+
+\usage{
+multi.menu(choices, title, header, allow.zero = TRUE)
+}
+
+\arguments{
+ \item{choices}{Character vector of labels for choices}
+ \item{title}{Title printed before menu}
+ \item{header}{Character vector of length 2 giving column titles}
+ \item{allow.zero}{Permit 0 as an acceptable response}
+}
+
+\description{
+\code{multi.menu} presents the user with a menu of choices
+labelled from 1 to the number of choices. The user may choose
+one or more options by entering a comma separated list. A range
+of values may also be specified using the ":" operator. Mixed
+expressions such as "1,3:5, 6" are permitted.
+
+If \code{allow.zero} is set to TRUE, one can select `0' to exit
+without choosing an item.
+}
+
+\value{
+Numeric vector giving the numbers of the options selected, or
+0 if no selection is made.
+}
+
+\author{Martyn Plummer}
+
+\seealso{
+ \code{\link{menu}}.
+}
+\keyword{utilities}
diff --git a/man/nchain.Rd b/man/nchain.Rd
new file mode 100644
index 0000000..1eb4a09
--- /dev/null
+++ b/man/nchain.Rd
@@ -0,0 +1,36 @@
+\name{nchain}
+\alias{niter}
+\alias{nvar}
+\alias{nchain}
+\title{Dimensions of MCMC objects}
+
+\usage{
+niter(x)
+nvar(x)
+nchain(x)
+}
+
+\arguments{
+ \item{x}{An \code{mcmc} or \code{mcmc.list} object}
+}
+
+\value{
+ A numeric vector of length 1:
+}
+
+\description{
+ These functions give the dimensions of an MCMC object
+
+ \describe{
+ \item{niter(x)}{returns the number of iterations.}
+ \item{nvar(x)}{returns the number of variables.}
+ \item{chain(x)}{returns the number of parallel chains.}
+ }
+}
+
+\seealso{
+ \code{\link{mcmc}},
+ \code{\link{mcmc.list}},
+}
+
+\keyword{ts}
diff --git a/man/plot.mcmc.Rd b/man/plot.mcmc.Rd
new file mode 100644
index 0000000..4a92491
--- /dev/null
+++ b/man/plot.mcmc.Rd
@@ -0,0 +1,32 @@
+\name{plot.mcmc}
+\alias{plot.mcmc}
+\title{Summary plots of mcmc objects}
+
+\usage{
+\method{plot}{mcmc}(x, trace = TRUE, density = TRUE, smooth = TRUE, bwf,
+ auto.layout = TRUE, ask = TRUE, ...)
+}
+
+\arguments{
+ \item{x}{an object of class \code{mcmc} or \code{mcmc.list}}
+ \item{trace}{Plot trace of each variable}
+ \item{density}{Plot density estimate of each variable}
+ \item{smooth}{Draw a smooth line through trace plots}
+ \item{bwf}{Bandwidth function for density plots}
+ \item{auto.layout}{Automatically generate output format}
+ \item{ask}{Prompt user before each page of plots}
+ \item{...}{Further arguments}
+}
+
+\description{
+ \code{plot.mcmc} summarizes an mcmc or mcmc.list object
+ with a trace of the sampled output and a density estimate
+ for each variable in the chain.
+}
+
+\author{Martyn Plummer}
+
+\seealso{
+ \code{\link{densplot}}, \code{\link{traceplot}}.
+}
+\keyword{hplot}
diff --git a/man/raftery.diag.Rd b/man/raftery.diag.Rd
new file mode 100644
index 0000000..5e53b8a
--- /dev/null
+++ b/man/raftery.diag.Rd
@@ -0,0 +1,88 @@
+\name{raftery.diag}
+\alias{raftery.diag}
+%\alias{print.raftery.diag}
+\title{Raftery and Lewis's diagnostic}
+
+\usage{
+raftery.diag(data, q=0.025, r=0.005, s=0.95, converge.eps=0.001)
+}
+
+\arguments{
+ \item{data}{an \code{mcmc} object}
+ \item{q}{the quantile to be estimated.}
+ \item{r}{the desired margin of error of the estimate.}
+ \item{s}{the probability of obtaining an estimate in the interval (q-r,q+r).}
+ \item{converge.eps}{Precision required for estimate of time to convergence.}
+}
+
+\description{
+ \code{raftery.diag} is a run length control diagnostic based on a
+ criterion of accuracy of estimation of the quantile \code{q}. It is
+ intended for use on a short pilot run of a Markov chain. The number
+ of iterations required to estimate the quantile \eqn{q} to within an
+ accuracy of +/- \eqn{r} with probability \eqn{p} is calculated. Separate
+ calculations are performed for each variable within each chain.
+
+ If the number of iterations in \code{data} is too small,
+ an error message is printed indicating the minimum length of
+ pilot run. The minimum length is the required sample size for a
+ chain with no correlation between consecutive samples. Positive
+ autocorrelation will increase the required sample size above this
+ minimum value. An estimate \code{I} (the `dependence factor') of the
+ extent to which autocorrelation inflates the required sample size
+ is also provided. Values of \code{I} larger than 5 indicate strong
+ autocorrelation which may be due to a poor choice of starting value,
+ high posterior correlations or `stickiness' of the MCMC algorithm.
+
+ The number of `burn in' iterations to be discarded at the beginning
+ of the chain is also calculated.
+}
+
+\value{
+ A list with class \code{raftery.diag}. A print method is available
+ for objects of this class. the contents of the list are
+ \item{tspar}{The time series parameters of \code{data}}
+ \item{params}{A vector containing the parameters \code{r}, \code{s}
+ and \code{q}}
+ \item{Niters}{The number of iterations in \code{data}}
+ \item{resmatrix}{A 3-d array containing the results: \eqn{M} the
+ length of "burn in", \eqn{N} the required sample size, \eqn{Nmin}
+ the minimum sample size based on zero autocorrelation and
+ \eqn{I = (M+N)/Nmin} the "dependence factor"}
+}
+
+\section{Theory}{
+ The estimated sample size for variable U is based on the process \eqn{Z_t
+ = d(U_t <= u)} where \eqn{d} is the indicator function and u is the
+ qth quantile of U. The process \eqn{Z_t} is derived from the Markov
+ chain \code{data} by marginalization and truncation, but is not itself
+ a Markov chain. However, \eqn{Z_t} may behave as a Markov chain if
+ it is sufficiently thinned out. \code{raftery.diag} calculates the
+ smallest value of thinning interval \eqn{k} which makes the thinned
+ chain \eqn{Z^k_t} behave as a Markov chain. The required sample size is
+ calculated from this thinned sequence. Since some data is `thrown away'
+ the sample size estimates are conservative.
+
+ The criterion for the number of `burn in' iterations \eqn{m} to be
+ discarded, is that the conditional distribution of \eqn{Z^k_m}
+ given \eqn{Z_0} should be within \code{converge.eps} of the equilibrium
+ distribution of the chain \eqn{Z^k_t}.
+}
+
+\note{
+ \code{raftery.diag} is based on the FORTRAN program `gibbsit'
+ written by Steven Lewis, and available from the Statlib archive.
+}
+
+\references{
+ Raftery, A.E. and Lewis, S.M. (1992). One long run with diagnostics:
+ Implementation strategies for Markov chain Monte Carlo.
+ \emph{Statistical Science}, \bold{7}, 493-497.
+
+ Raftery, A.E. and Lewis, S.M. (1995). The number of iterations,
+ convergence diagnostics and generic Metropolis algorithms. \emph{In}
+ Practical Markov Chain Monte Carlo (W.R. Gilks, D.J. Spiegelhalter
+ and S. Richardson, eds.). London, U.K.: Chapman and Hall.
+}
+
+\keyword{htest}
diff --git a/man/read.and.check.Rd b/man/read.and.check.Rd
new file mode 100644
index 0000000..5d732b6
--- /dev/null
+++ b/man/read.and.check.Rd
@@ -0,0 +1,44 @@
+\name{read.and.check}
+\alias{read.and.check}
+\title{Read data interactively and check that it satisfies conditions}
+
+\usage{
+read.and.check(message = "", what = numeric(), lower, upper, answer.in,
+default)
+}
+
+\arguments{
+ \item{message}{message displayed before prompting for user input.}
+ \item{what}{the type of \code{what} gives the type of data to be read.}
+ \item{lower}{lower limit of input, for numeric input only.}
+ \item{upper}{upper limit of input, for numeric input only.}
+ \item{answer.in}{the input must correspond to one of the elements of the
+ vector \code{answer.in}, if supplied.}
+ \item{default}{value assumed if user enters a blank line.}
+}
+
+\description{
+ Input is read interactively and checked against conditions specified
+ by the arguments \code{what}, \code{lower}, \code{upper} and
+ \code{answer.in}. If the input does not satisfy all the conditions,
+ an appropriate error message is produced and the user is prompted
+ to provide input. This process is repeated until a valid input value
+ is entered.
+}
+
+\value{
+ The value of the valid input. When the \code{default} argument is
+ specified, a blank line is accepted as valid input and in this case
+ \code{read.and.check} returns the value of \code{default}.
+}
+
+\note{
+ Since the function does not return a value until it receives valid
+ input, it extensively checks the conditions for consistency before
+ prompting the user for input. Inconsistent conditions will cause
+ an error.
+}
+
+\author{Martyn Plummer}
+
+\keyword{utilities}
diff --git a/man/read.coda.Rd b/man/read.coda.Rd
new file mode 100644
index 0000000..6e97b9e
--- /dev/null
+++ b/man/read.coda.Rd
@@ -0,0 +1,52 @@
+\name{read.coda}
+\alias{read.coda}
+\alias{read.jags}
+\title{Read output files in CODA format}
+
+\usage{
+read.coda(output.file, index.file, start, end, thin, quiet=FALSE)
+read.jags(file = "jags.out", start, end, thin, quiet=FALSE)
+}
+
+\arguments{
+ \item{output.file}{The name of the file containing the monitored
+ output}
+ \item{index.file}{The name of the file containing the index, showing
+ which rows of the output file correspond to which variables}
+ \item{file}{For JAGS output, the name of the output file. The
+ extension ".out" may be omitted. There must be a corresponding
+ ".ind" file with the same file stem.}
+ \item{start}{First iteration of chain}
+ \item{end}{Last iteration of chain}
+ \item{thin}{Thinning interval for chain}
+ \item{quiet}{Logical flag. If true, a progress summary will be printed}
+}
+
+\description{
+ \code{read.coda} reads Markov Chain Monte Carlo output in
+ the CODA format produced by OpenBUGS and JAGS. By default, all
+ of the data in the file is read, but the arguments \code{start},
+ \code{end} and \code{thin} may be used to read a subset of the
+ data. If the arguments given to \code{start}, \code{end} or
+ \code{thin} are incompatible with the data, they are ignored.
+}
+
+\value{
+ An object of class \code{mcmc} containing a representation of
+ the data in the file.
+}
+
+\references{
+ Spiegelhalter DJ, Thomas A, Best NG and Gilks WR (1995).
+ \emph{BUGS: Bayesian inference Using Gibbs Sampling, Version 0.50.}
+ MRC Biostatistics Unit, Cambridge.
+}
+
+\author{Karen Vines, Martyn Plummer}
+
+\seealso{
+ \code{\link{mcmc}},
+ \code{\link{read.coda.interactive}},
+ \code{\link{read.openbugs}}.
+}
+\keyword{file}
diff --git a/man/read.coda.interactive.Rd b/man/read.coda.interactive.Rd
new file mode 100644
index 0000000..d8d020e
--- /dev/null
+++ b/man/read.coda.interactive.Rd
@@ -0,0 +1,39 @@
+\name{read.coda.interactive}
+\alias{read.coda.interactive}
+\title{Read CODA output files interactively}
+
+\usage{
+read.coda.interactive()
+}
+
+\description{
+ \code{read.coda.interactive} reads Markov Chain Monte Carlo output
+ in the format produced by the classic BUGS program. No arguments are
+ required. Instead, the user is prompted for the required information.
+}
+
+\value{
+ An object of class \code{mcmc.list} containing a representation of
+ the data in one or more BUGS output files.
+}
+
+\references{
+ Spiegelhalter DJ, Thomas A, Best NG and Gilks WR (1995).
+ \emph{BUGS: Bayesian inference Using Gibbs Sampling, Version 0.50.}
+ MRC Biostatistics Unit, Cambridge.
+}
+
+\note{
+ This function is normally called by the \code{codamenu} function,
+ but can also be used on a stand-alone basis.
+}
+
+\author{Nicky Best, Martyn Plummer}
+
+\seealso{
+ \code{\link{mcmc}},
+ \code{\link{mcmc.list}},
+ \code{\link{read.coda}},
+ \code{\link{codamenu}}.
+}
+\keyword{file}
diff --git a/man/read.openbugs.Rd b/man/read.openbugs.Rd
new file mode 100644
index 0000000..5000f5b
--- /dev/null
+++ b/man/read.openbugs.Rd
@@ -0,0 +1,44 @@
+\name{read.openbugs}
+\alias{read.openbugs}
+\title{Read CODA output files produced by OpenBUGS}
+
+\usage{
+read.openbugs(stem="", start, end, thin, quiet=FALSE)
+}
+
+\arguments{
+ \item{stem}{Character string giving the stem for the output files.
+ OpenBUGS produces files with names "<stem>CODAindex.txt",
+ "<stem>CODAchain1.txt", "<stem>CODAchain2.txt", ...}
+ \item{start}{First iteration of chain}
+ \item{end}{Last iteration of chain}
+ \item{thin}{Thinning interval for chain}
+ \item{quiet}{Logical flag. If true, a progress summary will be printed}
+}
+
+\description{
+ \code{read.openbugs} reads Markov Chain Monte Carlo output in
+ the CODA format produced by OpenBUGS.
+
+ This is a convenience wrapper around the function \code{read.coda}
+ which allows you to read all the data output by OpenBUGS by
+ specifying only the file stem.
+}
+
+\value{
+ An object of class \code{mcmc.list} containing output from all
+ chains.
+}
+
+\references{
+ Spiegelhalter DJ, Thomas A, Best NG and Lunn D (2004).
+ \emph{WinBUGS User Manual, Version 2.0, June 2004},
+ MRC Biostatistics Unit, Cambridge.
+}
+
+\author{Martyn Plummer}
+
+\seealso{
+ \code{\link{read.coda}}.
+}
+\keyword{file}
diff --git a/man/spectrum0.Rd b/man/spectrum0.Rd
new file mode 100644
index 0000000..fb1ca46
--- /dev/null
+++ b/man/spectrum0.Rd
@@ -0,0 +1,63 @@
+\name{spectrum0}
+\alias{spectrum0}
+\title{Estimate spectral density at zero}
+\description{
+The spectral density at frequency zero is estimated by fitting a glm to
+the low-frequency end of the periodogram. \code{spectrum0(x)/length(x)}
+estimates the variance of \code{mean(x)}.
+
+}
+\usage{
+spectrum0(x, max.freq = 0.5, order = 1, max.length = 200)
+}
+\arguments{
+\item{x}{A time series.}
+\item{max.freq}{The glm is fitted on the frequency range (0, max.freq]}
+\item{order}{Order of the polynomial to fit to the periodogram.}
+
+\item{max.length}{The data \code{x} is aggregated if necessary by
+taking batch means so that the length of the series is less than
+\code{max.length}. If this is set to \code{NULL} no aggregation occurs.}
+}
+\details{
+The raw periodogram is calculated for the series \code{x} and a generalized
+linear model with family \code{Gamma} and log link is fitted to
+the periodogram.
+
+The linear predictor is a polynomial in terms of the frequency. The
+degree of the polynomial is determined by the parameter \code{order}.
+}
+\value{
+A list with the following values
+\item{spec}{The predicted value of the spectral density at frequency zero.}
+}
+\references{
+Heidelberger, P and Welch, P.D. A spectral method for confidence interval
+generation and run length control in simulations. Communications of the
+ACM, Vol 24, pp233-245, 1981.
+}
+\section{Theory}{
+Heidelberger and Welch (1991) observed that the usual non-parametric
+estimator of the spectral density, obtained by smoothing the periodogram,
+is not appropriate for frequency zero. They proposed an alternative
+parametric method which consisted of fitting a linear model to the
+log periodogram of the batched time series. Some technical problems
+with model fitting in their original proposal can be overcome by using
+a generalized linear model.
+
+Batching of the data, originally proposed in order to save space, has the
+side effect of flattening the spectral density and making a polynomial
+fit more reasonable. Fitting a polynomial of degree zero is equivalent
+to using the `batched means' method.
+}
+\note{
+The definition of the spectral density used here differs from that used by
+\code{spec.pgram}. We consider the frequency range to be between 0 and 0.5,
+not between 0 and \code{frequency(x)/2}.
+
+The model fitting may fail on chains with very high autocorrelation.
+}
+\seealso{
+ \code{\link{spectrum}}, \code{\link{spectrum0.ar}}, \code{\link{glm}}.
+}
+\keyword{ts}
diff --git a/man/spectrum0.ar.Rd b/man/spectrum0.ar.Rd
new file mode 100644
index 0000000..2fcdeeb
--- /dev/null
+++ b/man/spectrum0.ar.Rd
@@ -0,0 +1,34 @@
+\name{spectrum0.ar}
+\alias{spectrum0.ar}
+\title{Estimate spectral density at zero}
+\description{
+The spectral density at frequency zero is estimated by fitting an
+autoregressive model. \code{spectrum0(x)/length(x)} estimates the
+variance of \code{mean(x)}.
+}
+\usage{
+spectrum0.ar(x)
+}
+\arguments{
+\item{x}{A time series.}
+}
+\details{
+The \code{ar()} function to fit an autoregressive model to the time
+series x. For multivariate time series, separate models are fitted for
+each column. The value of the spectral density at zero is then given
+by a well-known formula.
+}
+\value{
+A list with the following values
+\item{spec}{The predicted value of the spectral density at frequency zero.}
+\item{order}{The order of the fitted model}
+}
+\note{
+The definition of the spectral density used here differs from that used by
+\code{spec.pgram}. We consider the frequency range to be between 0 and 0.5,
+not between 0 and \code{frequency(x)/2}.
+}
+\seealso{
+ \code{\link{spectrum}}, \code{\link{spectrum0}}, \code{\link{glm}}.
+}
+\keyword{ts}
diff --git a/man/summary.mcmc.Rd b/man/summary.mcmc.Rd
new file mode 100644
index 0000000..47b3784
--- /dev/null
+++ b/man/summary.mcmc.Rd
@@ -0,0 +1,36 @@
+\name{summary.mcmc}
+\alias{summary.mcmc}
+\alias{summary.mcmc.list}
+%\alias{print.summary.mcmc}
+%\alias{print.summary.mcmc.list}
+\title{Summary statistics for Markov Chain Monte Carlo chains}
+
+\usage{
+\method{summary}{mcmc}(object, quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975), ...)
+}
+
+\arguments{
+ \item{object}{an object of class \code{mcmc} or \code{mcmc.list}}
+ \item{quantiles}{a vector of quantiles to evaluate for each variable}
+ \item{...}{a list of further arguments}
+}
+
+\description{
+ \code{summary.mcmc} produces two sets of summary statistics for
+ each variable:
+
+ Mean, standard deviation, naive standard error of the mean
+ (ignoring autocorrelation of the chain) and time-series standard
+ error based on an estimate of the spectral density at 0.
+
+ Quantiles of the sample distribution using the \code{quantiles}
+ argument.
+}
+
+\author{Martyn Plummer}
+
+\seealso{
+ \code{\link{mcmc}},
+ \code{\link{mcmc.list}}.
+}
+\keyword{univar}
diff --git a/man/thin.Rd b/man/thin.Rd
new file mode 100644
index 0000000..e1d33f2
--- /dev/null
+++ b/man/thin.Rd
@@ -0,0 +1,28 @@
+\name{thin}
+\alias{thin}
+\title{Thinning interval}
+
+\usage{
+thin(x, ...)
+}
+
+\arguments{
+ \item{x}{a regular time series}
+ \item{...}{a list of arguments}
+ }
+
+\description{
+ \code{thin} returns the interval between successive
+ values of a time series. \code{thin(x)} is equivalent
+ to \code{1/frequency(x)}.
+
+ This is a generic function. Methods have been implemented
+ for \code{mcmc} objects.
+}
+
+\author{Martyn Plummer}
+
+\seealso{
+ \code{\link{time}}.
+}
+\keyword{ts}
diff --git a/man/time.mcmc.Rd b/man/time.mcmc.Rd
new file mode 100644
index 0000000..4339e33
--- /dev/null
+++ b/man/time.mcmc.Rd
@@ -0,0 +1,36 @@
+\name{time.mcmc}
+\alias{time.mcmc}
+\alias{start.mcmc}
+\alias{end.mcmc}
+\alias{frequency.mcmc}
+\alias{thin.mcmc}
+\alias{time.mcmc.list}
+\alias{start.mcmc.list}
+\alias{end.mcmc.list}
+\alias{thin.mcmc.list}
+\title{Time attributes for mcmc objects}
+
+\usage{
+\method{time}{mcmc}(x, ...)
+\method{start}{mcmc}(x, ...)
+\method{end}{mcmc}(x, ...)
+\method{thin}{mcmc}(x, ...)
+}
+
+\arguments{
+ \item{x}{an \code{mcmc} or \code{mcmc.list} object}
+ \item{...}{extra arguments for future methods}
+}
+
+\description{
+ These are methods for mcmc objects for the generic time
+ series functions.
+}
+
+\seealso{
+ \code{\link{time}},
+ \code{\link{start}},
+ \code{\link{frequency}},
+ \code{\link{thin}}.
+}
+\keyword{ts}
diff --git a/man/traceplot.Rd b/man/traceplot.Rd
new file mode 100644
index 0000000..38a1db7
--- /dev/null
+++ b/man/traceplot.Rd
@@ -0,0 +1,29 @@
+\name{traceplot}
+\alias{traceplot}
+\title{Trace plot of MCMC output}
+
+\usage{traceplot(x, smooth = TRUE, col, type, ylab, ...)}
+
+\arguments{
+ \item{x}{An \code{mcmc} or \code{mcmc.list} object}
+ \item{smooth}{draw smooth line through trace plot}
+ \item{col}{graphical parameter (see \code{par})}
+ \item{type}{graphical parameter (see \code{par})}
+ \item{ylab}{graphical parameter (see \code{par})}
+ \item{...}{further graphical parameters}
+}
+
+\description{
+Displays a plot of iterations \emph{vs.} sampled values for each variable
+in the chain, with a separate plot per variable.
+}
+
+\note{
+You can call this function directly, but it is more usually
+called by the \code{plot.mcmc} function.
+}
+
+\seealso{
+ \code{\link{densplot}}, \code{\link{plot.mcmc}}.
+}
+\keyword{hplot}
diff --git a/man/varnames.Rd b/man/varnames.Rd
new file mode 100644
index 0000000..74da4a1
--- /dev/null
+++ b/man/varnames.Rd
@@ -0,0 +1,37 @@
+\name{varnames}
+\alias{varnames}
+\alias{varnames<-}
+\alias{chanames}
+\alias{chanames<-}
+\title{Named dimensions of MCMC objects}
+
+\usage{
+varnames(x, allow.null=TRUE)
+chanames(x, allow.null=TRUE)
+varnames(x) <- value
+chanames(x) <- value
+}
+
+\arguments{
+\item{x}{an \code{mcmc} or \code{mcmc.list} object}
+\item{allow.null}{Logical argument that determines whether the function may return NULL}
+\item{value}{A character vector, or NULL}
+}
+
+\value{
+ A character vector , or NULL.
+}
+
+\description{
+ \code{varnames()} returns the variable names and \code{chanames}
+ returns the chain names, or NULL if these are not set.
+
+ If \code{allow.null = FALSE} then \code{NULL} values will be
+ replaced with canonical names.
+}
+
+\seealso{
+ \code{\link{mcmc}},
+ \code{\link{mcmc.list}}.
+}
+\keyword{manip}
diff --git a/man/window.mcmc.Rd b/man/window.mcmc.Rd
new file mode 100644
index 0000000..a2e589d
--- /dev/null
+++ b/man/window.mcmc.Rd
@@ -0,0 +1,36 @@
+\name{window.mcmc}
+\alias{window.mcmc}
+\alias{window.mcmc.list}
+\title{Time windows for mcmc objects}
+
+\usage{
+\method{window}{mcmc}(x, start, end, thin, ...)
+}
+
+\arguments{
+ \item{x}{an mcmc object}
+ \item{start}{the first iteration of interest}
+ \item{end}{the last iteration of interest}
+ \item{thin}{the required interval between successive samples}
+ \item{...}{futher arguments for future methods}
+ }
+
+\description{
+ \code{window.mcmc} is a method for \code{mcmc} objects which is
+ normally called by the generic function \code{window}
+
+ In addition to the generic parameters, \code{start} and \code{end}
+ the additional parameter \code{thin} may be used to thin out the
+ Markov chain. Setting thin=k selects every kth iteration starting
+ with the first. Note that the value of \code{thin} is \emph{absolute}
+ not relative. The value supplied given to the parameter \code{thin}
+ must be a multiple of \code{thin(x)}.
+
+ Values of \code{start}, \code{end} and \code{thin} which are inconsistent
+ with \code{x} are ignored, but a warning message is issued.
+}
+
+\seealso{
+ \code{\link{window}}, \code{\link{thin}}.
+}
+\keyword{ts}
diff --git a/misc/change.tfoption.Rd b/misc/change.tfoption.Rd
new file mode 100644
index 0000000..7e91728
--- /dev/null
+++ b/misc/change.tfoption.Rd
@@ -0,0 +1,29 @@
+\name{change.tfoption}
+\alias{change.tfoption}
+\title{Interactively modify logical codamenu options.}
+
+\usage{
+change.tfoption(string, option)
+}
+
+\arguments{
+\item{string}{Prompt string to print}
+\item{option}{String giving the name of an option to change}
+}
+
+\description{
+Interactively changes logical codamenu options using the
+\code{read.yesno} function.
+}
+
+\note{
+This is a utility function for the codamenu system, and is not
+meant to be called directly by the user.
+}
+
+\seealso{
+ \code{\link{coda.options}},
+ \code{\link{read.yesno}}.
+}
+
+\keyword{utilities}
diff --git a/misc/codamenu2.Rd b/misc/codamenu2.Rd
new file mode 100644
index 0000000..ee0c64f
--- /dev/null
+++ b/misc/codamenu2.Rd
@@ -0,0 +1,40 @@
+\name{codamenu2}
+\alias{codamenu.main}
+\alias{codamenu.anal}
+\alias{codamenu.diags}
+\alias{codamenu.diags.autocorr}
+\alias{codamenu.diags.crosscorr}
+\alias{codamenu.diags.heidel}
+\alias{codamenu.diags.raftery}
+\alias{codamenu.diags.gelman}
+\alias{codamenu.diags.geweke}
+\alias{codamenu.options}
+\alias{codamenu.options.data}
+\alias{codamenu.options.diag}
+\alias{codamenu.options.gelman}
+\alias{codamenu.options.geweke.bin}
+\alias{codamenu.options.geweke.win}
+\alias{codamenu.options.heidel}
+\alias{codamenu.options.plot}
+\alias{codamenu.options.plot.kernel}
+\alias{codamenu.options.plot.ps}
+\alias{codamenu.options.raftery}
+\alias{codamenu.options.stats}
+\alias{tidy.up}
+\alias{codamenu.ps}
+\alias{codamenu.output.header}
+
+\title{Submenus for the coda package}
+
+\description{
+These functions are not meant to be called directly by the user, but
+should be accessed via the \code{codamenu} function.
+
+All submenus are called from the \code{codamenu.main} function.
+Each submenu returns a string giving the name of the next submenu to be
+called, or "quit" if the user decides to quit. This system avoids the
+recursive calling of submenus which would accumulate `dead' environments.
+}
+
+\author{Martyn Plummer}
+\keyword{utilities}
diff --git a/misc/read.yesno.Rd b/misc/read.yesno.Rd
new file mode 100644
index 0000000..1542e79
--- /dev/null
+++ b/misc/read.yesno.Rd
@@ -0,0 +1,29 @@
+\name{read.yesno}
+\alias{read.yesno}
+\title{Read logical values from the command line}
+
+\usage{
+ read.yesno(string, default=TRUE)
+}
+
+\arguments{
+ \item{string}{A string to prompt the user with}
+ \item{default}{Logical value giving the default answer}
+}
+
+\value{
+ A single logical value.
+}
+
+\description{
+ This is a utility function for interactively reading a logical
+ value from the command line.
+
+ The prompt string is printed, along a hint about the default answer.
+ Answers "yes", "YES", "no", "NO" or any partial match are accepted.
+ If the user enters a blank line, the default value is returned.
+}
+
+\author{Martyn Plummer}
+
+\keyword{utilities}
diff --git a/misc/set.mfrow.Rd b/misc/set.mfrow.Rd
new file mode 100644
index 0000000..57b2ee0
--- /dev/null
+++ b/misc/set.mfrow.Rd
@@ -0,0 +1,26 @@
+\name{set.mfrow}
+\alias{set.mfrow}
+\title{Set layout for graphical plots}
+
+\usage{
+ set.mfrow (Nchains = 1, Nparms = 1, nplots = 1, sepplot = FALSE)
+}
+
+\arguments{
+ \item{Nchains}{Number of chains to plot}
+ \item{Nparms}{Number of variables in each chain}
+ \item{nplots}{Number of plots to produce}
+ \item{sepplot}{Logical value controlling whether separate plots are created for each chain}
+}
+
+\description{
+ This utility function is used by plotting functions in the coda package
+ to determine the default layout, unless overriden by the \code{user.layout}
+ option.
+}
+
+\note{
+ This is not meant to be called directly by the user.
+}
+
+\keyword{dplot}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-coda.git
More information about the debian-science-commits
mailing list