[r-cran-mcmcpack] 60/90: Imported Upstream version 1.1-4
Andreas Tille
tille at debian.org
Fri Dec 16 09:07:48 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-mcmcpack.
commit cb749f54df28793678358b54a2083cf981248da8
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 08:07:21 2016 +0100
Imported Upstream version 1.1-4
---
COPYING | 278 -----------------------------------------------
DESCRIPTION | 10 +-
HISTORY | 16 +++
MD5 | 175 +++++++++++++++++++++++++++++
R/MCMCbinaryChange.R | 112 ++++++++++---------
R/MCMCoprobit.R | 2 +-
R/MCMCpoissonChange.R | 177 ++++++++++++++++--------------
R/MCMCprobitChange.R | 150 +++++++++++++------------
data/PErisk.rda | Bin 1968 -> 4653 bytes
data/SupremeCourt.rda | Bin 471 -> 3896 bytes
man/MCMCbinaryChange.Rd | 3 +-
man/MCMCoprobit.Rd | 8 +-
man/MCMCpoissonChange.Rd | 5 +-
man/MCMCprobitChange.Rd | 6 +-
man/plotChangepoint.Rd | 0
man/plotState.Rd | 0
src/MCMCbinaryChange.cc | 34 +++---
src/MCMChlogit.cc | 2 +-
src/MCMChpoisson.cc | 2 +-
src/MCMChregress.cc | 6 +-
src/MCMClogit.cc | 11 +-
src/MCMCoprobit.cc | 15 +--
src/MCMCoprobitChange.cc | 46 ++++----
src/MCMCpoisson.cc | 141 ++++++++++++------------
src/MCMCpoissonChange.cc | 18 ++-
src/MCMCprobit.cc | 51 +++++----
src/MCMCprobitChange.cc | 22 ++--
src/MCMCregress.cc | 88 +++++++++------
28 files changed, 685 insertions(+), 693 deletions(-)
diff --git a/COPYING b/COPYING
deleted file mode 100644
index 727ef8f..0000000
--- a/COPYING
+++ /dev/null
@@ -1,278 +0,0 @@
- GNU GENERAL PUBLIC LICENSE
- Version 2, June 1991
-
- Copyright (C) 1989, 1991 Free Software Foundation, Inc.
- 675 Mass Ave, Cambridge, MA 02139, 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.
diff --git a/DESCRIPTION b/DESCRIPTION
index 87286f6..aa4033b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,9 +1,9 @@
Package: MCMCpack
-Version: 1.1-1
-Date: 2011-7-8
+Version: 1.1-4
+Date: 2011-9-21
Title: Markov chain Monte Carlo (MCMC) Package
Author: Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
-Maintainer: Andrew D. Martin <admartin at wustl.edu>
+Maintainer: Jong Hee Park <jhp at uchicago.edu>
Depends: R (>= 2.10.0), coda (>= 0.11-3), MASS, stats
Description: This package contains functions to perform Bayesian
inference using posterior simulation for a number of
@@ -18,6 +18,6 @@ Description: This package contains functions to perform Bayesian
License: GPL-2
SystemRequirements: gcc (>= 4.0)
URL: http://mcmcpack.wustl.edu
-Packaged: 2011-07-08 17:10:47 UTC; adm
+Packaged: 2011-09-26 13:42:11 UTC; jongheepark
Repository: CRAN
-Date/Publication: 2011-07-23 13:02:20
+Date/Publication: 2011-09-26 14:23:46
diff --git a/HISTORY b/HISTORY
index f105cc2..8ec5fba 100644
--- a/HISTORY
+++ b/HISTORY
@@ -1,6 +1,22 @@
//
// Changes and Bug Fixes
//
+1.1-3 to 1.1-4
+ * fixed an error in MCMChregress: NQ has to be replaced by NP (line 223)
+ * in MCMChregress, MCMChpoisson, MCMChlogit, R_CheckUserInterrupt(void) is replaced by R_CheckUserInterrupt()
+ * fixed an error in the marginal likelihood computation
+ in MCMCproit.cc and MCMCregress.cc (thanks to Sid Chib)
+ * changed B0 as prior precision in MCMCprobitChange, MCMCoprobitChange and MCMCpoissonChange
+ * the acceptance rate is reported if verbose > 0
+ * added the reporting of marginal likelihood components if verbose > 0
+
+1.0-11 to 1.1-2
+ * fixed an error in MCMCoprobit.R
+ * changed the documentation of MCMCoprobit.
+ * Added m=0 option to all changepoint models. Now all changepoint models allow users
+ to set m=0 for marginal likelihood computation.
+ * Fixed a bug in MCMCbinaryChange. I don't know why but log(dbeta) in Scythe often
+ generates NaN results while lndbeta1() does not.
1.0-10 to 1.0-11
* updated CITATION file in inst/ to coincide with JStatSoft publication
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..ae0a22f
--- /dev/null
+++ b/MD5
@@ -0,0 +1,175 @@
+c1f2fefa751f851fe593cc8bfed0c466 *DESCRIPTION
+510bf0941463c026351c1697e7516b8b *HISTORY
+c795674882accfc42204f5c71cb40180 *LICENSE
+3910bb9b28917e4e2473a2c83c4d24fa *NAMESPACE
+1eaf24e8d468ac0b624a9aa90cded306 *R/BayesFactors.R
+626c7d968d4f56fe081b595cd59d4976 *R/MCMCSVDreg.R
+c93b534213a07cec26b71eaa48f02295 *R/MCMCbinaryChange.R
+2736b2f4759f92e4e1d49d92fd0811e8 *R/MCMCdynamicEI.R
+87107048db50664880297536b41ce8e7 *R/MCMCdynamicIRT1d-b.R
+b0810a4be00979118db6e89ea9f2be46 *R/MCMCdynamicIRT1d.R
+f99ecb9f321a578b1cdfa9e5b6e163c6 *R/MCMCfactanal.R
+dcd83d013877c39ace4ca33c4357c779 *R/MCMChierBetaBinom.R
+065bd9bb18b1bc5d6011137361dc3412 *R/MCMChierEI.R
+4b95a74aa87e670a88106d4661b6e833 *R/MCMChlogit.R
+e2b64d567aa2e22410251a753f0b55c6 *R/MCMChpoisson.R
+0a759c18d9f3b280dbe6090c92e20e35 *R/MCMChregress.R
+a4db065464be685b498d341ecc3c2109 *R/MCMCirt1d.R
+78b6d6361c4e9b3056226d66c50a3479 *R/MCMCirtHier1d.R
+4b433d0cace650f7304ff9a94e69ad4a *R/MCMCirtKd.R
+7582e3ad855f5ece89f1dfbcf4a2a856 *R/MCMCirtKdHet.R
+0293dbca0de9cb96b203f1739cda2c02 *R/MCMCirtKdRob.R
+ba21d6d4a4dcaaa56b6784b749bcbc07 *R/MCMClogit.R
+7a83fa7762a8e777a5cacc21c14f2498 *R/MCMCmetrop1R.R
+33f53d5be86c9d714dd7b68a27285c2f *R/MCMCmixfactanal.R
+a6508dd38804705e750995ca27070cde *R/MCMCmnl.R
+4b7cb2b7fa8eb3d0435c55dfdcbed40e *R/MCMCoprobit.R
+6b63259803b2890eae6d5778592ed7dd *R/MCMCoprobitChange.R
+be0928d3cd74e5485b1e26a9f0a15c59 *R/MCMCordfactanal.R
+d438fc9dabd9993da5d4a45a3df68c39 *R/MCMCpoisson.R
+5bf87e428e3835f6f73d13fc432b36e5 *R/MCMCpoissonChange.R
+b98bdefca76d3316c664bcb2ba2bcc38 *R/MCMCprobit.R
+31fdcb514ab53ba258a581536fcbda60 *R/MCMCprobitChange.R
+e27dcf535f23c25018bb778b61640f22 *R/MCMCquantreg.R
+4b0822eb63ca0687af1c99ddc8f9c01d *R/MCMCregress.R
+720834bbb59d6c4ce7c610acd22558be *R/MCMCtobit.R
+4d2218bcdb4715ac9d791f26d5469f9b *R/MCmodels.R
+3a98b4d919a265b46c29d14310bb90d4 *R/SSVSquantreg.R
+555c87acdcc2f8572b14151e29022e20 *R/SSVSquantregsummary.R
+fdfcfe9909cadaf15c47bbb604be2257 *R/automate.R
+188a57ae832e3cc3bd447271d03a6996 *R/btsutil.R
+5aca34f825d2ec38f932094b6bbe601b *R/distn.R
+18daa9fee800c7d8d5827405738531e2 *R/hidden-hmodels.R
+0f904468b8ce8fc5c3869f87d11da905 *R/hidden.R
+f799711e5d3fb7140cbc0995360339b3 *R/procrust.R
+2e00b2b1593d22b5d1c579e92b294106 *R/scythe.R
+c67ac5238de55cdc8c49d017e0226ed9 *R/tomog.R
+6377da3453982f578df4a886f673122c *R/utility.R
+a4b896c744c34c6cece2779d541b538d *R/zzz.R
+e42f762c27e26e0ace372c03c47ac4e5 *README
+46e955baeb55718a380a11fb34d2ea67 *cleanup
+e09d1c5ac6250f832af90566a9e565fb *configure
+21e1a006fb388ed39daf16e6d2ec9549 *configure.ac
+e17202d22ce5ab5ab7b01e3b248a4008 *data/Nethvote.rda
+bfb9323ec59f28654e261f570bd1fb9a *data/PErisk.rda
+93e636ac58b966f9a74315f51be07d25 *data/Rehnquist.rda
+accb7b5b46207a02e49de2702a6faff4 *data/Senate.rda
+f5d5cc1cd9f9f1699e249946f8767179 *data/SupremeCourt.rda
+a5cfcf73e21c03eeaf4f3baa5c492c14 *inst/CITATION
+2b42855343b5f0c51c7e9411aef57509 *man/BayesFactor.Rd
+293ff17f852d0e60dd1aac1f419353d2 *man/MCMCSVDreg.Rd
+d479bf1112276aef579b481f05eb3a3f *man/MCMCbinaryChange.Rd
+5078ab3b9d86c7165909791acf1fc886 *man/MCMCdynamicEI.Rd
+102a94ced183cf023da7696f948bb82f *man/MCMCdynamicIRT1d.Rd
+29e89e27c6b11794c5dc6f50618375e0 *man/MCMCfactanal.Rd
+5dabac071a3722205135ae2cbae0b3d9 *man/MCMChierEI.Rd
+00a29263d43f0e2006de60d4800f5f7f *man/MCMChlogit.Rd
+b894aa07a0e68f9ffdeb3b09bab0aee2 *man/MCMChpoisson.Rd
+8d9c22e8d965d60581ec837803790c80 *man/MCMChregress.Rd
+52f2cb9bc1cde5e38d35aa5d9858be4a *man/MCMCirt1d.Rd
+336682b062ff43261f0e75a4b51a34d8 *man/MCMCirtHier1d.Rd
+0d6b55f1f6062483b894fd4f4e9cecfd *man/MCMCirtKd.Rd
+19142a1e2a1a9217001a2d19b377b4ce *man/MCMCirtKdHet.Rd
+1bcb750570641123acbbc149bdad93af *man/MCMCirtKdRob.Rd
+1ee97a34364e07ec83542d78fce4d823 *man/MCMClogit.Rd
+8c4a033b4e5fca56699e9cbd96a56222 *man/MCMCmetrop1R.Rd
+751e56edabc12da2d8246f7c34cd7613 *man/MCMCmixfactanal.Rd
+b7dfcf274126008379b0ee903f134b25 *man/MCMCmnl.Rd
+cf07e2c9f307215c8fe841f5381b62f8 *man/MCMCoprobit.Rd
+ac737a1e0cc404c1ac8dc91e2e870725 *man/MCMCoprobitChange.Rd
+a491294e97dd0b7aa5d0de525542ea8f *man/MCMCordfactanal.Rd
+03f2e93e6876a4f0e87b3469bfe38d2f *man/MCMCpoisson.Rd
+3d38651d02ff76f5d7430bdccf56c1d2 *man/MCMCpoissonChange.Rd
+1917e296040b80c92525e5cb8ad02e71 *man/MCMCprobit.Rd
+0b79ed083594613481da272097d0d359 *man/MCMCprobitChange.Rd
+2f9caf8983e405f912eb5b25f29784eb *man/MCMCquantreg.Rd
+df339d8d9dce904214b8d3d120e88182 *man/MCMCregress.Rd
+c314e36afdd60449cf4f8ec08e00f80d *man/MCMCtobit.Rd
+a96eacae6b64827e40505661fb636cd4 *man/MCbinomialbeta.Rd
+193401fb83df78594bb965c7bbac968f *man/MCmultinomdirichlet.Rd
+4f2164ebeaaf83bde6a7f6c78b718d32 *man/MCnormalnormal.Rd
+ecc40a241e849aacfa5ee3f17b8dd97d *man/MCpoissongamma.Rd
+c04a5372b5a75cf67ef78d7e94d121e9 *man/Nethvote.Rd
+1aae2198de60fb0e4d568b0f2f4d994c *man/PErisk.Rd
+034a4e9a54b8eff7dddcee5d0a923a02 *man/PostProbMod.Rd
+aa99f8aa2b507bc8e1465a2b6e3993bd *man/QRSSVSplot.Rd
+594c242058a0e3f9f2b8b5d4067c5b72 *man/QRSSVSsummary.Rd
+3ad4d7cb7b7b9fea628a16b8a9ab4053 *man/Rehnquist.Rd
+43eae86dc52b891ed0a9fc9d1c6887bb *man/SSVSquantreg.Rd
+e6849fc07eb2dd45dcc07300e7cd149d *man/Senate.Rd
+03bbb0adb292bf12b04fe466816d20a6 *man/SupremeCourt.Rd
+1d8a9684290a7afeb7202f3797f5da71 *man/choicevar.Rd
+50b0fec7baf798a609432b1578f65b41 *man/dirichlet.Rd
+5203f2019ccb8b2dbf032c1689335ec3 *man/dtomog.Rd
+b6f23a0e12fa9bb16f571261b265921f *man/invgamma.Rd
+0051db062affc35be76c3ef5591dac82 *man/iwishart.Rd
+629023d17bc9299449eec0b74a73340d *man/mptable.Rd
+96ab1db8e6a5fb7c6640f06dd95a162c *man/noncenhypergeom.Rd
+170bc2168de62e08f7157a8cbc635e06 *man/plotChangepoint.Rd
+c64abc923f2c6f56a167a815c881ffb1 *man/plotState.Rd
+6a88d877e88b5abfd7a66d659c9b4e6a *man/procrust.Rd
+41feb17dda3b00cffa2f6176ba584b74 *man/readscythe.Rd
+254168bb5258fd00a7b766fb7d3bbfd9 *man/tomog.Rd
+f5ba8eb700c67f188505a4585dd39f6b *man/topmodels.Rd
+a4c560360ba912bb342a429e434bcc96 *man/vech.Rd
+8699162f8e39a03d96dd633866f8c4ee *man/wishart.Rd
+9d56915bf90694106df6586cec1bada0 *man/writescythe.Rd
+ee8f580a79a520bbdf458410b49aae2a *man/xpnd.Rd
+1c68060796310df44a049e66f10b14f7 *src/MCMCSVDreg.cc
+6474ba3567cea0fc10d61bbdc35c649a *src/MCMCbinaryChange.cc
+c6353855b63e402d7e60d248b22e6f50 *src/MCMCdynamicEI.cc
+583345a23401484a6ba5f0b0b206c408 *src/MCMCdynamicIRT1d-b.cc
+5caf9c1724b4e023483fe45ad247f2f2 *src/MCMCdynamicIRT1d.cc
+b3323f678776ada7c1fdf0f400be12e0 *src/MCMCfactanal.cc
+2745b72494da22ef8e0f62553011bbc9 *src/MCMCfcds.h
+8e5750dd27556aa3cc6a5ba4858889ca *src/MCMChierBetaBinom.cc
+d4f36358e358d81d2d432c8cf23cc33d *src/MCMChierEI.cc
+4221e7cbdb9fbd6a5b6ce16955e88fb5 *src/MCMChlogit.cc
+ab7fc8d513b7c66d9b777af3fda19133 *src/MCMChpoisson.cc
+ea2d9f60541c419765dde67824ca9139 *src/MCMChregress.cc
+fa7404ccc4ce053d13b2e0b9a86819b8 *src/MCMCirt1d.cc
+abdf98660af4bed11d89e04074d9f1d1 *src/MCMCirtHier1d.cc
+ac9b6d3f588b46cf009816bf028ee5ed *src/MCMCirtKdHet.cc
+01541fda7e59e9653da81371de726eb6 *src/MCMCirtKdRob.cc
+82710c136632704b70c71cc8ee5dca5c *src/MCMClogit.cc
+8e8ab807995eedcd108f11ad639475a9 *src/MCMClogituserprior.cc
+d886c33fcacd18ba4435f4171159824f *src/MCMCmetrop1R.cc
+af5765bb03517864b5d29b8c4b594ad2 *src/MCMCmixfactanal.cc
+eaab5804f229b660e958f6aa2cf728ba *src/MCMCmnl.h
+97c3e9b9e2d034ff9bd12fd1e1c81466 *src/MCMCmnlMH.cc
+9f0c3eea4d8ecafe8ab931c12d1616f9 *src/MCMCmnlslice.cc
+a78375184abf7c97bf7fbf3f7614cad3 *src/MCMCoprobit.cc
+387f517e5df0bde2b483c8986ce10983 *src/MCMCoprobitChange.cc
+6addbe81efb51cb4b8401e9a6931f782 *src/MCMCordfactanal.cc
+6bfed80bb7c109537461e6381837b6c4 *src/MCMCpoisson.cc
+628013b775615c38951a2d931e1c1fe4 *src/MCMCpoissonChange.cc
+84b94f6e67b026a915a23fb85bdaf00b *src/MCMCprobit.cc
+8a2f4a79c6bc951310074f5757041424 *src/MCMCprobitChange.cc
+a6a01c3906a8b92b9f1a12cfa92a1a63 *src/MCMCprobitres.cc
+806af936660056856d4b908f8fb0cfa4 *src/MCMCquantreg.cc
+32d7032c2fda6bdd5a8177fd9858ffe4 *src/MCMCregress.cc
+ef5566f824887d9bbef5b76f1b64d839 *src/MCMCrng.h
+f28dc81e54a6ca54a2202a0d6f52ed40 *src/MCMCtobit.cc
+2add70b59d5e52089b08311d62f5011e *src/Makevars
+30978b84a8789fd5a22a0c606b94c02f *src/Makevars.in
+6ba1ab7daac8f7149663686fd7abfd53 *src/SSVSquantreg.cc
+514000a986da5e497a7474e17379c626 *src/algorithm.h
+da3e44a3ae5138017ebd112e8e7af5c4 *src/datablock.h
+d489c387b86b46548f3f149fa6778734 *src/defs.h
+2acfadfc5d2ac2a532dc37e9c1f694b1 *src/distributions.h
+8e2a5e1e884ef5dd8486920ce3646f30 *src/error.h
+176c95d1fc3f3f7073c5c49df98a3abb *src/ide.h
+d9bbcd86bcee5ba918f3d368bd270712 *src/la.h
+ef7b879e29de6e8e11527700115686d6 *src/lapack.h
+cfc7f8d03a981a3453f81f2bd67d02c5 *src/lecuyer.cc
+12d602243cc945e693e56174f31874a9 *src/lecuyer.h
+9ccd363601ec9907b2b6226355445ada *src/matrix.h
+7ab16dbb6757a5de37dc7de9867a3c10 *src/matrix_bidirectional_iterator.h
+f924b218cd6eaed80846ed0cba534b03 *src/matrix_forward_iterator.h
+3dc8b8e303daa7ffdf3e2d29b2fa43d9 *src/matrix_random_access_iterator.h
+b0f0c5a5c5197df32d18b77af2b85957 *src/mersenne.h
+563652e8eec5ba3bd78ad6951a3ca27a *src/optimize.h
+1df60939513617beba99fe82eb7374c9 *src/rng.h
+b48baf77b93da895a04154ef15c8683f *src/rtmvnorm.h
+8aa95ee02342c8863d2d03b0e722aa4a *src/smath.h
+c3770569055925fec0945088672f5b61 *src/stat.h
diff --git a/R/MCMCbinaryChange.R b/R/MCMCbinaryChange.R
old mode 100644
new mode 100755
index b720031..23a26a6
--- a/R/MCMCbinaryChange.R
+++ b/R/MCMCbinaryChange.R
@@ -28,8 +28,6 @@
n <- nrow(y)
ns <- m+1
- ## prior for transition matrix
- A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
## get marginal likelihood argument
marginal.likelihood <- match.arg(marginal.likelihood)
@@ -41,60 +39,70 @@
chib <- 1
}
- ## starting values
- Pstart <- check.P(P.start, m=m, a=a, b=b)
- phistart <- check.theta(phi.start, ns, y, min=0, max=1)
-
nstore <- mcmc/thin
+ if (m == 0){
+ b0 <- c0/(c0 + d0)
+ B0 <- c0*d0/(c0 + d0)^2*(c0 + d0 + 1)
+ output <- MCMCprobit(y~1, burnin = burnin, mcmc = mcmc,
+ thin = thin, verbose = verbose, b0 = b0, B0 = B0,
+ marginal.likelihood = marginal.likelihood)
+ attr(output, "y") <- y
+ }
+ else {
+ ## prior for transition matrix
+ A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
+ Pstart <- check.P(P.start, m=m, a=a, b=b)
+ phistart <- check.theta(phi.start, ns, y, min=0, max=1)
+
- ## call C++ code to draw sample
- posterior <- .C("MCMCbinaryChange",
- phiout = as.double(rep(0.0, nstore*ns)),
- Pout = as.double(rep(0.0, nstore*ns*ns)),
- psout = as.double(rep(0.0, n*ns)),
- sout = as.double(rep(0.0, nstore*n)),
- Ydata = as.double(y),
- Yrow = as.integer(nrow(y)),
- Ycol = as.integer(ncol(y)),
- m = as.integer(m),
- burnin = as.integer(burnin),
- mcmc = as.integer(mcmc),
- thin = as.integer(thin),
- verbose = as.integer(verbose),
- lecuyer=as.integer(lecuyer),
- seedarray=as.integer(seed.array),
- lecuyerstream=as.integer(lecuyer.stream),
- phistart = as.double(phistart),
- Pstart = as.double(Pstart),
- a = as.double(a),
- b = as.double(b),
- c0 = as.double(c0),
- d0 = as.double(d0),
- A0data = as.double(A0),
- logmarglikeholder = as.double(0.0),
- chib = as.integer(chib))
+ ## call C++ code to draw sample
+ posterior <- .C("MCMCbinaryChange",
+ phiout = as.double(rep(0.0, nstore*ns)),
+ Pout = as.double(rep(0.0, nstore*ns*ns)),
+ psout = as.double(rep(0.0, n*ns)),
+ sout = as.double(rep(0.0, nstore*n)),
+ Ydata = as.double(y),
+ Yrow = as.integer(nrow(y)),
+ Ycol = as.integer(ncol(y)),
+ m = as.integer(m),
+ burnin = as.integer(burnin),
+ mcmc = as.integer(mcmc),
+ thin = as.integer(thin),
+ verbose = as.integer(verbose),
+ lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+ phistart = as.double(phistart),
+ Pstart = as.double(Pstart),
+ a = as.double(a),
+ b = as.double(b),
+ c0 = as.double(c0),
+ d0 = as.double(d0),
+ A0data = as.double(A0),
+ logmarglikeholder = as.double(0.0),
+ chib = as.integer(chib))
- ## get marginal likelihood if Chib95
- if (marginal.likelihood == "Chib95"){
- logmarglike <- posterior$logmarglikeholder
- }
-
- ## pull together matrix and build MCMC object to return
- phi.holder <- matrix(posterior$phiout, nstore, )
- P.holder <- matrix(posterior$Pout, nstore, )
- s.holder <- matrix(posterior$sout, nstore, )
- ps.holder <- matrix(posterior$psout, n, )
+ ## get marginal likelihood if Chib95
+ if (marginal.likelihood == "Chib95"){
+ logmarglike <- posterior$logmarglikeholder
+ }
+
+ ## pull together matrix and build MCMC object to return
+ phi.holder <- matrix(posterior$phiout, nstore, )
+ P.holder <- matrix(posterior$Pout, nstore, )
+ s.holder <- matrix(posterior$sout, nstore, )
+ ps.holder <- matrix(posterior$psout, n, )
- output <- mcmc(data=phi.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
- varnames(output) <- paste("phi.", 1:ns, sep = "")
- attr(output,"title") <- "MCMCbinaryChange Posterior Sample"
- attr(output, "y") <- y
- attr(output, "m") <- m
- attr(output, "call") <- cl
- attr(output, "logmarglike") <- logmarglike
- attr(output, "prob.state") <- ps.holder/nstore
- attr(output, "s.store") <- s.holder
+ output <- mcmc(data=phi.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
+ varnames(output) <- paste("phi.", 1:ns, sep = "")
+ attr(output,"title") <- "MCMCbinaryChange Posterior Sample"
+ attr(output, "y") <- y
+ attr(output, "m") <- m
+ attr(output, "call") <- cl
+ attr(output, "logmarglike") <- logmarglike
+ attr(output, "prob.state") <- ps.holder/nstore
+ attr(output, "s.store") <- s.holder
+ }
return(output)
-
}## end of MCMC function
diff --git a/R/MCMCoprobit.R b/R/MCMCoprobit.R
index ef7551b..0af2139 100644
--- a/R/MCMCoprobit.R
+++ b/R/MCMCoprobit.R
@@ -36,7 +36,7 @@
mt <- terms(formula, data=data)
if(missing(data)) data <- sys.frame(sys.parent())
mf <- match.call(expand.dots = FALSE)
- mf$burnin <- mf$mcmc <- mf$b0 <- mf$B0 <- NULL
+ mf$burnin <- mf$mcmc <- mf$b0 <- mf$B0 <- mf$a0 <- mf$A0 <- NULL
mf$thin <- mf$... <- mf$tune <- mf$tdf <- mf$verbose <- mf$seed <- NULL
mf$beta.start <- mf$mcmc.method <- NULL
mf$drop.unused.levels <- TRUE
diff --git a/R/MCMCpoissonChange.R b/R/MCMCpoissonChange.R
index e525757..15e28a0 100644
--- a/R/MCMCpoissonChange.R
+++ b/R/MCMCpoissonChange.R
@@ -36,9 +36,6 @@
lecuyer.stream <- seeds[[3]]
if(!is.na(seed)) set.seed(seed)
- ## prior
- A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
-
if (k==1){
if (is.na(c0)||is.na(d0))
stop("You have to prior for lambda (c0 and d0) when there is no covariate.\n")
@@ -57,86 +54,104 @@
if (marginal.likelihood == "Chib95"){
chib <- 1
}
-
- ## get initial values of tau from observed y
- Pstart <- check.P(P.start, m, a=a, b=b)
- betastart <- beta.change.start(beta.start, ns, k, formula, family=poisson, data)
- if (k == 1){
- betastart <- exp(betastart)
+ if (m == 0){
+ if (marginal.likelihood == "Chib95"){
+ output <- MCMCpoisson(formula, burnin = burnin, mcmc = mcmc,
+ thin = thin, verbose = verbose,
+ b0 = b0, B0 = B0,
+ marginal.likelihood = "Laplace")
+ cat("\n Chib95 method is not yet available for m = 0. Laplace method is used instead.")
+ }
+ else {
+ output <- MCMCpoisson(formula, burnin = burnin, mcmc = mcmc,
+ thin = thin, verbose = verbose,
+ b0 = b0, B0 = B0)
+ }
}
- taustart <- tau.initial(y, tot.comp)
- componentstart <- round(runif(tot.comp, 1, 5))
+ else {
+ ## prior
+ A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
- ## normal mixture weights
- wr <- c(0.2294, 0.2590, 0.2480, 0.1525, 0.0472)
- mr <- c(0.0982, -1.5320, -0.7433, 0.8303, -3.1428)
- sr <- sqrt(c(0.2401, 1.1872, 0.3782, 0.1920, 3.2375))
-
- ## call C++ code to draw sample
- posterior <- .C("MCMCpoissonChange",
- betaout = as.double(rep(0.0, nstore*ns*k)),
- Pout = as.double(rep(0.0, nstore*ns*ns)),
- psout = as.double(rep(0.0, n*ns)),
- sout = as.double(rep(0.0, nstore*n)),
- Ydata = as.double(y),
- Yrow = as.integer(nrow(y)),
- Ycol = as.integer(ncol(y)),
- Xdata = as.double(X),
- Xrow = as.integer(nrow(X)),
- Xcol = as.integer(ncol(X)),
- m = as.integer(m),
- burnin = as.integer(burnin),
- mcmc = as.integer(mcmc),
- thin = as.integer(thin),
- verbose = as.integer(verbose),
- betastart = as.double(betastart),
- Pstart = as.double(Pstart),
- taustart = as.double(taustart),
- componentstart = as.double(componentstart),
- a = as.double(a),
- b = as.double(b),
- c0 = as.double(c0),
- d0 = as.double(d0),
- lecuyer = as.integer(lecuyer),
- seedarray = as.integer(seed.array),
- lecuyerstream = as.integer(lecuyer.stream),
- b0data = as.double(b0),
- B0data = as.double(B0),
- A0data = as.double(A0),
- logmarglikeholder = as.double(0.0),
- loglikeholder = as.double(0.0),
- wrin = as.double(wr),
- mrin = as.double(mr),
- srin = as.double(sr),
- chib = as.integer(chib),
-
- PACKAGE="MCMCpack")
-
- ## get marginal likelihood if Chib95
- if (marginal.likelihood == "Chib95"){
- logmarglike <- posterior$logmarglikeholder
- loglike <- posterior$loglikeholder
- }
-
- ## pull together matrix and build MCMC object to return
- beta.holder <- matrix(posterior$betaout, nstore, )
- P.holder <- matrix(posterior$Pout, nstore, )
- s.holder <- matrix(posterior$sout, nstore, )
- ps.holder <- matrix(posterior$psout, n, )
+ ## get initial values of tau from observed y
+ Pstart <- check.P(P.start, m, a=a, b=b)
+ betastart <- beta.change.start(beta.start, ns, k, formula, family=poisson, data)
+ if (k == 1){
+ betastart <- exp(betastart)
+ }
+ taustart <- tau.initial(y, tot.comp)
+ componentstart <- round(runif(tot.comp, 1, 5))
+
+ ## normal mixture weights
+ wr <- c(0.2294, 0.2590, 0.2480, 0.1525, 0.0472)
+ mr <- c(0.0982, -1.5320, -0.7433, 0.8303, -3.1428)
+ sr <- sqrt(c(0.2401, 1.1872, 0.3782, 0.1920, 3.2375))
- output <- mcmc(data=beta.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
- varnames(output) <- sapply(c(1:ns), function(i) { paste(xnames, "_regime", i, sep = "")})
- attr(output, "title") <- "MCMCpoissonChange Posterior Sample"
- attr(output, "formula") <- formula
- attr(output, "y") <- y
- attr(output, "X") <- X
- attr(output, "m") <- m
- attr(output, "call") <- cl
- attr(output, "logmarglike") <- logmarglike
- attr(output, "loglike") <- loglike
- attr(output, "prob.state") <- ps.holder/nstore
- attr(output, "s.store") <- s.holder
- attr(output, "P.store") <- P.holder
+ ## call C++ code to draw sample
+ posterior <- .C("MCMCpoissonChange",
+ betaout = as.double(rep(0.0, nstore*ns*k)),
+ Pout = as.double(rep(0.0, nstore*ns*ns)),
+ psout = as.double(rep(0.0, n*ns)),
+ sout = as.double(rep(0.0, nstore*n)),
+ Ydata = as.double(y),
+ Yrow = as.integer(nrow(y)),
+ Ycol = as.integer(ncol(y)),
+ Xdata = as.double(X),
+ Xrow = as.integer(nrow(X)),
+ Xcol = as.integer(ncol(X)),
+ m = as.integer(m),
+ burnin = as.integer(burnin),
+ mcmc = as.integer(mcmc),
+ thin = as.integer(thin),
+ verbose = as.integer(verbose),
+ betastart = as.double(betastart),
+ Pstart = as.double(Pstart),
+ taustart = as.double(taustart),
+ componentstart = as.double(componentstart),
+ a = as.double(a),
+ b = as.double(b),
+ c0 = as.double(c0),
+ d0 = as.double(d0),
+ lecuyer = as.integer(lecuyer),
+ seedarray = as.integer(seed.array),
+ lecuyerstream = as.integer(lecuyer.stream),
+ b0data = as.double(b0),
+ B0data = as.double(B0),
+ A0data = as.double(A0),
+ logmarglikeholder = as.double(0.0),
+ loglikeholder = as.double(0.0),
+ wrin = as.double(wr),
+ mrin = as.double(mr),
+ srin = as.double(sr),
+ chib = as.integer(chib),
+
+ PACKAGE="MCMCpack")
+
+ ## get marginal likelihood if Chib95
+ if (marginal.likelihood == "Chib95"){
+ logmarglike <- posterior$logmarglikeholder
+ loglike <- posterior$loglikeholder
+ }
+
+ ## pull together matrix and build MCMC object to return
+ beta.holder <- matrix(posterior$betaout, nstore, )
+ P.holder <- matrix(posterior$Pout, nstore, )
+ s.holder <- matrix(posterior$sout, nstore, )
+ ps.holder <- matrix(posterior$psout, n, )
+
+ output <- mcmc(data=beta.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
+ varnames(output) <- sapply(c(1:ns), function(i) { paste(xnames, "_regime", i, sep = "")})
+ attr(output, "title") <- "MCMCpoissonChange Posterior Sample"
+ attr(output, "formula") <- formula
+ attr(output, "y") <- y
+ attr(output, "X") <- X
+ attr(output, "m") <- m
+ attr(output, "call") <- cl
+ attr(output, "logmarglike") <- logmarglike
+ attr(output, "loglike") <- loglike
+ attr(output, "prob.state") <- ps.holder/nstore
+ attr(output, "s.store") <- s.holder
+ attr(output, "P.store") <- P.holder
+ }
return(output)
-
+
}
diff --git a/R/MCMCprobitChange.R b/R/MCMCprobitChange.R
index 0eb8fdf..72f44f9 100644
--- a/R/MCMCprobitChange.R
+++ b/R/MCMCprobitChange.R
@@ -41,82 +41,96 @@
mvn.prior <- form.mvn.prior(b0, B0, k)
b0 <- mvn.prior[[1]]
B0 <- mvn.prior[[2]]
- A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
-
-
- ## initial values
- Pstart <- check.P(P.start, m, a=a, b=b)
- betastart <- beta.change.start(beta.start, ns, k, formula, family=binomial, data)
-
+
chib <- 0
if (marginal.likelihood == "Chib95"){
chib <- 1
}
-
- ## call C++ code to draw sample
- posterior <- .C("MCMCprobitChange",
- betaout = as.double(rep(0.0, nstore*ns*k)),
- Pout = as.double(rep(0.0, nstore*ns*ns)),
- psout = as.double(rep(0.0, n*ns)),
- sout = as.double(rep(0.0, nstore*n)),
- Ydata = as.double(y),
- Yrow = as.integer(nrow(y)),
- Ycol = as.integer(ncol(y)),
- Xdata = as.double(X),
- Xrow = as.integer(nrow(X)),
- Xcol = as.integer(ncol(X)),
- m = as.integer(m),
- burnin = as.integer(burnin),
- mcmc = as.integer(mcmc),
- thin = as.integer(thin),
- verbose = as.integer(verbose),
-
- lecuyer=as.integer(lecuyer),
- seedarray=as.integer(seed.array),
- lecuyerstream=as.integer(lecuyer.stream),
- betastart = as.double(betastart),
- Pstart = as.double(Pstart),
-
- a = as.double(a),
- b = as.double(b),
- b0data = as.double(b0),
- B0data = as.double(B0),
- A0data = as.double(A0),
- logmarglikeholder = as.double(0.0),
- loglikeholder = as.double(0.0),
- chib = as.integer(chib))
-
- ## get marginal likelihood if Chib95
- if (chib==1){
- logmarglike <- posterior$logmarglikeholder
- loglike <- posterior$loglikeholder
+ if (m == 0){
+ if (marginal.likelihood == "Chib95"){
+ output <- MCMCprobit(formula=Y~X-1, burnin = burnin, mcmc = mcmc,
+ thin = thin, verbose = verbose,
+ b0 = b0, B0 = B0,
+ marginal.likelihood = "Laplace")
+ cat("\n Chib95 method is not yet available for m = 0. Laplace method is used instead.")
+ }
+ else {
+ output <- MCMCprobit(formula=Y~X-1, burnin = burnin, mcmc = mcmc,
+ thin = thin, verbose = verbose,
+ b0 = b0, B0 = B0)
+ }
+ attr(output, "y") <- y
}
else{
- logmarglike <- loglike <- 0
+ A0 <- trans.mat.prior(m=m, n=n, a=a, b=b)
+ Pstart <- check.P(P.start, m, a=a, b=b)
+ betastart <- beta.change.start(beta.start, ns, k, formula, family=binomial, data)
+ ## call C++ code to draw sample
+ posterior <- .C("MCMCprobitChange",
+ betaout = as.double(rep(0.0, nstore*ns*k)),
+ Pout = as.double(rep(0.0, nstore*ns*ns)),
+ psout = as.double(rep(0.0, n*ns)),
+ sout = as.double(rep(0.0, nstore*n)),
+ Ydata = as.double(y),
+ Yrow = as.integer(nrow(y)),
+ Ycol = as.integer(ncol(y)),
+ Xdata = as.double(X),
+ Xrow = as.integer(nrow(X)),
+ Xcol = as.integer(ncol(X)),
+ m = as.integer(m),
+ burnin = as.integer(burnin),
+ mcmc = as.integer(mcmc),
+ thin = as.integer(thin),
+ verbose = as.integer(verbose),
+
+ lecuyer=as.integer(lecuyer),
+ seedarray=as.integer(seed.array),
+ lecuyerstream=as.integer(lecuyer.stream),
+
+ betastart = as.double(betastart),
+ Pstart = as.double(Pstart),
+
+ a = as.double(a),
+ b = as.double(b),
+ b0data = as.double(b0),
+ B0data = as.double(B0),
+ A0data = as.double(A0),
+ logmarglikeholder = as.double(0.0),
+ loglikeholder = as.double(0.0),
+ chib = as.integer(chib))
+
+ ## get marginal likelihood if Chib95
+ if (chib==1){
+ logmarglike <- posterior$logmarglikeholder
+ loglike <- posterior$loglikeholder
+ }
+ else{
+ logmarglike <- loglike <- 0
+ }
+
+ ## pull together matrix and build MCMC object to return
+ beta.holder <- matrix(posterior$betaout, nstore, ns*k)
+ P.holder <- matrix(posterior$Pout, nstore, )
+ s.holder <- matrix(posterior$sout, nstore, )
+ ps.holder <- matrix(posterior$psout, n, )
+
+ output <- mcmc(data=beta.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
+ varnames(output) <- sapply(c(1:ns),
+ function(i){
+ paste(c(xnames), "_regime", i, sep = "")
+ })
+ attr(output, "title") <- "MCMCprobitChange Posterior Sample"
+ attr(output, "formula") <- formula
+ attr(output, "y") <- y
+ attr(output, "X") <- X
+ attr(output, "m") <- m
+ attr(output, "call") <- cl
+ attr(output, "logmarglike") <- logmarglike
+ attr(output, "loglike") <- loglike
+ attr(output, "prob.state") <- ps.holder/nstore
+ attr(output, "s.store") <- s.holder
}
-
- ## pull together matrix and build MCMC object to return
- beta.holder <- matrix(posterior$betaout, nstore, ns*k)
- P.holder <- matrix(posterior$Pout, nstore, )
- s.holder <- matrix(posterior$sout, nstore, )
- ps.holder <- matrix(posterior$psout, n, )
-
- output <- mcmc(data=beta.holder, start=burnin+1, end=burnin + mcmc, thin=thin)
- varnames(output) <- sapply(c(1:ns),
- function(i){
- paste(c(xnames), "_regime", i, sep = "")
- })
- attr(output, "title") <- "MCMCprobitChange Posterior Sample"
- attr(output, "formula") <- formula
- attr(output, "y") <- y
- attr(output, "X") <- X
- attr(output, "m") <- m
- attr(output, "call") <- cl
- attr(output, "logmarglike") <- logmarglike
- attr(output, "loglike") <- loglike
- attr(output, "prob.state") <- ps.holder/nstore
- attr(output, "s.store") <- s.holder
return(output)
}## end of MCMC function
diff --git a/data/PErisk.rda b/data/PErisk.rda
index d60ae32..7688d70 100644
Binary files a/data/PErisk.rda and b/data/PErisk.rda differ
diff --git a/data/SupremeCourt.rda b/data/SupremeCourt.rda
index 3b07c6c..1f1d1e7 100644
Binary files a/data/SupremeCourt.rda and b/data/SupremeCourt.rda differ
diff --git a/man/MCMCbinaryChange.Rd b/man/MCMCbinaryChange.Rd
index 24d5dfd..e9fcbe0 100755
--- a/man/MCMCbinaryChange.Rd
+++ b/man/MCMCbinaryChange.Rd
@@ -116,13 +116,14 @@
y3 <- rbinom(120, 1, true.phi[3])
y <- as.ts(c(y1, y2, y3))
+ model0 <- MCMCbinaryChange(y, m=0, c0=2, d0=2, mcmc=1000, burnin=1000, verbose=500, marginal.likelihood = "Chib95")
model1 <- MCMCbinaryChange(y, m=1, c0=2, d0=2, mcmc=1000, burnin=1000, verbose=500, marginal.likelihood = "Chib95")
model2 <- MCMCbinaryChange(y, m=2, c0=2, d0=2, mcmc=1000, burnin=1000, verbose=500, marginal.likelihood = "Chib95")
model3 <- MCMCbinaryChange(y, m=3, c0=2, d0=2, mcmc=1000, burnin=1000, verbose=500, marginal.likelihood = "Chib95")
model4 <- MCMCbinaryChange(y, m=4, c0=2, d0=2, mcmc=1000, burnin=1000, verbose=500, marginal.likelihood = "Chib95")
model5 <- MCMCbinaryChange(y, m=5, c0=2, d0=2, mcmc=1000, burnin=1000, verbose=500, marginal.likelihood = "Chib95")
- print(BayesFactor(model1, model2, model3, model4, model5))
+ print(BayesFactor(model0, model1, model2, model3, model4, model5))
## plot two plots in one screen
par(mfrow=c(attr(model2, "m") + 1, 1), mai=c(0.4, 0.6, 0.3, 0.05))
diff --git a/man/MCMCoprobit.Rd b/man/MCMCoprobit.Rd
index 6913fee..5859650 100644
--- a/man/MCMCoprobit.Rd
+++ b/man/MCMCoprobit.Rd
@@ -72,17 +72,17 @@ MCMCoprobit(formula, data = parent.frame(), burnin = 1000, mcmc = 10000,
Default value of 0 is equivalent to an improper uniform prior on
\eqn{\beta}{beta}. }
- \item{a0}{The prior mean of \eqn{\alpha}{alpha}. This can either be a
+ \item{a0}{The prior mean of \eqn{\gamma}{gamma}. This can either be a
scalar or a column vector with dimension equal to the number of
betas. If this takes a scalar value, then that value will serve as
the prior mean for all of the betas.}
- \item{A0}{The prior precision of \eqn{\alpha}{alpha}. This can either be a
+ \item{A0}{The prior precision of \eqn{\gamma}{gamma}. This can either be a
scalar or a square matrix with dimensions equal to the number of
betas. If this takes a scalar value, then that value times an
- identity matrix serves as the prior precision of \eqn{\beta}{beta}.
+ identity matrix serves as the prior precision of \eqn{\gamma}{gamma}.
Default value of 0 is equivalent to an improper uniform prior on
- \eqn{\beta}{beta}. }
+ \eqn{\gamma}{gamma}. }
\item{mcmc.method}{Can be set to either "Cowles" (default) or "AC" to perform
posterior sampling of cutpoints based on Cowles (1996) or Albert and Chib (2001) respectively.}
diff --git a/man/MCMCpoissonChange.Rd b/man/MCMCpoissonChange.Rd
index 9990a1a..ff94937 100644
--- a/man/MCMCpoissonChange.Rd
+++ b/man/MCMCpoissonChange.Rd
@@ -143,6 +143,9 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
formula = y ~ x1
## fit multiple models with a varying number of breaks
+ model0 <- MCMCpoissonChange(formula, m=0,
+ mcmc = 1000, burnin = 1000, verbose = 500,
+ b0 = rep(0, 2), B0 = 5*diag(2), marginal.likelihood = "Chib95")
model1 <- MCMCpoissonChange(formula, m=1,
mcmc = 1000, burnin = 1000, verbose = 500,
b0 = rep(0, 2), B0 = 5*diag(2), marginal.likelihood = "Chib95")
@@ -160,7 +163,7 @@ Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park. 2011.
b0 = rep(0, 2), B0 = 5*diag(2), marginal.likelihood = "Chib95")
## find the most reasonable one
- print(BayesFactor(model1, model2, model3, model4, model5))
+ print(BayesFactor(model0, model1, model2, model3, model4, model5))
## draw plots using the "right" model
par(mfrow=c(attr(model2, "m") + 1, 1), mai=c(0.4, 0.6, 0.3, 0.05))
diff --git a/man/MCMCprobitChange.Rd b/man/MCMCprobitChange.Rd
index 9702c96..c0963eb 100644
--- a/man/MCMCprobitChange.Rd
+++ b/man/MCMCprobitChange.Rd
@@ -146,9 +146,9 @@ y3 <- rbinom(100, 1, true.phi3)
Y <- as.ts(c(y1, y2, y3))
## fit multiple models with a varying number of breaks
-out0 <- MCMCprobit(formula=Y~X-1, data=parent.frame(),
- mcmc=1000, burnin=1000, thin=1, verbose=1000,
- b0 = 0, B0 = 10, marginal.likelihood = c("Laplace"))
+out0 <- MCMCprobitChange(formula=Y~X-1, data=parent.frame(), m=0,
+ mcmc=1000, burnin=1000, thin=1, verbose=1000,
+ b0 = 0, B0 = 10, a = 1, b = 1, marginal.likelihood = c("Chib95"))
out1 <- MCMCprobitChange(formula=Y~X-1, data=parent.frame(), m=1,
mcmc=1000, burnin=1000, thin=1, verbose=1000,
b0 = 0, B0 = 10, a = 1, b = 1, marginal.likelihood = c("Chib95"))
diff --git a/man/plotChangepoint.Rd b/man/plotChangepoint.Rd
old mode 100644
new mode 100755
diff --git a/man/plotState.Rd b/man/plotState.Rd
old mode 100644
new mode 100755
diff --git a/src/MCMCbinaryChange.cc b/src/MCMCbinaryChange.cc
index 87d8a55..0689688 100644
--- a/src/MCMCbinaryChange.cc
+++ b/src/MCMCbinaryChange.cc
@@ -99,7 +99,7 @@ Matrix<> binary_state_sampler(rng<RNGTYPE>& stream,
template <typename RNGTYPE>
void MCMCbinaryChange_impl(rng<RNGTYPE>& stream, const Matrix<>& Y,
Matrix<>& phi, Matrix<>& P, const Matrix<>& A0,
- double m, double c0, double d0,
+ const double m, const double c0, const double d0,
unsigned int burnin, unsigned int mcmc, unsigned int thin,
unsigned int verbose, bool chib,
Matrix<>& phi_store, Matrix<>& P_store,
@@ -140,9 +140,9 @@ void MCMCbinaryChange_impl(rng<RNGTYPE>& stream, const Matrix<>& Y,
}
double c1 = addY[j] + c0;
double d1 = addN[j] - addY[j] + d0;
- phi[j] = stream.rbeta(c1, d1);
-
+ phi[j] = stream.rbeta(c1, d1);
}
+
//////////////////////
// 3. Sample P
//////////////////////
@@ -169,8 +169,6 @@ void MCMCbinaryChange_impl(rng<RNGTYPE>& stream, const Matrix<>& Y,
++count;
}
-
-
if(verbose > 0 && iter % verbose == 0){
Rprintf("\n\n MCMCbinaryChange iteration %i of %i", (iter+1), tot_iter);
for (int j = 0;j<ns; ++j){
@@ -197,6 +195,7 @@ void MCMCbinaryChange_impl(rng<RNGTYPE>& stream, const Matrix<>& Y,
// phi
//////////////////////
Matrix<> density_phi(nstore, ns);
+ // Matrix<> density_phi_j(ns, 1);
for (int iter = 0; iter<nstore; ++iter){
Matrix<> addY(ns, 1);
@@ -204,19 +203,18 @@ void MCMCbinaryChange_impl(rng<RNGTYPE>& stream, const Matrix<>& Y,
for (int j = 0; j<ns; ++j){
for (int i = 0; i<n; ++i){
- if (s_store(iter, i) == (j+1)) {
- addN[j] = addN[j] + 1;
- addY[j] = addY[j] + Y[i];
- }
+ if (s_store(iter, i) == (j+1)) {
+ addN[j] = addN[j] + 1;
+ addY[j] = addY[j] + Y[i];
+ }
}
double c1 = addY[j] + c0;
double d1 = addN[j] - addY[j] + d0;
- density_phi(iter, j) = dbeta(phi_st[j], c1, d1);
+ density_phi(iter, j) = dbeta(phi_st[j], c1, d1);
}
-
}
double pdf_phi = log(prod(meanc(density_phi)));
-
+
//////////////////////
// P
//////////////////////
@@ -289,16 +287,24 @@ void MCMCbinaryChange_impl(rng<RNGTYPE>& stream, const Matrix<>& Y,
Matrix<> density_P_prior(ns, 1);
for (int j=0; j<ns ; ++j){
- density_phi_prior[j] = log(dbeta(phi_st[j], c0, d0));
+ density_phi_prior[j] = scythe::lndbeta1(phi_st[j], c0, d0);
}
for (int j =0; j< (ns-1); ++j){
- density_P_prior[j] = log(dbeta(P_st(j,j), A0(j,j), A0(j,j+1)));
+ density_P_prior[j] = scythe::lndbeta1(P_st(j,j), A0(j,j), A0(j,j+1));
}
// compute marginal likelihood
double logprior = sum(density_phi_prior) + sum(density_P_prior);
logmarglike = (loglike + logprior) - (pdf_phi + pdf_P);
+ if(verbose > 0){
+ Rprintf("\nlogmarglike = %10.5f\n", logmarglike);
+ Rprintf("loglike = %10.5f\n", loglike);
+ Rprintf("log_prior = %10.5f\n", logprior);
+ Rprintf("log_phi = %10.5f\n", pdf_phi);
+ Rprintf("log_P = %10.5f\n", pdf_P);
+ }
+
}
}
////////////////////////////////////////////
diff --git a/src/MCMChlogit.cc b/src/MCMChlogit.cc
index 49f2391..250f0a3 100644
--- a/src/MCMChlogit.cc
+++ b/src/MCMChlogit.cc
@@ -442,7 +442,7 @@ extern "C"{
//////////////////////////////////////////////////
// User interrupt
- void R_CheckUserInterrupt(void); // allow user interrupts
+ R_CheckUserInterrupt(); // allow user interrupts
} // end MCMC loop
diff --git a/src/MCMChpoisson.cc b/src/MCMChpoisson.cc
index fb4ad62..57eac31 100644
--- a/src/MCMChpoisson.cc
+++ b/src/MCMChpoisson.cc
@@ -430,7 +430,7 @@ extern "C"{
//////////////////////////////////////////////////
// User interrupt
- void R_CheckUserInterrupt(void); // allow user interrupts
+ R_CheckUserInterrupt(); // allow user interrupts
} // end MCMC loop
diff --git a/src/MCMChregress.cc b/src/MCMChregress.cc
index dfe42f9..b784c14 100644
--- a/src/MCMChregress.cc
+++ b/src/MCMChregress.cc
@@ -219,8 +219,8 @@ extern "C"{
// invVi, sum_V and sum_v
// see http://en.wikipedia.org/wiki/Woodbury_matrix_identity with A=(1/V_run)*Id
- Matrix<double> sum_V(NQ,NQ);
- Matrix<double> sum_v(NQ,1);
+ Matrix<double> sum_V(NP,NP);
+ Matrix<double> sum_v(NP,1);
for (int k=0; k<NGROUP; k++) {
sum_V += (1/V_run)*cpXk_arr[k]-pow(1/V_run,2)*tXWk_arr[k]*invpd(invpd(Vb_run)+tWk_arr[k]*(1/V_run)*Wk_arr[k])*tWXk_arr[k];
sum_v += (1/V_run)*tXYk_arr[k]-pow(1/V_run,2)*tXWk_arr[k]*invpd(invpd(Vb_run)+tWk_arr[k]*(1/V_run)*Wk_arr[k])*tWYk_arr[k];
@@ -349,7 +349,7 @@ extern "C"{
//////////////////////////////////////////////////
// User interrupt
- void R_CheckUserInterrupt(void); // allow user interrupts
+ R_CheckUserInterrupt(); // allow user interrupts
} // end MCMC loop
diff --git a/src/MCMClogit.cc b/src/MCMClogit.cc
index 16396cb..74453d8 100644
--- a/src/MCMClogit.cc
+++ b/src/MCMClogit.cc
@@ -119,11 +119,12 @@ MCMClogit_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
R_CheckUserInterrupt(); // allow user interrupts
}// end MCMC loop
-
- Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
- Rprintf("The Metropolis acceptance rate for beta was %3.5f",
- static_cast<double>(accepts) / static_cast<double>(tot_iter));
- Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+ if (verbose > 0){
+ Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+ Rprintf("The Metropolis acceptance rate for beta was %3.5f",
+ static_cast<double>(accepts) / static_cast<double>(tot_iter));
+ Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+ }
}
extern "C"{
diff --git a/src/MCMCoprobit.cc b/src/MCMCoprobit.cc
index aaef110..8fa9b7c 100644
--- a/src/MCMCoprobit.cc
+++ b/src/MCMCoprobit.cc
@@ -350,12 +350,12 @@ void MCMCoprobit_impl (rng<RNGTYPE>& stream, const int * Y,
R_CheckUserInterrupt(); // allow user interrupts
}// end of MCMC
-
- Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
- Rprintf("The Metropolis acceptance rate for beta was %3.5f",
- static_cast<double>(accepts) / static_cast<double>(tot_iter));
- Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
-
+ if (verbose > 0){
+ Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+ Rprintf("The Metropolis acceptance rate for beta was %3.5f",
+ static_cast<double>(accepts) / static_cast<double>(tot_iter));
+ Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+ }
result = storemat;
}
@@ -390,7 +390,8 @@ extern "C"{
const Matrix <> b0(*b0row, *b0col, b0data);
const Matrix <> B0(*B0row, *B0col, B0data);
const Matrix <> alpha_prior_mean(*a0row, *a0col, a0data);
- const Matrix <> alpha_prior_var(*A0row, *A0col, A0data);
+ const Matrix <> alpha_prior_prec(*A0row, *A0col, A0data);
+ const Matrix <> alpha_prior_var = invpd(alpha_prior_prec);
const Matrix<> tune(*tunerow, *tunecol, tunedata);
Matrix<> storagematrix;
diff --git a/src/MCMCoprobitChange.cc b/src/MCMCoprobitChange.cc
index 14ec4cb..7784fe4 100644
--- a/src/MCMCoprobitChange.cc
+++ b/src/MCMCoprobitChange.cc
@@ -270,11 +270,11 @@ void MCMCoprobitChange_impl(rng<RNGTYPE>& stream,
Matrix<> XpX = t(Xj)*Xj;
Matrix<> XpZ = t(Xj)*Zj;
Matrix<> XpY = t(Xj)*yj;
- Matrix<> Bn = invpd(B0inv + XpX/Sigma[0]);
- Matrix<> bn = Bn*(B0inv*b0 + XpY/Sigma[0]);
+ Matrix<> Bn = invpd(B0 + XpX/Sigma[0]);
+ Matrix<> bn = Bn*(B0*b0 + XpY/Sigma[0]);
beta_linear(j,_) = stream.rmvnorm(bn, Bn);
- Matrix<> Bn2 = invpd(B0inv + XpX);
- Matrix<> bn2 = Bn2*(B0inv*b0 + XpZ);
+ Matrix<> Bn2 = invpd(B0 + XpX);
+ Matrix<> bn2 = Bn2*(B0*b0 + XpZ);
beta(j,_) = stream.rmvnorm(bn2, Bn2);
beta_count_storage[j] = beta_count;
}
@@ -476,11 +476,11 @@ void MCMCoprobitChange_impl(rng<RNGTYPE>& stream,
Matrix<> XpX = t(Xj)*Xj;
Matrix<> XpZ = t(Xj)*Zj;
Matrix<> XpY = t(Xj)*yj;
- Matrix<> Bn = invpd(B0inv + XpX/Sigma[0]);
- Matrix<> bn = Bn*(B0inv*b0 + XpY/Sigma[0]);
+ Matrix<> Bn = invpd(B0 + XpX/Sigma[0]);
+ Matrix<> bn = Bn*(B0*b0 + XpY/Sigma[0]);
beta_linear(j,_) = stream.rmvnorm(bn, Bn);
- Matrix<> Bn2 = invpd(B0inv + XpX);
- Matrix<> bn2 = Bn2*(B0inv*b0 + XpZ);
+ Matrix<> Bn2 = invpd(B0 + XpX);
+ Matrix<> bn2 = Bn2*(B0*b0 + XpZ);
beta(j,_) = stream.rmvnorm(bn2, Bn2);
beta_count_storage[j] = beta_count;
}
@@ -547,11 +547,11 @@ void MCMCoprobitChange_impl(rng<RNGTYPE>& stream,
Matrix<> XpX = t(Xj)*Xj;
Matrix<> XpZ = t(Xj)*Zj;
Matrix<> XpY = t(Xj)*yj;
- Matrix<> Bn = invpd(B0inv + XpX/Sigma[0]);
- Matrix<> bn = Bn*(B0inv*b0 + XpY/Sigma[0]);
+ Matrix<> Bn = invpd(B0 + XpX/Sigma[0]);
+ Matrix<> bn = Bn*(B0*b0 + XpY/Sigma[0]);
beta_linear(j,_) = stream.rmvnorm(bn, Bn);
- Matrix<> Bn2 = invpd(B0inv + XpX);
- Matrix<> bn2 = Bn2*(B0inv*b0 + XpZ);
+ Matrix<> Bn2 = invpd(B0 + XpX);
+ Matrix<> bn2 = Bn2*(B0*b0 + XpZ);
beta(j,_) = stream.rmvnorm(bn2, Bn2);
beta_count_storage[j] = beta_count;
density_beta(iter, j) = exp(lndmvn(t(beta_st(j,_)), bn2, Bn2));
@@ -635,7 +635,7 @@ void MCMCoprobitChange_impl(rng<RNGTYPE>& stream,
density_P[ns-1] = 1; //
for (int j=0; j<ns ; ++j){
- density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0);
+ density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0inv);
}
for (int j =0; j< (ns-1); ++j){
@@ -646,18 +646,14 @@ void MCMCoprobitChange_impl(rng<RNGTYPE>& stream,
double logprior = sum(density_beta_prior) + sum(density_P_prior) + density_gamma_prior;
logmarglike = (loglike + logprior) - (pdf_beta + pdf_P + pdf_gamma);
- /*
- Rprintf("\n ----------------- marginal likelihood outputs ----------------- \n");
- Rprintf("\n logmarglike %10.5f", logmarglike, "\n");
- Rprintf("\n loglike %10.5f", loglike, "\n");
- Rprintf("\n pdf_beta %10.5f", pdf_beta, "\n");
- Rprintf("\n pdf_gamma %10.5f", pdf_gamma, "\n");
- Rprintf("\n pdf_P %10.5f", pdf_P, "\n");
- Rprintf("\n logprior %10.5f", logprior, "\n");
- Rprintf("\n density_beta_prior %10.5f", sum(density_beta_prior), "\n");
- Rprintf("\n density_gamma_prior %10.5f", density_gamma_prior, "\n");
- Rprintf("\n density_P_prior %10.5f", sum(density_P_prior), "\n");
- */
+ if (verbose >0 ){
+ Rprintf("\nlogmarglike = %10.5f\n", logmarglike);
+ Rprintf("loglike = %10.5f\n", loglike);
+ Rprintf("log_prior = %10.5f\n", logprior);
+ Rprintf("log_beta = %10.5f\n", pdf_beta);
+ Rprintf("log_P = %10.5f\n", pdf_P);
+ Rprintf("log_gamma = %10.5f\n", pdf_gamma);
+ }
} // end of marginal likelihood
}//end
diff --git a/src/MCMCpoisson.cc b/src/MCMCpoisson.cc
index 8b3851b..c5478d2 100644
--- a/src/MCMCpoisson.cc
+++ b/src/MCMCpoisson.cc
@@ -74,69 +74,70 @@ static double poisson_logpost(const Matrix<>& Y,
*/
template <typename RNGTYPE>
void MCMCpoisson_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
- const Matrix<>& X, const Matrix<>& tune, Matrix<>& beta, const Matrix<>& b0,
- const Matrix<>& B0, const Matrix<>& V, unsigned int burnin,
- unsigned int mcmc, unsigned int thin, unsigned int verbose, Matrix<>& result) {
-
- // define constants
- const unsigned int tot_iter = burnin + mcmc; // total number iterations
- const unsigned int nstore = mcmc / thin; // number of draws to store
- const unsigned int k = X.cols();
+ const Matrix<>& X, const Matrix<>& tune, Matrix<>& beta, const Matrix<>& b0,
+ const Matrix<>& B0, const Matrix<>& V, unsigned int burnin,
+ unsigned int mcmc, unsigned int thin, unsigned int verbose, Matrix<>& result) {
+
+ // define constants
+ const unsigned int tot_iter = burnin + mcmc; // total number iterations
+ const unsigned int nstore = mcmc / thin; // number of draws to store
+ const unsigned int k = X.cols();
+
+ // storage matrix or matrices
+ Matrix<> storemat(nstore, k);
+
+ // proposal parameters
+ const Matrix<> propV = tune * invpd(B0 + invpd(V)) * tune;
+ const Matrix<> propC = cholesky(propV) ;
- // storage matrix or matrices
- Matrix<> storemat(nstore, k);
+ double logpost_cur = poisson_logpost(Y, X, beta, b0, B0);
+
+ // MCMC loop
+ int count = 0;
+ int accepts = 0;
+ for (unsigned int iter = 0; iter < tot_iter; ++iter){
- // proposal parameters
- const Matrix<> propV = tune * invpd(B0 + invpd(V)) * tune;
- const Matrix<> propC = cholesky(propV) ;
+ // sample beta
+ const Matrix<> beta_can = gaxpy(propC, stream.rnorm(k,1,0,1), beta);
+ const double logpost_can = poisson_logpost(Y,X,beta_can, b0, B0);
+ const double ratio = ::exp(logpost_can - logpost_cur);
- double logpost_cur = poisson_logpost(Y, X, beta, b0, B0);
+ if (stream.runif() < ratio){
+ beta = beta_can;
+ logpost_cur = logpost_can;
+ ++accepts;
+ }
- // MCMC loop
- int count = 0;
- int accepts = 0;
- for (unsigned int iter = 0; iter < tot_iter; ++iter){
-
- // sample beta
- const Matrix<> beta_can = gaxpy(propC, stream.rnorm(k,1,0,1), beta);
- const double logpost_can = poisson_logpost(Y,X,beta_can, b0, B0);
- const double ratio = ::exp(logpost_can - logpost_cur);
-
- if (stream.runif() < ratio){
- beta = beta_can;
- logpost_cur = logpost_can;
- ++accepts;
- }
-
- // store values in matrices
- if (iter >= burnin && (iter % thin==0)){
- storemat(count,_) = beta;
- ++count;
- }
-
- // print output to stdout
- if(verbose > 0 && iter % verbose == 0){
- Rprintf("\n\nMCMCpoisson iteration %i of %i \n", (iter+1), tot_iter);
- Rprintf("beta = \n");
- for (unsigned int j=0; j<k; ++j)
- Rprintf("%10.5f\n", beta[j]);
- Rprintf("Metropolis acceptance rate for beta = %3.5f\n\n",
- static_cast<double>(accepts) /
- static_cast<double>(iter+1));
- }
-
- R_CheckUserInterrupt(); // allow user interrupts
- }// end MCMC loop
- result = storemat;
-
- Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
- Rprintf("The Metropolis acceptance rate for beta was %3.5f",
- static_cast<double>(accepts) / static_cast<double>(tot_iter));
- Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+ // store values in matrices
+ if (iter >= burnin && (iter % thin==0)){
+ storemat(count,_) = beta;
+ ++count;
+ }
+
+ // print output to stdout
+ if(verbose > 0 && iter % verbose == 0){
+ Rprintf("\n\nMCMCpoisson iteration %i of %i \n", (iter+1), tot_iter);
+ Rprintf("beta = \n");
+ for (unsigned int j=0; j<k; ++j)
+ Rprintf("%10.5f\n", beta[j]);
+ Rprintf("Metropolis acceptance rate for beta = %3.5f\n\n",
+ static_cast<double>(accepts) /
+ static_cast<double>(iter+1));
+ }
+
+ R_CheckUserInterrupt(); // allow user interrupts
+ }// end MCMC loop
+
+ result = storemat;
+ if (verbose > 0){
+ Rprintf("\n\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+ Rprintf("The Metropolis acceptance rate for beta was %3.5f",
+ static_cast<double>(accepts) / static_cast<double>(tot_iter));
+ Rprintf("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n");
+ }
}
-
extern "C"{
-void MCMCpoisson(double *sampledata, const int *samplerow,
+ void MCMCpoisson(double *sampledata, const int *samplerow,
const int *samplecol, const double *Ydata,
const int *Yrow, const int *Ycol, const double *Xdata,
const int *Xrow, const int *Xcol, const int *burnin,
@@ -150,23 +151,23 @@ void MCMCpoisson(double *sampledata, const int *samplerow,
const double *Vdata, const int *Vrow, const int *Vcol) {
// pull together Matrix objects
- const Matrix <> Y(*Yrow, *Ycol, Ydata);
- const Matrix <> X(*Xrow, *Xcol, Xdata);
- const Matrix <> tune(*tunerow, *tunecol, tunedata);
- Matrix <> beta(*betastartrow, *betastartcol,
+ const Matrix <> Y(*Yrow, *Ycol, Ydata);
+ const Matrix <> X(*Xrow, *Xcol, Xdata);
+ const Matrix <> tune(*tunerow, *tunecol, tunedata);
+ Matrix <> beta(*betastartrow, *betastartcol,
betastartdata);
- const Matrix <> b0(*b0row, *b0col, b0data);
- const Matrix <> B0(*B0row, *B0col, B0data);
- const Matrix <> V(*Vrow, *Vcol, Vdata);
+ const Matrix <> b0(*b0row, *b0col, b0data);
+ const Matrix <> B0(*B0row, *B0col, B0data);
+ const Matrix <> V(*Vrow, *Vcol, Vdata);
+
+ Matrix<> storagematrix;
+ MCMCPACK_PASSRNG2MODEL(MCMCpoisson_impl, Y, X, tune, beta, b0, B0,
+ V, *burnin, *mcmc, *thin, *verbose, storagematrix);
- Matrix<> storagematrix;
- MCMCPACK_PASSRNG2MODEL(MCMCpoisson_impl, Y, X, tune, beta, b0, B0,
- V, *burnin, *mcmc, *thin, *verbose, storagematrix);
-
const unsigned int size = *samplerow * *samplecol;
for (unsigned int i=0; i<size; ++i)
- sampledata[i] = storagematrix(i);
- }
+ sampledata[i] = storagematrix(i);
+ }
}
#endif
diff --git a/src/MCMCpoissonChange.cc b/src/MCMCpoissonChange.cc
index e9af4c5..791d891 100644
--- a/src/MCMCpoissonChange.cc
+++ b/src/MCMCpoissonChange.cc
@@ -646,8 +646,8 @@ void MCMCpoissonRegChangepoint_impl(rng<RNGTYPE>& stream,
for (int h = 0; h<nstate[j]; ++h){
Xwj(h, _) = Xj(h,_)*wi[h];
}
- Matrix<> Bn = invpd(B0inv + ::t(Xj)*Xwj);
- Matrix<> bn = Bn*gaxpy(B0inv, b0, -1*::t(Xj)*yj);
+ Matrix<> Bn = invpd(B0 + ::t(Xj)*Xwj);
+ Matrix<> bn = Bn*gaxpy(B0, b0, -1*::t(Xj)*yj);
beta(j,_) = stream.rmvnorm(bn, Bn);
}
@@ -760,8 +760,8 @@ void MCMCpoissonRegChangepoint_impl(rng<RNGTYPE>& stream,
for (int h = 0; h<nstate[j]; ++h){
Xwj(h, _) = Xj(h,_)*wi[h];
}
- Matrix<> Bn = invpd(B0inv + ::t(Xj)*Xwj);
- Matrix<> bn = Bn*gaxpy(B0inv, b0, -1*::t(Xj)*yj);
+ Matrix<> Bn = invpd(B0 + ::t(Xj)*Xwj);
+ Matrix<> bn = Bn*gaxpy(B0, b0, -1*::t(Xj)*yj);
density_beta(iter, j) = exp(lndmvn(::t(beta_st(j,_)), bn, Bn));
}
}
@@ -841,7 +841,7 @@ void MCMCpoissonRegChangepoint_impl(rng<RNGTYPE>& stream,
density_P[ns-1] = 1; //
for (int j=0; j<ns ; ++j){
- density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0);
+ density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0inv);
}
for (int j =0; j< (ns-1); ++j){
@@ -851,7 +851,13 @@ void MCMCpoissonRegChangepoint_impl(rng<RNGTYPE>& stream,
// compute marginal likelihood
const double logprior = sum(density_beta_prior) + sum(density_P_prior);
const double logmarglike = (loglike + logprior) - (pdf_beta + pdf_P);
-
+ if (verbose >0 ){
+ Rprintf("\nlogmarglike = %10.5f\n", logmarglike);
+ Rprintf("loglike = %10.5f\n", loglike);
+ Rprintf("log_prior = %10.5f\n", logprior);
+ Rprintf("log_beta = %10.5f\n", pdf_beta);
+ Rprintf("log_P = %10.5f\n", pdf_P);
+ }
logmarglikeholder[0] = logmarglike;
loglikeholder[0] = loglike;
diff --git a/src/MCMCprobit.cc b/src/MCMCprobit.cc
index fb9b2cf..71cf593 100644
--- a/src/MCMCprobit.cc
+++ b/src/MCMCprobit.cc
@@ -66,8 +66,8 @@ void MCMCprobit_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
// storage matrix or matrices
Matrix<> beta_store(nstore, k);
- Matrix<> Z_store(nstore, N);
-
+ Matrix<> bn_store(nstore, k);
+
// initialize Z
Matrix<> Z(N,1);
@@ -86,12 +86,14 @@ void MCMCprobit_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
// [beta|Z]
const Matrix<> XpZ = t(X) * Z;
- beta = NormNormregress_beta_draw(XpX, XpZ, b0, B0, 1.0, stream);
+ const Matrix<double> Bn = invpd(B0 + XpX);
+ const Matrix<double> bn = Bn*gaxpy(B0, b0, XpZ);
+ beta = stream.rmvnorm(bn, Bn);
// store values in matrices
if (iter >= burnin && (iter % thin==0)){
beta_store(count,_) = beta;
- Z_store(count,_) = Z;
+ bn_store(count,_) = bn;
++count;
}
@@ -109,18 +111,18 @@ void MCMCprobit_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
if(chib==1){
// Rprintf("\n Marginal Likelihood Computation Starts!\n");
- Matrix<double> beta_star = meanc(beta_store);
- const Matrix<double> Z_reduced(N, 1);
- double logbeta_sum = 0;
+ Matrix<double> beta_star(k, 1);
+ beta_star(_ ,0) = meanc(beta_store);
+ Matrix<double> bn_reduced(k, 1);
+ Matrix<double> density_beta(nstore, 1);
for (int iter = 0; iter<nstore; ++iter){
- Z_reduced(_,0) = Z_store(iter,_);
- const Matrix<double> XpZ = (::t(X)*Z_reduced);
- const Matrix<double> Bn = invpd(B0inv + XpX);
- const Matrix<double> bn = Bn*gaxpy(B0inv, b0, XpZ);
- logbeta_sum += lndmvn(beta_star, bn, Bn);
+ bn_reduced(_, 0) = bn_store(iter, _);
+ Matrix<double> bn_reduced1 = bn_store(iter, _);
+ const Matrix<double> Bn = invpd(B0 + XpX);
+ density_beta(iter) = ::exp(lndmvn(beta_star, bn_reduced, Bn));
}
- double logbeta = logbeta_sum/nstore;
-
+ double logbeta = log(mean(density_beta));
+
double loglike = 0.0;
Matrix<> eta = X * beta_star;
for (unsigned int i = 0; i < N; ++i) {
@@ -129,14 +131,23 @@ void MCMCprobit_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
}
// calculate log prior ordinate
- double logprior = lndmvn(beta_star, b0, B0inv);
+ double logprior = 0.0;
+ if (k == 1){
+ logprior = log(dnorm(beta_star(0), b0(0), sqrt(B0inv(0))));
+ }
+ else{
+ logprior = lndmvn(beta_star, b0, B0inv);
+ }
+ //
logmarglike = loglike + logprior - logbeta;
-
- // Rprintf("\n logmarglike %10.5f", logmarglike, "\n");
- // Rprintf("\n loglike %10.5f", loglike, "\n");
-
- }// end of marginal likelihood computation
+ if (verbose > 0){
+ Rprintf("\nlogmarglike = %10.5f\n", logmarglike);
+ Rprintf("loglike = %10.5f\n", loglike);
+ Rprintf("log_prior = %10.5f\n", logprior);
+ Rprintf("log_beta = %10.5f\n", logbeta);
+ }
+ }// end of marginal likelihood computation
result = beta_store;
}
diff --git a/src/MCMCprobitChange.cc b/src/MCMCprobitChange.cc
index ced6161..7cbf7f1 100644
--- a/src/MCMCprobitChange.cc
+++ b/src/MCMCprobitChange.cc
@@ -248,8 +248,8 @@ void MCMCprobitChange_impl(rng<RNGTYPE>& stream,
const Matrix<double> Xj = X((beta_count - nstate[j]), 0, (beta_count - 1), k-1);
const Matrix<double> XpX = (::t(Xj)*Xj);
const Matrix<double> XpZ = (::t(Xj)*Zj);
- const Matrix<double> Bn = invpd(B0inv + XpX);
- const Matrix<double> bn = Bn*gaxpy(B0inv, b0, XpZ);
+ const Matrix<double> Bn = invpd(B0 + XpX);
+ const Matrix<double> bn = Bn*gaxpy(B0, b0, XpZ);
density_beta(iter, j) = exp(lndmvn(::t(beta_st(j,_)), bn, Bn));
}
}
@@ -325,24 +325,24 @@ void MCMCprobitChange_impl(rng<RNGTYPE>& stream,
density_P[ns-1] = 1; //
for (int j=0; j<ns ; ++j){
- density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0);
+ density_beta_prior[j] = lndmvn(::t(beta_st(j,_)), b0, B0inv);
}
for (int j =0; j< (ns-1); ++j){
- density_P_prior[j] = log(dbeta(P_st(j,j), A0(j,j), A0(j,j+1)));
+ density_P_prior[j] = scythe::lndbeta1(P_st(j,j), A0(j,j), A0(j,j+1));
}
// compute marginal likelihood
double logprior = sum(density_beta_prior) + sum(density_P_prior);
logmarglike = (loglike + logprior) - (pdf_beta + pdf_P);
- // Rprintf("\n density_beta_prior %10.5f", density_beta_prior[0], "\n");
- // Rprintf("\n density_P_prior %10.5f", density_P_prior[0], "\n");
- Rprintf("\n logmarglike %10.5f", logmarglike, "\n");
- Rprintf("\n loglike %10.5f", loglike, "\n");
- // Rprintf("\n logprior %10.5f", logprior, "\n");
- // Rprintf("\n pdf_beta %10.5f", pdf_beta, "\n");
- // Rprintf("\n pdf_P %10.5f", pdf_P, "\n");
+ if (verbose >0 ){
+ Rprintf("\nlogmarglike = %10.5f\n", logmarglike);
+ Rprintf("loglike = %10.5f\n", loglike);
+ Rprintf("log_prior = %10.5f\n", logprior);
+ Rprintf("log_beta = %10.5f\n", pdf_beta);
+ Rprintf("log_P = %10.5f\n", pdf_P);
+ }
} // end of marginal likelihood
}//end extern "C"
diff --git a/src/MCMCregress.cc b/src/MCMCregress.cc
index 32eb701..bcad188 100644
--- a/src/MCMCregress.cc
+++ b/src/MCMCregress.cc
@@ -66,11 +66,11 @@ static double digamma(double theta, double a, double b) {
*/
template <typename RNGTYPE>
void MCMCregress_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
- const Matrix<>& X, Matrix<>& beta, const Matrix<>& b0,
- const Matrix<>& B0, double c0, double d0,
- unsigned int burnin, unsigned int mcmc, unsigned int thin,
- unsigned int verbose, bool chib,
- Matrix<>& result, double& logmarglike)
+ const Matrix<>& X, Matrix<>& beta, const Matrix<>& b0,
+ const Matrix<>& B0, double c0, double d0,
+ unsigned int burnin, unsigned int mcmc, unsigned int thin,
+ unsigned int verbose, bool chib,
+ Matrix<>& result, double& logmarglike, double& loglike)
{
// define constants and form cross-product matrices
const unsigned int tot_iter = burnin + mcmc; //total iterations
@@ -114,53 +114,67 @@ void MCMCregress_impl (rng<RNGTYPE>& stream, const Matrix<>& Y,
if (chib == 1) {
// marginal likelihood calculation stuff starts here
const double sigma2star = meanc(t(sigmamatrix))(0);
- double sigma2fcdsum = 0.0;
+ const double sigmastar = sqrt(sigma2star);
+ const Matrix<> betastar = t(meanc(t(betamatrix)));
+
+ // step 1
+ Matrix<> sigma2_density(tot_iter, 1);
// second set of Gibbs scans
- for (unsigned int iter = 0; iter < tot_iter; ++iter) {
- double sigma2 = NormIGregress_sigma2_draw (X, Y, beta, c0, d0,
- stream);
- beta = NormNormregress_beta_draw (XpX, XpY, b0, B0, sigma2,
- stream);
+ for (unsigned int iter = 0; iter < nstore; ++iter) {
+ // double sigma2 = sigmamatrix(iter);
+ // beta = NormNormregress_beta_draw (XpX, XpY, b0, B0, sigma2, stream);
+ beta = betamatrix (_, iter);
const Matrix<> e = gaxpy(X, (-1*beta), Y);
const Matrix<> SSE = crossprod (e);
const double c_post = (c0 + X.rows ()) * 0.5;
const double d_post = (d0 + SSE[0]) * 0.5;
- sigma2fcdsum += digamma(sigma2star, c_post, d_post);
-
- // print output to stdout
- if(verbose > 0 && iter % verbose == 0) {
- Rprintf("\n\nMCMCregress (reduced) iteration %i of %i \n",
- (iter+1), tot_iter);
- }
-
- R_CheckUserInterrupt(); // allow user interrupts
+ sigma2_density(iter) = ::exp(digamma(sigma2star, c_post, d_post));
+ R_CheckUserInterrupt();
} // end MCMC loop
+ double pdf_sigma2 = log(mean(sigma2_density));
- double sigma2fcdmean = sigma2fcdsum / static_cast<double>(tot_iter);
-
- const Matrix<> betastar = t(meanc(t(betamatrix)));
- const double sig2_inv = 1.0 / sigma2star;
- const Matrix<> sig_beta = invpd (B0 + XpX * sig2_inv);
- const Matrix<> betahat = sig_beta * gaxpy(B0, b0, XpY*sig2_inv);
- const double logbetafcd = lndmvn(betastar, betahat, sig_beta);
+ // step 2
+ const Matrix<> Bn = invpd (B0 + XpX /sigma2star);
+ const Matrix<> bn = Bn * gaxpy(B0, b0, XpY/sigma2star);
+ double pdf_beta = 0;
+ if (k == 1){
+ pdf_beta = log(dnorm(betastar(0), bn(0), sqrt(Bn(0))));
+ }
+ else {
+ pdf_beta = lndmvn(betastar, bn, Bn);
+ }
// calculate loglikelihood at (betastar, sigma2star)
- double sigmastar = sqrt(sigma2star);
Matrix<> eta = X * betastar;
- double loglike = 0.0;
+ double loglike_sum = 0.0;
for (unsigned int i = 0; i < X.rows(); ++i) {
- loglike += lndnorm(Y(i), eta(i), sigmastar);
+ loglike_sum += lndnorm(Y(i), eta(i), sigmastar);
}
+ loglike = loglike_sum;
// calculate log prior ordinate
- double logprior = log(digamma(sigma2star, c0/2.0, d0/2.0)) +
+ double logprior = 0;
+ if (k == 1){
+ logprior = log(digamma(sigma2star, c0/2.0, d0/2.0)) +
+ log(dnorm(betastar(0), b0(0), 1/sqrt(B0(0))));
+ }
+ else{
+ logprior = log(digamma(sigma2star, c0/2.0, d0/2.0)) +
lndmvn(betastar, b0, invpd(B0));
-
+ }
+
// put pieces together and print the marginal likelihood
- logmarglike = loglike + logprior - logbetafcd - log(sigma2fcdmean);
+ logmarglike = loglike + logprior - pdf_beta - pdf_sigma2;
+ if (verbose >0 ){
+ Rprintf("\nlogmarglike = %10.5f\n", logmarglike);
+ Rprintf("loglike = %10.5f\n", loglike);
+ Rprintf("log_prior = %10.5f\n", logprior);
+ Rprintf("log_beta = %10.5f\n", pdf_beta);
+ Rprintf("log_sigma2 = %10.5f\n", pdf_sigma2);
+ }
}
-
+
result = cbind(t(betamatrix), t(sigmamatrix));
} // end MCMCregress
@@ -176,7 +190,7 @@ extern "C" {
const int *b0row, const int *b0col,
const double *B0data, const int *B0row,
const int *B0col, const double *c0, const double *d0,
- double* logmarglikeholder, const int* chib)
+ double* logmarglikeholder, double* loglikeholder, const int* chib)
{
// pull together Matrix objects
Matrix<> Y(*Yrow, *Ycol, Ydata);
@@ -186,11 +200,13 @@ extern "C" {
Matrix<> B0(*B0row, *B0col, B0data);
double logmarglike;
+ double loglike;
Matrix<> storagematrix;
MCMCPACK_PASSRNG2MODEL(MCMCregress_impl, Y, X, betastart, b0, B0,
*c0, *d0, *burnin, *mcmc, *thin, *verbose,
- *chib, storagematrix, logmarglike);
+ *chib, storagematrix, logmarglike, loglike);
logmarglikeholder[0] = logmarglike;
+ loglikeholder[0] = loglike;
const unsigned int size = *samplerow * *samplecol;
for (unsigned int i = 0; i < size; ++i)
sampledata[i] = storagematrix(i);
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-mcmcpack.git
More information about the debian-science-commits
mailing list