# Update on using Poisson zero trick to model skewed normal distribution Classic List Threaded 12 messages Open this post in threaded view
|

## Update on using Poisson zero trick to model skewed normal distribution

 It took some time, but I was able to switch to JAGS and figure out how to use the Poisson zeros trick to model a distribution that isn't included in JAGS - the skewed normal distribution.  The code runs and returns a result. In order to make sure that it is working correctly, I've used rsn (random skewed normal) to create a data set with known parameters and then I see how closely the model comes to yielding those known parameters. I've used sample sizes of 50 and 500 and modified the priors and modeling parameters multiple times. (I've even "cheated" and set the initial values to the correct values, etc.) But the outputs do not seem close to me. My parameters are: Location = 0 Scale = 1 Shape = 0.5 Across multiple modeling runs, here are the parameter values I'm getting: Location: Average = 0.63 (Std = 0.29) Scale: Average = 0.48 (Std = 0.08) Shape: Average = 0.98 (Std = 0.12) Any thoughts? For example, are these values actually the best I should expect? Of course, it's possible that my equation is wrong, although I've checked it many times (and had someone else check it). So, I doubt this is the problem. But I'm willing to consider anything... Any/all advice/suggestions would be most welcome! Thanks! -Gwen
Open this post in threaded view
|

## Re: Update on using Poisson zero trick to model skewed normal distribution

 Administrator Without seeing the complete programming code, it is difficult to know what's really going on. But two general questions come to mind:First, are you sure that your formula for the skewed normal (used to specify the log-likelihood for the Poisson zero) is parameterized the same way as the rsn random value generator you used? If they are parameterized differently, the results won't match. Second, if the posterior distributions are skewed, then the means of the posterior won't necessarily match the generating values very well (I think). Try looking at the modes of the posteriors instead of the means. On Fri, Aug 9, 2013 at 8:52 AM, Gwendolyn Campbell [via Doing Bayesian Data Analysis] wrote: It took some time, but I was able to switch to JAGS and figure out how to use the Poisson zeros trick to model a distribution that isn't included in JAGS - the skewed normal distribution.  The code runs and returns a result. In order to make sure that it is working correctly, I've used rsn (random skewed normal) to create a data set with known parameters and then I see how closely the model comes to yielding those known parameters. I've used sample sizes of 50 and 500 and modified the priors and modeling parameters multiple times. (I've even "cheated" and set the initial values to the correct values, etc.) But the outputs do not seem close to me. My parameters are: Location = 0 Scale = 1 Shape = 0.5 Across multiple modeling runs, here are the parameter values I'm getting: Location: Average = 0.63 (Std = 0.29) Scale: Average = 0.48 (Std = 0.08) Shape: Average = 0.98 (Std = 0.12) Any thoughts? For example, are these values actually the best I should expect? Of course, it's possible that my equation is wrong, although I've checked it many times (and had someone else check it). So, I doubt this is the problem. But I'm willing to consider anything... Any/all advice/suggestions would be most welcome! Thanks! -Gwen If you reply to this email, your message will be added to the discussion below: http://doing-bayesian-data-analysis.12272.x6.nabble.com/Update-on-using-Poisson-zero-trick-to-model-skewed-normal-distribution-tp5000706.html To start a new topic under Doing Bayesian Data Analysis, email [hidden email] To unsubscribe from Doing Bayesian Data Analysis, click here. NAML
Open this post in threaded view
|

## Re: Update on using Poisson zero trick to model skewed normal distribution

 Thanks, John, for getting back to me so quickly.  I had thought of the issue of mis-matched parameters, but it looks like the pdf I used for the skewed normal distribution (from wikipedia) and the rsn function use the same parameters...  They certainly use the same parameter names and I aligned them in my code.  I don't think this is it.   I'm not 100% sure I understand your second point.  It may be that I didn't explain clearly how I got the numbers that I originally posted.  Using one data set, I ran the model code 15 times.  Each time I changed something about the modeling parameters - maybe the number of chains, maybe the initial values of the priors, etc.  Each time I ran it, it output a value for each parameter (location, scale and shape) of the skewed normal distribution that it converged on.  The averages that I reported were simply the averages across those 15 runs.  Does that make sense?  Is it relevant to your comment?  If not, then I'm sorry - let me know and I'll try to make sense out of your comment again.   Thanks, Gwen
Open this post in threaded view
|

## Re: Update on using Poisson zero trick to model skewed normal distribution

 In reply to this post by John K. Kruschke PS - Here's the code: (with alex=location, bob=scale and chip=shape) -------------------------------------- require(rjags)         rm(list=ls(all=TRUE)) #------------------------------------------------------------------------------ # THE MODEL. modelstring = " model { for (i in 1:50) {    zeros[i] ~ dpois(phi[i])    phi[i] <- log(2.5)+log(abs(pow(bob,2)/pow(chip,2)*exp(pow(chip,2))*(exp(pow(chip,2))-1)))+log(abs(bob-chip*(x[i]-alex)))+pow((-1/chip)*log(abs(1-((1/bob)*(chip)*(x[i]-alex))))-(alex-bob/chip*(exp((pow(chip,2))/2))-1),2)/(2*pow((pow(bob,2)/pow(chip,2)*exp(pow(chip,2))*(exp((pow(chip,2)))-1)),2))    }      alex ~ dunif(-2,2)      bob ~ dunif(0.01,5)      chip ~ dunif(0.01,5) } " # close quote for modelstring # Write model to a file, and send to BUGS: writeLines(modelstring,con="model.txt") #------------------------------------------------------------------------------ # THE DATA. # These data were created with parameter values [n=50, location=0, scale=1, shape=0.5]  dataList = list (    zeros = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,              0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,              0, 0, 0, 0, 0, 0, 0, 0, 0, 0),    x = c(-1.65534362, -1.33096199, -1.15361727, -1.14141452, -1.01578541,          -0.88231001, -0.86824978, -0.66212271, -0.59999216, -0.47238434,          -0.46603680, -0.38355513, -0.28087996, -0.19715436, -0.10505644,          0.01120428, 0.21384728, 0.22846082, 0.23366761, 0.26649210,          0.28282272, 0.32582633, 0.34812280, 0.42625461, 0.44375500,          0.54342404, 0.55071416, 0.57397426, 0.57490333, 0.66954006,          0.74401185, 0.76160107, 0.77146795, 0.78891622, 0.79762613,          0.92826248, 1.01350578, 1.07251810, 1.10103133, 1.24612441,          1.27405717, 1.28771869, 1.29966906, 1.35210939, 1.42059526,          1.50968752, 1.55416321, 1.55942491, 1.69682208, 1.92283079)  )  initsList = list(    alex = c(0),    bob = c(1),    chip = c(0.5)  ) #------------------------------------------------------------------------------ # RUN THE CHAINS parameters = c ("alex", "bob", "chip") adaptSteps = 1000 burnInSteps = 1000 nChains = 7 numSavedSteps = 50000 thinSteps = 1 nIter = ceiling ( (numSavedSteps * thinSteps) / nChains) jagsModel = jags.model( "model.txt" , data=dataList , inits=initsList ,                         n.chains=nChains , n.adapt=adaptSteps ) cat("Burning in MCMC chain...\n") update(jagsModel , n.iter=burnInSteps) cat("Sampling final MCMC chain...\n") codaSamples = coda.samples ( jagsModel , variable.names = parameters, n.iter=nIter, thin = thinSteps) #------------------------------------------------------------------------------ # EXAMINE THE RESULTS mcmcChain = as.matrix( codaSamples ) chipSample = mcmcChain[,"chip"] bobSample = mcmcChain[,"bob"] alexSample = mcmcChain[,"alex"] summary (codaSamples)
Open this post in threaded view
|

## Re: Update on using Poisson zero trick to model skewed normal distribution

Open this post in threaded view
|

## Re: Update on using Poisson zero trick to model skewed normal distribution

 Hey John, I have already tried it with a sample size of 500 (without success), but I'm happy to do it with 1,000.  I'll give it a try and see what happens.  :) And thanks for the simpler versions of the formula.   AND I'll re-visit your previous email and figure out what's going on.  ;) I guess the thing that is troubling me is that our actual data sets contain fewer than 20 points each - so if it requires 1,000 data points to be able to get to the parameters of the underlying populations, then how much confidence can we have in any results we gain from our data?  :( Thanks again, Gwen
Open this post in threaded view
|

## Re: Update on using Poisson zero trick to model skewed normal distribution

Open this post in threaded view
|

## Re: Update on using Poisson zero trick to model skewed normal distribution

 Hey John, Success!  Your version of the equation - using both functions (dnorm & phi) worked and matched the 3 parameters for a data set with 1,000 points AND a data set for 500 points.  It didn't do as well on a data set with only 50 points - as could be expected...   We have 5 sets of data, each with just shy of 20 data points, but at the moment we are particularly interested in comparing 3 of them.   We've run the one-way anova with nonhomogeneous variance JAGS code (assuming normal distribution) - but we don't think our data are distributed normally and want to run it again using this poisson zeros trick to allow a skewed normal distribution.   Now that I am confident that I've got a working formula - thank you! - I'll get to work incorporating it into the code.  And I'll see about using a hierarchical approach, as you suggested.   Thanks again for all your detailed help - you may not remember, but you gave me quite a lot of help a couple of years ago when I was first trying to get R and BUGS installed and working together on my laptop.  I don't know very many people who are so willing to give so much of their time to help strangers with newbie problems like mine...  (One of our colleagues took one of your 2 week summer courses and that's how we learned about the Bayesian approach to analyzing statistical data in the first place.)  So, thanks again!  And I'll let you know what happens.  :) Take care, Gwen
