Merge branch 'master' of https://github.com/stan-dev/docs

Author: mitzimorris

.gitattributes

@@ -0,0 +1 @@
+src/functions-reference/functions_index.qmd linguist-generated=true

.github/ISSUE_TEMPLATE.md

@@ -8,6 +8,3 @@ Describe the issue as clearly as possible.
 
 #### Additional Information:
 Provide any additional information here.
-
-#### Current Version:
-v2.18.0

src/bibtex/all.bib

@@ -1845,3 +1845,52 @@ @article{Timonen+etal:2023:ODE-PSIS
   pages = {e614} 
 }
 
+@article{Vehtari+etal:2024:PSIS,
+  author  = {Aki Vehtari and Daniel Simpson and Andrew Gelman and Yuling Yao and Jonah Gabry},
+  title   = {Pareto smoothed importance sampling},
+  journal = {Journal of Machine Learning Research},
+  year    = {2024},
+  volume  = {25},
+  number  = {72},
+  pages   = {1--58}
+}
+
+@article{Gelman:etal:2020:workflow,
+  title={Bayesian workflow},
+  author={Gelman, Andrew and Vehtari, Aki and Simpson, Daniel and Margossian, Charles C and Carpenter, Bob and Yao, Yuling and Kennedy, Lauren and Gabry, Jonah and B{\"u}rkner, Paul-Christian and Modr{\'a}k, Martin},
+  journal={arXiv preprint arXiv:2011.01808},
+  year={2020}
+}
+
+@article{Magnusson+etal:2024:posteriordb,
+  title={posteriordb: Testing, benchmarking and developing {Bayesian} inference algorithms},
+  author={Magnusson, M{\aa}ns and Torgander, Jakob and B{\"u}rkner, Paul-Christian and Zhang, Lu and Carpenter, Bob and Vehtari, Aki},
+  journal={arXiv preprint arXiv:2407.04967},
+  year={2024}
+
+@article{egozcue+etal:2003,
+  title={Isometric logratio transformations for compositional data analysis},
+  author={Egozcue, Juan Jos{\'e} and Pawlowsky-Glahn, Vera and Mateu-Figueras, Gl{\`o}ria and Barcelo-Vidal, Carles},
+  journal={Mathematical Geology},
+  volume={35},
+  number={3},
+  pages={279--300},
+  year={2003}
+}
+
+@book{filzmoser+etal:2018,
+  title={Geometrical properties of compositional data},
+  author={Filzmoser, Peter and Hron, Karel and Templ, Matthias},
+  booktitle={Applied Compositional Data Analysis: With Worked Examples in R},
+  pages={35--68},
+  year={2018},
+  publisher={Springer}
+}
+
+@misc{seyboldt:2024,
+  author="Seyboldt, Adrian",
+  title="Add ZeroSumNormal distribution",
+  note="pyro-ppl GitHub repository issue \#1751",
+  year = "2024",
+  url ="https://github.com/pyro-ppl/numpyro/pull/1751#issuecomment-1980569811"
+}

src/cmdstan-guide/installation.qmd

@@ -27,10 +27,8 @@ Python bindings to CmdStan seamlessly. Additionally, it provides the
 command `cmdstan_model` to activate the CmdStan makefile from anywhere.
 
 _Note_: This requires that conda has been installed already on your machine.
-You can either install [miniconda](https://docs.conda.io/en/latest/miniconda.html), a free, minimal installer for conda
-or you can get the full [Anaconda](https://docs.continuum.io/anaconda/) system
-which provides graphical installer wizards for [MacOS](https://www.anaconda.com/products/distribution#macos)
-and [Windows](https://www.anaconda.com/products/distribution#windows) users.
+We recommend using the [miniforge](https://github.com/conda-forge/miniforge)
+distribution. 
 
 We recommend installing CmdStan in a new conda environment:
 
@@ -228,7 +226,7 @@ On Linux and macOS:
 
 # default output written to file `output.csv`,
 # default num_samples is 1000, output file should have approx. 1050 lines
-> ls -l output.csv
+> wc -l output.csv
 
 # run the `bin/stansummary utility to summarize parameter estimates
 > bin/stansummary output.csv
@@ -353,36 +351,39 @@ or cluster administrator to install these tools for you.
 
 #### MacOS
 
-To install a C++ development
-environment on a Mac, use Apple's Xcode development environment
-https://developer.apple.com/xcode/.
-
-From the [Xcode home page](https://developer.apple.com/xcode/)
-`View in Mac App Store`.
-
-- From the App Store, click `Install`, enter an Apple ID, and wait
-for Xcode to finish installing.
--  Open the Xcode application, click top-level menu `Preferences`,
-click top-row button `Downloads`, click button for
-`Components`, click on the `Install` button to the right of
-the `Command Line Tools` entry, then wait for it to finish
-installing.
-- Click the top-level menu item `Xcode`, then click item `Quit
-Xcode` to quit.
-
-To test, open the Terminal application and enter:
+To check if you already already have an appropriate toolchain
+installed, open the Terminal application and enter:
 ```
 clang++ --version
 make --version
 ```
 
-If you have installed XCode, but don't have `make`, you can install the
-XCode command-line tools via command:
+If either of these commands prints the message 
+`command not found`, you will need to install Xcode's
+command line tools.
+
+Open the Terminal application and enter:
+
 ```
 xcode-select --install
 ```
 
+Select "Install" in the window that opens.
+
+After the installation completes, you can double check that
+installation was successful by reopening the Terminal and 
+running:
+```
+clang++ --version
+make --version
+```
+
+You can read more about Xcode on its site: 
+[https://developer.apple.com/xcode/](https://developer.apple.com/xcode/)
 
+We don't recommend trying to use the GNU C++ compiler, available via Homebrew,
+based on the number of reports of installation difficulties from Mac users on GitHub
+as well as the Stan forums.
 
 #### Windows {#windows}
 

src/functions-reference/functions_index.qmd

@@ -783,8 +783,10 @@ calculations, but the result is likely to be reduced acceptance
 probabilities and less efficient sampling.
 
 The rounding functions cannot be used as indices to arrays because
-they return real values.  Stan may introduce integer-valued versions
-of these in the future, but as of now, there is no good workaround.
+they return real values. For operations over `data` or in the
+`generated quantities` block, the
+[`to_int()` function](integer-valued_basic_functions.qmd#casting-functions)
+ can be used.
 
 <!-- R; floor; (T x); -->
 \index{{\tt \bfseries floor }!{\tt (T x): R}|hyperpage}
@@ -1601,8 +1603,12 @@ The logarithm of one minus the natural exponentiation of x
 Return the natural logarithm of the difference of the natural
 exponentiation of x and the natural exponentiation of y. \begin{equation*}
 \mathrm{log\_diff\_exp}(x,y) = \begin{cases} \log(\exp(x)-\exp(y)) &
-\text{if } x > y \\[6pt] \textrm{NaN} & \text{otherwise} \end{cases}
+\text{if } +\infty > x \ge y \\[6pt]
+\textrm{NaN} & \text{otherwise} \end{cases}
 \end{equation*}
+
+When x is equal to y, `log_diff_exp(x, y)` returns $-\infty$, consistent with `log(0)` returning $-\infty$. This includes the case in which x and y are both equal to $-\infty$, which corresponds to `log(0 - 0)` because `exp(negative_infinity())` returns 0.
+
 {{< since 2.0 >}}
 
 <!-- R; log_diff_exp; (T1 x, T2 y); -->
@@ -1624,11 +1630,27 @@ proportion theta, defined by \begin{eqnarray*}
 \lambda_1, \ \log(1 - \theta) + \lambda_2\right). \end{eqnarray*}
 {{< since 2.6 >}}
 
-<!-- R; log_mix; (T1 theta, T2 lp1, T3 lp2); -->
-\index{{\tt \bfseries log\_mix }!{\tt (T1 theta, T2 lp1, T3 lp2): real}|hyperpage}
+<!-- R; log_mix; (T1 thetas, T2 lps); -->
+\index{{\tt \bfseries log\_mix }!{\tt (T1 thetas, T2 lps): real}|hyperpage}
+
+`R` **`log_mix`**`(T1 thetas, T2 lps)`<br>\newline
+
+Calculates the log mixture density given `thetas`,
+mixing proportions which should be between 0 and 1 and sum to 1,
+and `lps`, log densities.
+The `lps` variable must be either a 1-d container of the same
+length as `thetas`, or an array of such.
+
+\begin{eqnarray*}
+\mathrm{log\_mix}(\theta, \lambda)
+& = & \log \!\left( \sum_{n=1}^N \theta_n * \exp(\lambda_n) \right) \\[3pt]
+& = & \mathrm{log\_sum\_exp}\!\left(\log(\theta) + \lambda\right).
+\end{eqnarray*}
+
+This is a generalization of the above signature of three arguments to
+more than two densities.
+For example, `log_mix(lambda, lp1, lp2) == log_mix({lambda, 1 - lambda}, {lp1, lp2})`.
 
-`R` **`log_mix`**`(T1 theta, T2 lp1, T3 lp2)`<br>\newline
-Vectorized implementation of the `log_mix` function
 {{< since 2.26 >}}
 
 <!-- R; log_sum_exp; (T1 x, T2 y); -->

src/functions-reference/real-valued_basic_functions.qmd

@@ -586,3 +586,36 @@ The log Poisson probability mass of `y` given the log-rate `alpha + x * beta`.
 The log Poisson probability mass of `y` given the log-rate `alpha + x * beta`
 dropping constant additive terms.
 {{< since 2.25 >}}
+
+## Beta negative binomial distribution {#beta-neg-binomial}
+
+### Probability mass function
+
+If $r \in \mathbb{R}^+$, $\alpha \in \mathbb{R}^+$, and $\beta \in \mathbb{R}^+$, then for $n \in \mathbb{N}$, \begin{equation*}
+\text{BetaNegBinomial}(n|r,\alpha,\beta) = \frac {\Gamma (n+r )}{n!\;\Gamma (r )}
+\frac {\mathrm {B} (\beta+n,\alpha +r )}{\mathrm {B} (\beta,\alpha )}. \end{equation*}
+
+### Distribution statement
+
+`n ~ ` **`beta_neg_binomial`**`(r,alpha,beta)`
+
+Increment target log probability density with `beta_neg_binomial_lupmf(n | r, alpha, beta)`.
+{{< since 2.36 >}}
+<!-- real; beta_neg_binomial ~; -->
+\index{{\tt \bfseries beta\_neg\_binomial }!sampling statement|hyperpage}
+
+### Stan functions
+
+<!-- real; beta_neg_binomial_lpmf; (ints n | reals r, reals alpha, reals beta); -->
+\index{{\tt \bfseries beta\_neg\_binomial\_lpmf }!{\tt (ints n \textbar\ reals r, reals alpha, reals beta): real}|hyperpage}
+
+`real` **`beta_neg_binomial_lpmf`**`(ints n | reals r, reals alpha, reals beta)`<br>\newline
+The log beta negative binomial probability mass of `n` given parameters `r`, `alpha` and `beta`.
+{{< since 2.36 >}}
+
+<!-- real; beta_neg_binomial_lupmf; (ints n | reals r, reals alpha, reals beta); -->
+\index{{\tt \bfseries beta\_neg\_binomial\_lupmf }!{\tt (ints n \textbar\ reals r, reals alpha, reals beta): real}|hyperpage}
+
+`real` **`beta_neg_binomial_lupmf`**`(ints n | reals r, reals alpha, reals beta)`<br>\newline
+The log beta negative binomial probability mass of `n` given parameters `r`, `alpha` and `beta` dropping constant additive terms.
+{{< since 2.36 >}}

src/functions-reference/unbounded_discrete_distributions.qmd

@@ -149,9 +149,9 @@ any of the following.
 
 ```
 int, real, complex, vector, simplex, unit_vector,
-ordered, positive_ordered, row_vector, matrix,
-cholesky_factor_corr, cholesky_factor_cov,
-corr_matrix, cov_matrix, array
+sum_to_zero_vector, ordered, positive_ordered,
+row_vector, matrix, cholesky_factor_corr,
+cholesky_factor_cov, corr_matrix, cov_matrix, array
 ```
 
 The following built in functions are also reserved and
@@ -810,9 +810,9 @@ In addition to single integer indexes, as described in
 [the language indexing section](#language-indexing.section), Stan supports multiple indexing.
 Multiple indexes can be integer arrays of indexes, lower
 bounds, upper bounds, lower and upper bounds, or simply shorthand for
-all of the indexes.  If the upper bound is smaller than the lower bound, 
-the range is empty (unlike, e.g., in R). The upper bound and lower bound can be 
-expressions that evaluate to integer. A complete list of index types is 
+all of the indexes.  If the upper bound is smaller than the lower bound,
+the range is empty (unlike, e.g., in R). The upper bound and lower bound can be
+expressions that evaluate to integer. A complete list of index types is
 given in the following table.
 
 ##### Indexing Options Table {- #index-types-table}
@@ -1078,6 +1078,7 @@ the following table shows the mapping from types to their primitive types.
    | `vector`               | `vector`             |
    | `simplex`              | `vector`             |
    | `unit_vector`          | `vector`             |
+   | `sum_to_zero_vector`   | `vector`             |
    | `ordered`              | `vector`             |
    | `positive_ordered`     | `vector`             |
    | `row_vector`           | `row_vector`         |
@@ -1378,7 +1379,7 @@ model {
 }
 ```
 
-Algebraically, 
+Algebraically,
 [the distribution statement](statements.qmd#distribution-statements.section)
 in the model could be reduced to
 

src/reference-manual/expressions.qmd

@@ -24,6 +24,7 @@
 
 <identifier> ::= IDENTIFIER
                | TRUNCATE
+               | JACOBIAN
 
 <decl_identifier> ::= <identifier>
                     | <reserved_word>
@@ -57,10 +58,13 @@
                   | POSITIVEORDERED
                   | SIMPLEX
                   | UNITVECTOR
+                  | SUMTOZERO
                   | CHOLESKYFACTORCORR
                   | CHOLESKYFACTORCOV
                   | CORRMATRIX
                   | COVMATRIX
+                  | STOCHASTICCOLUMNMATRIX
+                  | STOCHASTICROWMATRIX
                   | PRINT
                   | REJECT
                   | FATAL_ERROR
@@ -165,11 +169,16 @@
                  | POSITIVEORDERED LBRACK <expression> RBRACK
                  | SIMPLEX LBRACK <expression> RBRACK
                  | UNITVECTOR LBRACK <expression> RBRACK
+                 | SUMTOZERO LBRACK <expression> RBRACK
                  | CHOLESKYFACTORCORR LBRACK <expression> RBRACK
                  | CHOLESKYFACTORCOV LBRACK <expression> [COMMA <expression>]
                    RBRACK
                  | CORRMATRIX LBRACK <expression> RBRACK
                  | COVMATRIX LBRACK <expression> RBRACK
+                 | STOCHASTICCOLUMNMATRIX LBRACK <expression> COMMA
+                   <expression> RBRACK
+                 | STOCHASTICROWMATRIX LBRACK <expression> COMMA <expression>
+                   RBRACK
 
 <type_constraint> ::= [LABRACK <range> RABRACK]
                     | LABRACK <offset_mult> RABRACK

src/reference-manual/grammar.txt

@@ -4,7 +4,7 @@ pagetitle: Pathfinder
 
 # Pathfinder
 
-Stan supports the Pathfinder algorithm @zhang_pathfinder:2022.
+Stan supports the Pathfinder algorithm [@zhang_pathfinder:2022].
 Pathfinder is a variational method for approximately
 sampling from differentiable log densities.  Starting from a random
 initialization, Pathfinder locates normal approximations to the target
@@ -22,6 +22,31 @@ the problem of L-BFGS getting stuck at local optima or in saddle points on plate
 Compared to ADVI and short dynamic HMC runs, Pathfinder
 requires one to two orders of magnitude fewer log density and gradient
 evaluations, with greater reductions for more challenging posteriors.
-While the evaluations in @zhang_pathfinder:2022 found that
-single-path and multi-path Pathfinder outperform ADVI for most of the models in the PosteriorDB evaluation set,
+While the evaluations by @zhang_pathfinder:2022 found that
+single-path and multi-path Pathfinder outperform ADVI for most of the models in the PosteriorDB [@Magnusson+etal:2024:posteriordb] evaluation set,
 we recognize the need for further experiments on a wider range of models.
+
+## Diagnosing Pathfinder
+
+Pathfinder diagnoses the accuracy of the approximation by computing the density ratio of the true posterior and 
+the approximation and using Pareto-$\hat{k}$ diagnostic [@Vehtari+etal:2024:PSIS] to assess whether these ratios can
+be used to improve the approximation via resampling. The
+normalization for the posterior can be  estimated reliably [@Vehtari+etal:2024:PSIS, Section 3], which is the
+first requirement for reliable resampling.  If estimated Pareto-$\hat{k}$ for the ratios is smaller than 0.7,
+there is still need to further diagnose reliability of importance sampling estimate for all quantities of interest [@Vehtari+etal:2024:PSIS, Section 2.2]. If estimated Pareto-$\hat{k}$ is larger than 0.7, then the 
+estimate for the normalization is unreliable and any Monte Carlo estimate may have a big error. The resampled draws
+can still contain some useful information about the location and shape of the posterior which can be used in early
+parts of Bayesian workflow [@Gelman:etal:2020:workflow].
+
+## Using Pathfinder for initializing MCMC
+
+If estimated Pareto-$\hat{k}$ for the ratios is smaller than 0.7, the resampled posterior draws are almost as
+good for initializing MCMC as would independent draws from the posterior be. If estimated Pareto-$\hat{k}$ for the 
+ratios is larger than 0.7, the Pathfinder draws are not reliable for posterior inference directly, but they are still 
+very likely better for initializing MCMC than random draws from an arbitrary pre-defined distribution (e.g. uniform from 
+-2 to 2 used by Stan by default). If Pareto-$\hat{k}$ is larger than 0.7, it is likely that one of the ratios is much bigger
+than others and the default resampling with replacement would produce copies of one unique draw. For initializing several
+Markov chains, it is better to use resampling without replacement to guarantee unique initialization for each chain. At the
+moment Stan allows turning off the resampling completely, and then the resampling without replacement can be done outside of
+Stan.
+