sampling from a mixture of truncated normal distributions
314 views
Skip to first unread message
Joseph Catanzarite
Sep 10, 2016, 8:27:40 PM9/10/16
to Stan users mailing list
Hi,
I have a problem where I want to sample from a distribution that is loosely speaking a "two-sided Gaussian" -- it is composed of the left half of a gaussian of width sigma2 and the right half of a gaussian of width sigma 1. The right half is to be scaled so that its peak is the same as that of the left. So the distribution looks like a Gaussian that is asymmetric about its peak.
I think a way to do this is to use a mixture of two truncated Gaussians with mixture fractions equal to sigma1/(sigma1+sigma2) and sigma2/(sigma1+sigma2)
The distribution for the left side is y ~ normal(mu,sigma1) T[,mu], so that it has support only on the left side of the peak
The distribution for the right side is y ~ normal(mu,sigma2) T[mu,], so that it has support only on the right side of the peak
I implemented the mixture by drawing a uniform random variable called mixingParam; if it's less than sigma1/(sigma1+sigma2), I draw from the first distribution and otherwise I draw from the second.
The complication is that y and mu and sigma are really vectors.
After reading the sections in the Stan manual about truncated distributions and mixtures, I implemented this in my Stan model.
The PyStan code is included below: I've included the entire Stan model, but the relevant part is in the parameters and model blocks.
The model compiles, but it doesn't sample, returning the error below.
Any help to correctly implement this sampling would be appreciated!
Cheers,
Joe Catanzarite
from pystan import StanModel
stan_model_occ = """
functions {
real[,] population_model_insolation(int nxGrid, int nRadiusGrid,
real lnf0, real alpha, real beta, real[] x_grid, real[] planet_radius_grid,
real[] planet_radius_rng, real[] x_rng) {
/* Population inference with an independent power law model
Using modified code above that computes
completeness in the parameter space of [ insolation , planet radius ],
A double power law model for the population.
:param nxGrid: size of insolation grid
:param nRadiusGrid: size of radius grid
:param lnf0: normalization,
:param beta: exponent of insolation power law,
:param alpha: exponent of radius power law
:param x_grid: //1D grid
:param planet_radius_grid: //1D grid
:param x_rng: // lower and upper boundary values of insolation grid
:param radius_rng: // lower and upper boundary values of radius grid
*/
// Declare internal variables
real powerLawNumberDensity[nRadiusGrid,nxGrid];
// Loop over phase space bins, compute number density in each bin,
// from the power law using the input lnf0, alpha and beta parameters
for (j in 1:nxGrid)
for (i in 1:nRadiusGrid) {
powerLawNumberDensity[i,j] <- ( exp( lnf0 ) * ( alpha+1 ) * ( beta+1 ) * ( planet_radius_grid[i]^alpha ) * ( x_grid[j]^beta ) ) / ( ( planet_radius_rng[2]^( alpha+1 ) - ( planet_radius_rng[1]^( alpha+1 ) ) ) * ( ( x_rng[2]^( beta+1 ) ) - ( x_rng[1]^( beta+1 ) ) ) );
}
return powerLawNumberDensity;
}
real lnlike(real alpha, real beta, real lnf0, int nKois, real[] planet_radius_rng,
real[] x_rng, real[,] vol, real[] x_grid, real[] planet_radius_grid,
real[] koi_x, real[] koi_rps_true, int nRadiusGrid, int nxGrid, real e, real[,] new_completeness ) {
/* log likelihood function
:param alpha:
:param beta:
:param lnf0:
:param nKois
:param planet_radius_rng
:param x_rng
:param vol:
:param x_grid:
:param planet_radius_grid:
:param koi_x:
:param koi_rps_true:
:param nRadiusGrid
:param nxGrid
:param e
:param new_completeness
*/
// Declare internal variables
real pop[nRadiusGrid,nxGrid];
real pop1[nRadiusGrid,nxGrid];
real pop2[nRadiusGrid,nxGrid];
real populationDensity[nRadiusGrid-1,nxGrid-1];
real sumExpectedCounts;
real logPop[nKois];
real lnLikelihood;
real sumLogPop;
real lnLikelihood0;
real popModelInsol0[nRadiusGrid,nxGrid];
real popModelInsol[nKois,nKois];
// Compute the population model for a grid over [ insolation, radius ] parameter space
popModelInsol0 <- population_model_insolation(nxGrid, nRadiusGrid, lnf0, alpha, beta, x_grid, planet_radius_grid, planet_radius_rng, x_rng);
// Compute the expected density of detections per bin in the [ insolation, radius ] parameter space,
// based on the population model and the completeness
for (i in 1:nRadiusGrid)
for (j in 1: nxGrid) {
//sum_completeness_0[i,j] <- new_completeness[j,i];
pop[i,j] <- popModelInsol0[i,j] * new_completeness[j,i];
}
// Initialize
sumExpectedCounts <- 0.0;
sumLogPop <-0.0;
for (i in 1:nRadiusGrid-1)
for (j in 1: nxGrid-1) {
// Compute the number density of detections on a grid of bin centers
// Figured out the proper syntax for translation of the Python line below into Stan
// pop <- 0.5*(pop[:-1, :-1] + pop[1:, 1:]) is python syntax from DFM notebook.
// Note that indices start with 1 in Stan
pop1[i,j] <- pop[i,j];
pop2[i,j] <- pop[i+1,j+1];
populationDensity[i,j] <- 0.5*(pop1[i,j]+pop2[i,j]);
// Integrate the planet counts over the phase space volume
sumExpectedCounts <- sumExpectedCounts + (populationDensity[i,j]*vol[j,i]);
}
// Population model evaluated over the list of detected planets
//
popModelInsol <- population_model_insolation(nKois, nKois, lnf0, alpha, beta, koi_x, koi_rps_true, planet_radius_rng, x_rng);
for (k in 1:nKois) {
logPop[k] <- log( popModelInsol[k,k] );
sumLogPop <- sumLogPop + logPop[k];
}
// Combine the two terms to form the log-likelihood function for the Inhomogeneous Poisson Process
lnLikelihood <- sumLogPop - sumExpectedCounts;
// Catch error: if lnLikelihood has bad value, return a very small number
lnLikelihood0 <- if_else( is_inf(lnLikelihood) || is_nan(lnLikelihood) , -1e15, lnLikelihood );
return lnLikelihood0;
}
}
// end function block
data {
// Declare precomputed variables as data
int<lower=1> nKois;
int<lower=1> nRadiusGrid;
int<lower=1> nxGrid;
real x_rng[2];
real planet_radius_rng[2];
real x_grid[nxGrid];
real planet_radius_grid[nRadiusGrid];
real vol[nxGrid-1,nRadiusGrid-1];
real koi_x[nKois];
real koi_rps_obs[nKois];
real koi_rps_err1[nKois];
real koi_rps_err2[nKois];
real koi_rps_obs_lower[nKois];
real koi_rps_obs_upper[nKois];
real e;
real new_completeness[nxGrid,nRadiusGrid];
}
parameters {
// For v1: only compute population level parameter posteriors.
real<lower=-10,upper=5> alpha;
real<lower=-5,upper=5> beta;
real<lower=-5,upper=5> lnf0;
# Ad hoc lower bound constraint on planet radius
real<lower=0.5> koi_rps_true[nKois];
# Mixing parameter
real<lower=0,upper=1> mixingParam;
}
model {
// Hyperpriors
alpha ~ uniform(-10,5);
beta ~ uniform(-5,5);
lnf0 ~ uniform(-5,5);
mixingParam ~ uniform(0,1);
// Sample the planet radius errors
for (n in 1:nKois)
// Draw from either left-side Gaussian with prob f1 or right-side Gaussian with prob f2
// Where mixing weights are
// f1 = koi_rps_err2[n]/(koi_rps_err1[n] + koi_rps_err2[n])
// f2 = koi_rps_err1[n]/(koi_rps_err1[n] + koi_rps_err2[n])
if(mixingParam < koi_rps_err2[n]/(koi_rps_err1[n] + koi_rps_err2[n]) )
// Draw from left-side Gaussian, of width koi_rps_err2
koi_rps_true[n] ~ normal(koi_rps_obs[n],koi_rps_err2[n]) T[,koi_rps_obs[n]];
else
// Draw from right-side Gaussian, of width koi_rps_err1
koi_rps_true[n] ~ normal(koi_rps_obs[n],koi_rps_err1[n]) T[koi_rps_obs[n],];
//koi_rps_true[n] ~ normal(koi_rps_obs[n], koi_rps_err2[n]);
// Sample using lnLikelihood function
increment_log_prob(lnlike(alpha, beta, lnf0, nKois, planet_radius_rng, x_rng, vol, x_grid, planet_radius_grid, koi_x, koi_rps_true, nRadiusGrid, nxGrid, e, new_completeness));
}
"""
# Compiled Stan Model
sm = StanModel(model_code=stan_model_occ)
data = {'new_completeness':new_completeness, 'nKois':nKois,'nRadiusGrid':nRadiusGrid, 'nxGrid':nxGrid,
'x_rng':x_rng, 'planet_radius_rng':planet_radius_rng,'x_grid':x_grid,
'planet_radius_grid':planet_radius_grid,
'vol':vol,'koi_x':koi_x, 'koi_rps_obs':koi_rps, 'koi_rps_err1':koi_rps_err1, 'koi_rps_err2':abs(koi_rps_err2),'e':e,
'koi_rps_obs_lower':koi_rps_obs_lower,'koi_rps_obs_upper':koi_rps_obs_upper}
init = [{'lnf0':-0.1,'alpha':-2.2,'beta':-0.5}]
# Set total number of samples per chain and number of chains
# number of production samples is half of the total number of samples;
# the other half are warmup samples.
nTotalPerChain = 500
nSamplesPerChain = int(round(nTotalPerChain/2))
nChains = 5
fit = sm.sampling(data=data, iter=nTotalPerChain, chains=nChains, n_jobs=-1)
# Return a dictionary of arrays of posterior samples
la = fit.extract(permuted=True)
alpha = la['alpha']
beta = la['beta']
lnf0 = la['lnf0']
a = fit.extract(permuted=False)
print(fit)
I have a problem where I want to sample from a distribution that is loosely speaking a "two-sided Gaussian" -- it is composed of the left half of a gaussian of width sigma2 and the right half of a gaussian of width sigma 1. The right half is to be scaled so that its peak is the same as that of the left. So the distribution looks like a Gaussian that is asymmetric about its peak.
I think a way to do this is to use a mixture of two truncated Gaussians with mixture fractions equal to sigma1/(sigma1+sigma2) and sigma2/(sigma1+sigma2)
The distribution for the left side is y ~ normal(mu,sigma1) T[,mu], so that it has support only on the left side of the peak
The distribution for the right side is y ~ normal(mu,sigma2) T[mu,], so that it has support only on the right side of the peak
I implemented the mixture by drawing a uniform random variable called mixingParam; if it's less than sigma1/(sigma1+sigma2), I draw from the first distribution and otherwise I draw from the second.
The complication is that y and mu and sigma are really vectors.
After reading the sections in the Stan manual about truncated distributions and mixtures, I implemented this in my Stan model.
The PyStan code is included below: I've included the entire Stan model, but the relevant part is in the parameters and model blocks.
The model compiles, but it doesn't sample, returning the error below.
Any help to correctly implement this sampling would be appreciated!
Cheers,
Joe Catanzarite
--------------------------------------------------------------------------- RuntimeError Traceback (most recent call last) <ipython-input-28-9777225edb81> in <module>() 25 nSamplesPerChain = int(round(nTotalPerChain/2)) 26 nChains = 5 ---> 27 fit = sm.sampling(data=data, iter=nTotalPerChain, chains=nChains, n_jobs=-1) 28 #fit = sm.sampling(data=data, iter=2000, chains=5, init=init, n_jobs=-1) 29 #fit = pystan.stan(model_code=stan_model_occ, data=data, iter=100, chains=2, n_jobs=-1, verbose=True); /Users/jcatanza/anaconda/lib/python2.7/site-packages/pystan/model.pyc in sampling(self, data, pars, chains, iter, warmup, thin, seed, init, sample_file, diagnostic_file, verbose, algorithm, control, n_jobs, **kwargs) 719 call_sampler_args = izip(itertools.repeat(data), args_list, itertools.repeat(pars)) 720 call_sampler_star = self.module._call_sampler_star --> 721 ret_and_samples = _map_parallel(call_sampler_star, call_sampler_args, n_jobs) 722 samples = [smpl for _, smpl in ret_and_samples] 723 /Users/jcatanza/anaconda/lib/python2.7/site-packages/pystan/model.pyc in _map_parallel(function, args, n_jobs) 79 n_jobs = None 80 pool = multiprocessing.Pool(processes=n_jobs) ---> 81 map_result = pool.map(function, args) 82 pool.close() 83 pool.join() /Users/jcatanza/anaconda/lib/python2.7/multiprocessing/pool.pyc in map(self, func, iterable, chunksize) 249 ''' 250 assert self._state == RUN --> 251 return self.map_async(func, iterable, chunksize).get() 252 253 def imap(self, func, iterable, chunksize=1): /Users/jcatanza/anaconda/lib/python2.7/multiprocessing/pool.pyc in get(self, timeout) 565 return self._value 566 else: --> 567 raise self._value 568 569 def _set(self, i, obj): RuntimeError: Rejecting initial value: Log probability evaluates to log(0), i.e. negative infinity. Stan can't start sampling from this initial value. Rejecting initial value:==========================================
from pystan import StanModel
stan_model_occ = """
functions {
real[,] population_model_insolation(int nxGrid, int nRadiusGrid,
real lnf0, real alpha, real beta, real[] x_grid, real[] planet_radius_grid,
real[] planet_radius_rng, real[] x_rng) {
/* Population inference with an independent power law model
Using modified code above that computes
completeness in the parameter space of [ insolation , planet radius ],
A double power law model for the population.
:param nxGrid: size of insolation grid
:param nRadiusGrid: size of radius grid
:param lnf0: normalization,
:param beta: exponent of insolation power law,
:param alpha: exponent of radius power law
:param x_grid: //1D grid
:param planet_radius_grid: //1D grid
:param x_rng: // lower and upper boundary values of insolation grid
:param radius_rng: // lower and upper boundary values of radius grid
*/
// Declare internal variables
real powerLawNumberDensity[nRadiusGrid,nxGrid];
// Loop over phase space bins, compute number density in each bin,
// from the power law using the input lnf0, alpha and beta parameters
for (j in 1:nxGrid)
for (i in 1:nRadiusGrid) {
powerLawNumberDensity[i,j] <- ( exp( lnf0 ) * ( alpha+1 ) * ( beta+1 ) * ( planet_radius_grid[i]^alpha ) * ( x_grid[j]^beta ) ) / ( ( planet_radius_rng[2]^( alpha+1 ) - ( planet_radius_rng[1]^( alpha+1 ) ) ) * ( ( x_rng[2]^( beta+1 ) ) - ( x_rng[1]^( beta+1 ) ) ) );
}
return powerLawNumberDensity;
}
real lnlike(real alpha, real beta, real lnf0, int nKois, real[] planet_radius_rng,
real[] x_rng, real[,] vol, real[] x_grid, real[] planet_radius_grid,
real[] koi_x, real[] koi_rps_true, int nRadiusGrid, int nxGrid, real e, real[,] new_completeness ) {
/* log likelihood function
:param alpha:
:param beta:
:param lnf0:
:param nKois
:param planet_radius_rng
:param x_rng
:param vol:
:param x_grid:
:param planet_radius_grid:
:param koi_x:
:param koi_rps_true:
:param nRadiusGrid
:param nxGrid
:param e
:param new_completeness
*/
// Declare internal variables
real pop[nRadiusGrid,nxGrid];
real pop1[nRadiusGrid,nxGrid];
real pop2[nRadiusGrid,nxGrid];
real populationDensity[nRadiusGrid-1,nxGrid-1];
real sumExpectedCounts;
real logPop[nKois];
real lnLikelihood;
real sumLogPop;
real lnLikelihood0;
real popModelInsol0[nRadiusGrid,nxGrid];
real popModelInsol[nKois,nKois];
// Compute the population model for a grid over [ insolation, radius ] parameter space
popModelInsol0 <- population_model_insolation(nxGrid, nRadiusGrid, lnf0, alpha, beta, x_grid, planet_radius_grid, planet_radius_rng, x_rng);
// Compute the expected density of detections per bin in the [ insolation, radius ] parameter space,
// based on the population model and the completeness
for (i in 1:nRadiusGrid)
for (j in 1: nxGrid) {
//sum_completeness_0[i,j] <- new_completeness[j,i];
pop[i,j] <- popModelInsol0[i,j] * new_completeness[j,i];
}
// Initialize
sumExpectedCounts <- 0.0;
sumLogPop <-0.0;
for (i in 1:nRadiusGrid-1)
for (j in 1: nxGrid-1) {
// Compute the number density of detections on a grid of bin centers
// Figured out the proper syntax for translation of the Python line below into Stan
// pop <- 0.5*(pop[:-1, :-1] + pop[1:, 1:]) is python syntax from DFM notebook.
// Note that indices start with 1 in Stan
pop1[i,j] <- pop[i,j];
pop2[i,j] <- pop[i+1,j+1];
populationDensity[i,j] <- 0.5*(pop1[i,j]+pop2[i,j]);
// Integrate the planet counts over the phase space volume
sumExpectedCounts <- sumExpectedCounts + (populationDensity[i,j]*vol[j,i]);
}
// Population model evaluated over the list of detected planets
//
popModelInsol <- population_model_insolation(nKois, nKois, lnf0, alpha, beta, koi_x, koi_rps_true, planet_radius_rng, x_rng);
for (k in 1:nKois) {
logPop[k] <- log( popModelInsol[k,k] );
sumLogPop <- sumLogPop + logPop[k];
}
// Combine the two terms to form the log-likelihood function for the Inhomogeneous Poisson Process
lnLikelihood <- sumLogPop - sumExpectedCounts;
// Catch error: if lnLikelihood has bad value, return a very small number
lnLikelihood0 <- if_else( is_inf(lnLikelihood) || is_nan(lnLikelihood) , -1e15, lnLikelihood );
return lnLikelihood0;
}
}
// end function block
data {
// Declare precomputed variables as data
int<lower=1> nKois;
int<lower=1> nRadiusGrid;
int<lower=1> nxGrid;
real x_rng[2];
real planet_radius_rng[2];
real x_grid[nxGrid];
real planet_radius_grid[nRadiusGrid];
real vol[nxGrid-1,nRadiusGrid-1];
real koi_x[nKois];
real koi_rps_obs[nKois];
real koi_rps_err1[nKois];
real koi_rps_err2[nKois];
real koi_rps_obs_lower[nKois];
real koi_rps_obs_upper[nKois];
real e;
real new_completeness[nxGrid,nRadiusGrid];
}
parameters {
// For v1: only compute population level parameter posteriors.
real<lower=-10,upper=5> alpha;
real<lower=-5,upper=5> beta;
real<lower=-5,upper=5> lnf0;
# Ad hoc lower bound constraint on planet radius
real<lower=0.5> koi_rps_true[nKois];
# Mixing parameter
real<lower=0,upper=1> mixingParam;
}
model {
// Hyperpriors
alpha ~ uniform(-10,5);
beta ~ uniform(-5,5);
lnf0 ~ uniform(-5,5);
mixingParam ~ uniform(0,1);
// Sample the planet radius errors
for (n in 1:nKois)
// Draw from either left-side Gaussian with prob f1 or right-side Gaussian with prob f2
// Where mixing weights are
// f1 = koi_rps_err2[n]/(koi_rps_err1[n] + koi_rps_err2[n])
// f2 = koi_rps_err1[n]/(koi_rps_err1[n] + koi_rps_err2[n])
if(mixingParam < koi_rps_err2[n]/(koi_rps_err1[n] + koi_rps_err2[n]) )
// Draw from left-side Gaussian, of width koi_rps_err2
koi_rps_true[n] ~ normal(koi_rps_obs[n],koi_rps_err2[n]) T[,koi_rps_obs[n]];
else
// Draw from right-side Gaussian, of width koi_rps_err1
koi_rps_true[n] ~ normal(koi_rps_obs[n],koi_rps_err1[n]) T[koi_rps_obs[n],];
//koi_rps_true[n] ~ normal(koi_rps_obs[n], koi_rps_err2[n]);
// Sample using lnLikelihood function
increment_log_prob(lnlike(alpha, beta, lnf0, nKois, planet_radius_rng, x_rng, vol, x_grid, planet_radius_grid, koi_x, koi_rps_true, nRadiusGrid, nxGrid, e, new_completeness));
}
"""
# Compiled Stan Model
sm = StanModel(model_code=stan_model_occ)
data = {'new_completeness':new_completeness, 'nKois':nKois,'nRadiusGrid':nRadiusGrid, 'nxGrid':nxGrid,
'x_rng':x_rng, 'planet_radius_rng':planet_radius_rng,'x_grid':x_grid,
'planet_radius_grid':planet_radius_grid,
'vol':vol,'koi_x':koi_x, 'koi_rps_obs':koi_rps, 'koi_rps_err1':koi_rps_err1, 'koi_rps_err2':abs(koi_rps_err2),'e':e,
'koi_rps_obs_lower':koi_rps_obs_lower,'koi_rps_obs_upper':koi_rps_obs_upper}
init = [{'lnf0':-0.1,'alpha':-2.2,'beta':-0.5}]
# Set total number of samples per chain and number of chains
# number of production samples is half of the total number of samples;
# the other half are warmup samples.
nTotalPerChain = 500
nSamplesPerChain = int(round(nTotalPerChain/2))
nChains = 5
fit = sm.sampling(data=data, iter=nTotalPerChain, chains=nChains, n_jobs=-1)
# Return a dictionary of arrays of posterior samples
la = fit.extract(permuted=True)
alpha = la['alpha']
beta = la['beta']
lnf0 = la['lnf0']
a = fit.extract(permuted=False)
print(fit)
Daniel Lee
Sep 10, 2016, 10:16:00 PM9/10/16
to stan-users mailing list
Hi Joe,
Does any Stan program run for you in PyStan? I wanted to make sure before trying to see if there's something inside you Stan program that is preventing it from sampling.
Daniel
--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Joseph Catanzarite
Sep 10, 2016, 10:23:46 PM9/10/16
to stan-...@googlegroups.com
Hi Dan
Good question! Yes the previous version of the code (which just samples rp_true from a normal distribution) compiles and runs fine.
Cheers,
Joe
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/fChtE7RDds8/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.
Daniel Lee
Sep 10, 2016, 10:37:13 PM9/10/16
to stan-users mailing list
Hi Joe,
Sorry -- I missed the relevant error message from the original post:
RuntimeError: Rejecting initial value: Log probability evaluates to log(0), i.e. negative infinity. Stan can't start sampling from this initial value. Rejecting initial value:
The first thing I'd do is remove the initialization. I haven't really needed to provide initial values for Stan when sampling. And, it'll help diagnose problems with multiple chains if you allow random initialization.
You might also want to check that there is support (probability > 0) for the parameter values. If not, then you might need constraints on your parameters.
Daniel
joseph catanzarite
Sep 11, 2016, 6:02:16 PM9/11/16
to stan-...@googlegroups.com, jcat...@gmail.com, Catanzarite, Joseph H (ARC-383D)[SETI INSTITUTE]
Thanks, Dan
I turned off initialization, but got the same result.
To simplify the question as far as possible:
is the code snippet below a valid way to use truncated gaussians
to
draw from an asymmetric gaussian distribution, centered on y_obs[n] and with width err2 on the left and err1 on the right?
The support for y_true[n] is real<lower=0>, and all values of y_obs[n] are positive:
for (n in 1:N)
// Draw from either left-side Gaussian with prob f1 or
right-side Gaussian with prob f2
// Where mixing weights are
// f1 = err2[n]/(err1[n] + err2[n])
// f2 = err1[n]/(err1[n] + err2[n])
mixingParam ~ uniform(0,1)
if(mixingParam < err2[n]/(err1[n] + err2[n]) )
// Draw from left-side Gaussian, of width err2
y_true[n] ~ normal(y_obs[n],err2[n]) T[,y_obs[n]];
else
// Draw from right-side Gaussian, of width err1
y_true[n] ~ normal(y_obs[n],err1[n]) T[y_obs[n],];
Cheers,
Joe
if(mixingParam < err2[n]/(err1[n] + err2[n]) )
// Draw from left-side Gaussian, of width err2
y_true[n] ~ normal(y_obs[n],err2[n]) T[,y_obs[n]];
else
// Draw from right-side Gaussian, of width err1
y_true[n] ~ normal(y_obs[n],err1[n]) T[y_obs[n],];
Cheers,
Joe
-- You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group. To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/fChtE7RDds8/unsubscribe. To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com. To post to this group, send email to stan-...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
Reply all
Reply to author
Forward
0 new messages