Open this post in threaded view
|

## Re: Update on using Poisson zero trick to model skewed normal distribution

Open this post in threaded view
|

## Re: Update on using Poisson zero trick to model skewed normal distribution

 In reply to this post by Gwendolyn Campbell As a learning exercise I am trying to recreate your Skew Normal fitting method, but have not been able to get reasonable fits on the randomly generated SN data. Could you look at my code below to see if my model specification is correct? I have not used the zeros trick before, but if I can get this to work, I may try expanding the model to a skew-T distribution to allow for higher kurtosis modeling. After specifying the data and zeros: ---  rsnData<-rsn(n=1000,xi=0,omega=1,alpha=0.5) ---  rsnZeros<-c(rep(0,length(rsnData))) I ran the program for 500000 samples (i tried 10 times this too, but no change...) The posterior plots of the sampled parameters are not even close to the true parameters. Any thoughts?  I would greatly appreciate some skilled advice! code: #################################################### # (with alex=location, bob=scale and chip=shape) --------------------------------------   require(rjags)     library(sn)       # rm(list=ls(all=TRUE)) #------------------------------------------------------------------------------ #old# phi[i] <- log(2.5)+log(abs(pow(bob,2)/pow(chip,2)*exp(pow(chip,2))*(exp(pow(chip,2))-1)))+log(abs(bob-chip*(x[i]-alex)))+pow((-1/chip)*log(abs(1-((1/bob)*(chip)*(x[i]-alex))))-(alex-bob/chip*(exp((pow(chip,2))/2))-1),2)/(2*pow((pow(bob,2)/pow(chip,2)*exp(pow(chip,2))*(exp((pow(chip,2)))-1)),2)) # THE MODEL. modelstring = " model { for (i in 1:50) { zeros[i] ~ dpois(phi[i]) phi[i] <- -log(2/chip)+-log(dnorm((x[i]-alex)/bob,0,1))+-log(phi(chip*(x[i]-alex)/bob)) } alex ~ dunif(-2,2) bob ~ dunif(0.01,5) chip ~ dunif(0.01,5) } " # close quote for modelstring # Write model to a file, and send to BUGS: writeLines(modelstring,con="model.txt") #------------------------------------------------------------------------------ # THE DATA. # These data were created with parameter values [n=50, location=0, scale=1, shape=0.5] # define data in workspace# rsnData<-rsn(n=1000,xi=0,omega=1,alpha=0.5) # define zeros in workspace# rsnZeros<-c(rep(0,length(rsnData))) dataList = list (   zeros = rsnZeros,   x = rsnData ) #cheated here with inits = to actual params initsList = list(   alex = c(0),   bob = c(1),   chip = c(0.5) ) #------------------------------------------------------------------------------ # RUN THE CHAINS parameters = c ("alex", "bob", "chip") adaptSteps = 1000 burnInSteps = 1000 nChains = 7 numSavedSteps = 500000 thinSteps = 1 nIter = ceiling ( (numSavedSteps * thinSteps) / nChains) jagsModel = jags.model( "model.txt" , data=dataList , inits=initsList ,                         n.chains=nChains , n.adapt=adaptSteps ) cat("Burning in MCMC chain...\n") update(jagsModel , n.iter=burnInSteps) cat("Sampling final MCMC chain...\n") codaSamples = coda.samples ( jagsModel , variable.names = parameters, n.iter=nIter, thin = thinSteps) #------------------------------------------------------------------------------ # EXAMINE THE RESULTS mcmcChain = as.matrix( codaSamples ) chipSample = mcmcChain[,"chip"] bobSample = mcmcChain[,"bob"] alexSample = mcmcChain[,"alex"] a.avg<-mean(alexSample) ; plot(density(alexSample)) b.avg<-mean(bobSample)  ; plot(density(bobSample)) c.avg<-mean(chipSample) ; plot(density(chipSample)) summary (codaSamples) # plot the random data with known parameters plot(density(rsnData), xlim=c(-10,20)) xPlot<- seq(-10,20, legnth=5000) ySn<- dsn(xPlot, a.avg, b.avg,c.avg) lines(ySn) #plot actual ySn.act<-dsn(xPlot,xi=0,omega=1,alpha=0.5) lines(ySn.act) ######################################
 John, I added the C-term as you suggested to ensure phi>0, but am still getting a poor fit  of the SN parameters.  The image below represents the JAGs results (n=5000000 samples), with the blue line representing the fit SN, the red as the 'true' distribution, and the black as the normalized density of the generated rsn data. The location parameter (alex) is close, but the other two parameter estimates do not even include the true values! 2. Quantiles for each variable:          2.5%      25%     50%     75%   97.5% alex -1.03514 -0.23274 0.16914 0.56997 1.35645 bob   2.94024  3.93393 4.40042 4.73211 4.97501 chip  0.01001  0.01006 0.01014 0.01029 0.01078 I am at a loss on what the issue could be.  Any additional thoughts? Thanks, Robert