Random samples from the prior of a hyper2 object

Uses Metropolis-Hastings to return random samples from the prior of a hyper2 object

rp(n, H, startp = NULL, fcm = NULL, fcv = NULL, SMALL = 1e-06, l=loglik, fillup=TRUE, ...)

Arguments

H: Object of class hyper2
n: Number of samples
startp: Starting value for the Markov chain, with default NULL being interpreted as starting from the evaluate
fcm,fcv: Constraints as for maxp()
SMALL: Notional small value for numerical stability
l: Log-likelihood function with default loglik()
fillup: Boolean, with default TRUE meaning to return a matrix with the fillup value added, and column names matching the pnames() of argument H
...: Further arguments, currently ignored

Details

Uses the implementation of Metropolis-Hastings from the MCE package to sample from the posterior PDF of a hyper2 object.

If the distribution is Dirichlet, use rdirichlet() to generate random observations: it is much faster, and produces serially independent samples. To return uniform samples, use rp_unif() (documented at dirichlet.Rd).

Value

Returns a matrix, each row being a unit-sum observation.

Author

Robin K. S. Hankin

Note

Function rp() a random sample from a given normalized likelihood function. To return a random likelihood function, use rhyper2().

File inst/ternaryplot_hyper2.Rmd shows how to use Ternary::ternaryPlot() with rp().

Examples

rp(10,icons)
#>              NB         L        PB       THC        OA       WAIS
#>  [1,] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.16666667
#>  [2,] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.16666667
#>  [3,] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.16666667
#>  [4,] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.16666667
#>  [5,] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.16666667
#>  [6,] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.16666667
#>  [7,] 0.1479843 0.1912383 0.2406049 0.1540501 0.1663641 0.09975832
#>  [8,] 0.1479843 0.1912383 0.2406049 0.1540501 0.1663641 0.09975832
#>  [9,] 0.1479843 0.1912383 0.2406049 0.1540501 0.1663641 0.09975832
#> [10,] 0.1479843 0.1912383 0.2406049 0.1540501 0.1663641 0.09975832

plot(loglik(rp(30,icons),icons),type='b')

Random samples from the prior of a `hyper2` object