The source code for this document, and the main manuscript, is available at https://github.com/richarddmorey/ConfidenceIntervalsFallacy.
We presented a situation where \(N=2\) observations were distributed uniformly:
\[ \begin{eqnarray*} y_i &\stackrel{iid}{\sim}& \mbox{Uniform}(\theta-5,\theta+5),\,i=1,\ldots,N \end{eqnarray*} \] and the goal is to estimate \(\theta\), the location of the submarine hatch. Without loss of generality we denote \(x_1\) as the smaller of the two observations. In the text, we considered five 50% confidence procedures; in this section, we give the details about the sampling distribution procedure and the Bayes procedure that were omitted from the text.
Consider the sample mean, \(\bar{y} = (y_1 + y_2)/2\). As the sum of two uniform deviates, it is a well-known fact that \(\bar{y}\) will have a triangular distribution with location \(\theta\) and minimum and maximum \(\theta-5\) and \(\theta+5\), respectively. This distribution is shown in Figure 1.
The sampling distribution of the mean \(\bar{x}\) in the submarine scenario. The shaded region represents the central 50% of the area. The unshaded triangle marked ‘a’ has area .25, and the standard deviation of this sampling distribution is about 2.04.It is desired to find the width of the base of the shaded region in Figure 1 such that it has an area of .5. To do this we first find the width of the base of the unshaded triangular area marked ‘a’ in Figure 1 such that the area of the triangle is .25. The corresponding unshaded triangle on the left side will also have area .25, which means that since the figure is a density, the shaded region must have the remaining area of .5. Elementary geometry will show that the width of the base of triangle ‘a’ is \(5/\sqrt{2}\), meaning that the distance between \(\theta\) and the altitude of triangle ‘a’ is \(5 - 5/\sqrt{2}\) or about 1.46m.
We can thus say that
\[ Pr(- (5 - 5/\sqrt{2}) < \bar{y} - \theta < 5 - 5/\sqrt{2} ) = .5 \]
which implies that, in repeated sampling,
\[ Pr(\bar{y} - (5 - 5/\sqrt{2}) < \theta < \bar{y} + (5 - 5/\sqrt{2}) ) = .5 \]
which defines the sampling distribution confidence procedure. This is an example of using \(\bar{y} - \theta\) as a pivotal quantity (Casella & Berger, 2002).
We can also derive the standard deviation of the sampling distribution of \(\bar{y}\), also called the standard error. It is defined as: \[ SE(\bar{y}) = \sqrt{V(\bar{y})} = \sqrt{\int_{-5}^{5}z^2p(z)\,dz} \] where \(p(z)\) is the triangular sampling distribution in Figure 1 centered around \(\theta=0\). Solving the integral yields \[ SE(\bar{y}) = \frac{5}{\sqrt{6}} \approx 2.04. \]
The posterior distribution is proportional to the likelihood times the prior. The likelihood is \[ p(y_1,y_2\mid\theta) \propto \prod_{i=1}^2 {\cal I}(\theta-5 < y_i < \theta+5); \] where \(\cal I\) is an indicator function. Note since this is the product of two indicator functions, it can only be nonzero when both indicator functions’ conditions are met; that is, when \(y_1+5\) and \(y_2+5\) are both greater than \(\theta\), and \(y_1-5\) and \(y_2-5\) are both less than \(\theta\). If the minimum of \(y_1+5\) and \(y_2+5\) is greater than \(\theta\), then so to must be the maximum. The likelihood thus can be rewritten \[ p(x_1,x_2\mid\theta) \propto {\cal I}(x_2 - 5 <\theta< x_1+5); \] where \(x_1\) and \(x_2\) are the minimum and maximum observations, respectively. If the prior for \(\theta\) is proportional to a constant, then the posterior is \[ p(\theta\mid x_1,x_2) \propto {\cal I}(x_2 - 5 <\theta< x_1+5), \] This posterior is a uniform distribution over all {} possible values of \(\theta\) (that is, all \(\theta\) values within 5 meters of all observations), has width \[ 10 - (x_{2} - x_{1}), \] and is centered around \(\bar{x}\). Because the posterior comprises all values of \(\theta\) the data have not ruled out – and is essentially just the classical likelihood – the width of this posterior can be taken as an indicator of the precision of the estimate of \(\theta\).
The middle 50% of the likelihood can be taken as a 50% objective Bayesian credible interval. Proof that this Bayesian procedure is also a confidence procedure is trivial and can be found in Welch (1939).
The submersible example was selected in part because it is so trivial; the confidence intervals and Bayesian credible intervals can be derived with very little effort. However, for more complicated problems, credible intervals can be more challenging to derive. Thankfully, modern Bayesian software tools make estimation of credible intervals in many problems as trivial as stating the problem along with priors on the parameters.
BUGS is a special language that allows users to define a model and prior. Using a software that interprets the BUGS language, such as JAGS (Plummer, 2003) or WinBUGS (Lunn, Thomas, Best, & Spiegelhalter, 2000), the model and prior are then combined with the data. The software then outputs samples from the posterior distribution for all the parameters, which can be used to create credible intervals.
A full explanation of how to use the BUGS language is beyond the scope of this supplement. Readers can find more information about using BUGS in Ntzoufras (2009)], (Lee & Wagenmakers, 2013), and many tutorials are available on the world wide web. Here, we show how to obtain a credible interval for the submersible Example 1 using JAGS; in a later section we show how to obtain a confidence interval for Example 2, \(\omega^2\) in ANOVA designs.
We first define the model in and prior in R using the BUGS language. Notice that this is simply stating the distributions of the data points, along with a prior for \(\theta\).
BUGS_model = "
model{
y1 ~ dunif(theta - 5, theta + 5)
y2 ~ dunif(theta - 5, theta + 5)
theta ~ dnorm( theta_mean, theta_precision)
}
"
We now define a list of values that will get passed to JAGS. y1
and y2
are the data observed values from the manuscript’s Figure 1A, and the prior we choose is an informative prior for demonstration.
for_JAGS = list( y1 = -4.5,
y2 = 4.5,
theta_mean = -2.5,
theta_precision = 1/10^2 )
Since precision is the reciprocal of variance, the prior on \(\theta\) corresponds to a Normal\((\mu=-2.5, \sigma=10)\) prior. All that remains is to load JAGS and combine the model information in BUGS_model
with the data in for_JAGS
, then to obtain samples from the posterior distribution.
# Load the rjags package to interface with JAGS
require( rjags )
## Loading required package: rjags
## Loading required package: coda
## Linked to JAGS 3.4.0
## Loaded modules: basemod,bugs
# Set initial value for the sampler
initial.values = list(theta = 0)
# Combine the model with the data
compiled_model = jags.model( file = textConnection(BUGS_model),
data = for_JAGS, inits = initial.values,
quiet = TRUE )
# Sample from the posterior distribution
posterior_samples = coda.samples( model = compiled_model,
variable.names = c("theta"),
n.iter = 100000 )
We can now plot the samples we obtained using the hist
function in R.
theta_samples = posterior_samples[[ 1 ]][ , "theta" ]
hist( theta_samples, breaks = 20, freq = FALSE )
Note the resemblance to Figure 5, bottom panel, in the manuscript. We use the summary
function on the samples to obtain a point estimate as well as quantiles of the posterior distribution, which can be used to form credible intervals. The 50% central credible interval is the interval between the 25th and 75th percentile.
summary(theta_samples)
##
## Iterations = 1001:101000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 1e+05
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## -0.0011830 0.2881537 0.0009112 0.0010825
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## -0.474832 -0.250590 -0.001707 0.248218 0.474519
In the manuscript, we compare Steiger’s (2004) confidence intervals for \(\omega^2\) to Bayesian highest posterior density (HPD) credible intervals. In this section we describe how the Bayesian HPD intervals were computed.
Consider a one-way design with \(J\) groups and \(N\) observations in each group. Let \(y_{ij}\) be the \(i\)th observation in the \(j\)th group. Also suppose that \[ y_{ij} \stackrel{indep.}{\sim} \mbox{Normal}(\mu_j, \sigma^2) \] where \(\mu_j\) is the population mean of the \(j\)th group and \(\sigma^2\) is the error variance. We assume a “non-informative” prior on parameters \(\boldsymbol\mu,\sigma^2\): \[ p(\mu_1,\ldots,\mu_J,\sigma^2) \propto (\sigma^2)^{-1}. \] This prior is flat on \((\mu_1,\ldots,\mu_J, \log\sigma^2)\). In application, it would be wiser to assume an informative prior on these parameters, in particular assuming a population over the \(\mu\) parameters or even the possibility that \(\mu_1 = \ldots = \mu_J = 0\) (Rouder, Morey, Speckman, & Province, 2012). However, for this manuscript we compare against a “non-informative” prior in order to show the differences between the confidence interval and the Bayesian result with “objective” priors.
Assuming the prior above, an elementary Bayesian calculation (Gelman, Carlin, Stern, & Rubin, 2004) reveals that \[ \sigma^2\mid\boldsymbol y \sim \mbox{Inverse Gamma}(J(N-1)/2, S/2) \] where \(S\) is the error sum-of-squares from the corresponding one-way ANOVA, and \[ \mu_j\mid\sigma^2, \boldsymbol y \stackrel{indep.}{\sim} \mbox{Normal}(\bar{x}_j, \sigma^2/N) \] where \(\mu_j\) and \(\bar{x}_j\) are the true and observed means for the \(j\)th group. Following Steiger (2004) we can define \[ \alpha_j = \mu_j - \frac{1}{J}\sum_{j=1}^J\mu_j \] as the deviation from the grand mean of the \(j\)th group, and \[ \begin{eqnarray*} \lambda &=& N\sum_{j=1}^J \left(\frac{\alpha}{\sigma}\right)^2\\ \omega^2 &=& \frac{\lambda}{\lambda + NJ}. \end{eqnarray*} \]
It is now straightforward to set up an MCMC sampler for \(\omega^2\). Let \(M\) be the number of MCMC iterations desired. We first sample \(M\) samples from the marginal posterior distribution of \(\sigma^2\), then sample the group means from the conditional posterior distribution for \(\mu_1,\ldots,\mu_J\). Using these posterior samples, \(M\) posterior samples for \(\lambda\) and \(\omega^2\) can be computed.
The following R function will sample from the marginal posterior distribution of \(\omega^2\):
## Assumes that data.frame y has two columns:
## $y is the dependent variable
## $grp is the grouping variable, as a factor
Bayes.posterior.omega2
## function (y, conf.level = 0.95, iterations = 10000)
## {
## J = nlevels(y$grp)
## N = nrow(y)/J
## aov.results = summary(aov(y ~ grp, data = y))
## SSE = aov.results[[1]][2, 2]
## sig2 = 1/rgamma(iterations, J * (N - 1)/2, SSE/2)
## lambda = matrix(NA, iterations)
## group.means = tapply(y$y, y$grp, mean)
## for (m in 1:iterations) {
## mu = rnorm(J, group.means, sqrt(sig2[m]/N))
## lambda[m] = N * sum((mu - mean(mu))^2/sig2[m])
## }
## mcmc(lambda/(lambda + N * J))
## }
The Bayes.posterior.omega2
function can be used to compute the posterior and HPD for the first example in the manuscript. The fake.data.F
function, defined in the R language in the file steiger.utility.R
(available with the manuscript source code at https://github.com/richarddmorey/ConfidenceIntervalsFallacy, generates a data set with a specified \(F\) statistic.
cl = .683 ## Confidence level corresponding to standard error
J = 3 ## Number of groups
N = 10 ## observations in a group
df1 = J - 1
df2 = J * (N - 1)
## F statistic from manuscript
Fstat = 0.1748638
set.seed(1)
y = fake.data.F(Fstat, df1, df2)
## Steiger confidence interval
steigerCI = steigerCI.omega2(Fstat,df1,df2, conf.level=cl)
samples.omega2 = Bayes.posterior.omega2(y, cl, 100000)
We can compute the Bayesian HPD interval with the HPDinterval
function in the R package coda
:
library(coda)
HPDinterval( samples.omega2, prob = cl )
## lower upper
## var1 5.219606e-06 0.08299102
## attr(,"Probability")
## [1] 0.683
Although the code above can be used to quickly sample \(\omega^2\) for any one-way design, it is not particularly generalizable for typical users. We can use the BUGS language for Bayesian modeling to create credible intervals in a way that is more accessible to the general user.
BUGS_model = "
model{
for( i in 1:NJ ){ # iterate over participants
# Error model for this observation
y[ i ] ~ dnorm( mu[ group[i] ], precision )
}
for( j in 1:J ){ # iterate over groups
# prior for group mean
mu[ j ] ~ dnorm( mean_mu, precision_mu )
# group mean's standardized squared deviation
# from overall mean
mu_dev_sq[ j ] <- pow( mu[ j ] - mean( mu ), 2 ) / variance
}
# BUGS uses the inverse of the variance (precision)
# instead of the variance
precision <- 1 / variance
# prior on error variance
variance ~ dgamma( a_variance, b_variance )
# Define our quantities of interest
lambda <- N * sum( mu_dev_sq )
omega2 <- lambda / ( lambda + N * J )
}
"
In the R code below, we define all the constants and the data needed for the analysis, including the prior parameters. These prior parameters were chosen to approximate the “non-informative” prior we used in the previous analysis. As we mentioned in the manuscript, we do not generally advise the use of such non-informative priors; these values are merely chosen for demonstration. In practice, reasonable values would be chosen to inform the analysis.
for_JAGS = list( y = y$y,
group = y$grp,
N = N,
J = J,
NJ = N*J,
mean_mu = 0,
precision_mu = 1e-6,
a_variance = 1e-6,
b_variance = 1e-6 )
The following code joins the model (BUGS_model
) with the data and defined constants (for_JAGS
) and 10,000 samples from the posterior distribution, outputting the samples of omega2
, the parameter of interest.
# Load the rjags package to interface with JAGS
require( rjags )
# Combine the model with the data
compiled_model = jags.model( file = textConnection(BUGS_model),
data = for_JAGS, quiet = TRUE )
# Sample from the posterior distribution
posterior_samples = coda.samples( model = compiled_model,
variable.names = c("omega2"),
n.iter = 10000 )
The object posterior_samples
now contains all posterior samples of the parameter \(\omega^2\). We can plot their histogram:
omega2_samples = posterior_samples[[ 1 ]][ , "omega2" ]
hist( omega2_samples, breaks = 20, freq = FALSE )
Note the close similarity between Figure 4 and Figure 3. We can do whatever we like with these samples; of particular interest would be a point estimate and credible interval. For the point estimate, we might select the posterior mean; for the credible interval, we can compute a highest-density region:
# Compute the posterior mean
mean( omega2_samples )
## [1] 0.06745069
# Compute the HDR
HPDinterval( omega2_samples , prob = cl )
## lower upper
## var1 2.832317e-05 0.08108649
## attr(,"Probability")
## [1] 0.683
Further useful information about the posterior can be obtained using the summary
function.
summary( omega2_samples )
##
## Iterations = 1001:11000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 10000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## 0.0674507 0.0594478 0.0005945 0.0005739
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## 0.002015 0.021750 0.051531 0.095971 0.223484
Casella, G., & Berger, R. L. (2002). Statistical inference. Pacific Grove, CA: Duxbury.
Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. (2004). Bayesian data analysis (2nd edition). London: Chapman; Hall.
Lee, M. D., & Wagenmakers, E.-J. (2013). Bayesian modeling for cognitive science: A practical course. Cambridge University Press.
Lunn, D., Thomas, A., Best, N., & Spiegelhalter, D. (2000). WinBUGS – a Bayesian modelling framework: Concepts, structure, and extensibility. Statistics and Computing, 10, 325–337.
Ntzoufras, I. (2009). Bayesian Modeling Using WinBUGS. Hoboken, New Jersey: Wiley.
Plummer, M. (2003). JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling. In Proceedings of the 3rd international workshop on distributed statistical computing.
Rouder, J. N., Morey, R. D., Speckman, P. L., & Province, J. M. (2012). Default Bayes factors for ANOVA designs. Journal of Mathematical Psychology, 56, 356–374. Retrieved from http://dx.doi.org/10.1016/j.jmp.2012.08.001
Steiger, J. H. (2004). Beyond the \(F\) test: Effect size confidence intervals and tests of close fit in the analysis of variance and contrast analysis. Psychological Methods, 9(2), 164–182.
Welch, B. L. (1939). On confidence limits and sufficiency, with particular reference to parameters of location. The Annals of Mathematical Statistics, 10(1), 58–69.


Who Rules InfoVis?
Data by Charles Perin and Sheelagh Carpendale. Info:
http://innovis.cpsc.ucalgary.ca/Research/WhoRulesInfovis
Originally published at
"Who Rules Infovis? unwrapping the Conference Organization"
, IEEE Conference on Visualization, Poster.


Disney Announces Five New Marvel Movies
by Cameron FrewThe Infinity Saga is over, but the show must go on: five new Marvel Cinematic Universe movies have release dates.
Advertisements
After Avengers: Endgame and Spider-Man: Far From Home closed more than 10 years of big-scale storytelling, for the first time in a decade, we don’t know where the MCU is leading.
However, with the launch of new streaming platform Disney+, the world belongs to the House of Mouse more than ever.
Advertisements
Adding to the list of upcoming superhero releases, Disney has set dates for five untitled Marvel movies, as per The Hollywood Reporter.
The unknown films are set for: October 7, 2022; February 17, 2023; May 5, 2023; July 28, 2023; and November 3, 2023. All things going to plan, that means we can expect four Marvel Studio films in 2022 and 2023.
The studio previously set May 6, 2022, for Black Panther 2, along with reserving February 18, 2022, and July 29, 2022, for other untitled Marvel projects.
What could these untitled films be? Well, James Gunn’s Guardians of the Galaxy Vol. 3 and Peyton Reed’s Ant-Man 3 are in the works currently. It’s also a relatively safe bet one of those dates is reserved for Captain Marvel 2, following the success of the first film.
Advertisements
There’s also Blade, the cult-classic vampire killer, who was announced to be returning to the big screen earlier this year. Mahershala Ali will play the hero, but there’s been no confirmation of when.
Not to mention the small matter of the Fantastic Four and X-Men – following Disney’s purchase of Fox, they’re currently in development in some capacity to join the MCU, along with Ryan Reynolds’ Deadpool.
Earlier this year, Marvel Entertainment Chief Creative Officer Kevin Feige unveiled Phase Four of the MCU. The upcoming films are: Black Widow (May 1, 2020); The Eternals (November 6, 2020); Shang-Chi and the Legend of the Ten Rings (February 12, 2021); Doctor Strange and the Multiverse of Madness (May 7, 2021); and Thor: Love and Thunder (November 5, 2021).
Advertisements
Five shows will also launch on Disney+, including Falcon and the Winter Soldier (Autumn 2020), WandaVision (Spring 2021), Loki (Spring 2021), What If…? (Summer 2021), and Hawkeye (Autumn 2021), as well as Moon Knight, She-Hulk and Ms. Marvel.
After graduating from Glasgow Caledonian University with an NCTJ and BCTJ-accredited Multimedia Journalism degree, Cameron ventured into the world of print journalism at The National, while also working as a freelance film journalist on the side, becoming an accredited Rotten Tomatoes critic in the process. He’s now left his Scottish homelands and took up residence at UNILAD as a journalist.


Welcome to the DreamFactory Admin Console!
This console is designed to allow you to create and manage nearly all of the back-end services of your project. You can create new applications, manage users, assign roles, create schemas, and/or connect to external services. You are then free to build a front-end, or client, application that utilizes those services. If you're ready, let's get started!
Where Will You Develop?
Some developers prefer to edit their HTML, CSS, and JavaScript code on their local computers, then deploy the package to DreamFactory Services Platform™. In certain situations, this option can make debugging simpler, and allows use of common desktop development tools. Other prefer to edit hosted application code directly from the Admin Console. This allows you to run and test your application immediately, on-the-fly. You can always switch back and forth between these methods as needed.
Your application's back-end has been successfully created. The next section explains how to obtain the latest DreamFactory Javascript SDK. There is also a README inside the SDK package with additional applicable information.
You can return to the Admin Console at any time to explore the other available features of your DSP. This includes, but is not limited to, managing users, creating and assigning roles, creating and managing schemata, and wiring up external services!
Download the SDK!
The DreamFactory Javascript SDK is a convenient wrapper for the entire DreamFactory Services Platform™ REST API. The SDK automatically contains your REST endpoint(s) and API key so you can get coding immediately. When a new service is provisioned from the Admin Console, the Javascript SDK will automatically be rebuilt to include the new service the next time your application is loaded.
Access to the SDK is provided by the Javascript df
object. The df
object is a property added to the DOM's
window
object. To see this in action, using Chrome's Developer Tools, open the console and start typing
window.df.apis
. This will display all the methods and properties of the df
object.
Your DSP's REST API accepts (and subsequently returns) data in either JSON or XML format. The details of switching input/output formats is explained in the API / SDK tab. That tab brings up the known and/or available services with detailed documentation, example code, and a live testing system. There are also examples for iOS and Android at http://www.dreamfactory.com/developers/examples .
Use of the SDK on the local file system requires your list of allowed hosts to be open (i.e. "*"). This setting is available in the CORS Access section under the Config tab.
Congratulations!
Your application has been successfully created. When you select the Apps tab from the left sidebar, you will see your application in the list. Each DSP instance comes with example applications to help get you started as well.
You can return to the Admin Console at any time to explore the other available features of your DSP. This includes, but is not limited to, managing users, creating and assigning roles, creating and managing schemata, and wiring up external services!
While exploring the Admin Console you'll be able to change application properties, preview your application, edit any hosted files, or even delete the application. In the file editor you will see your
index.html
file. This has example code to help you get started calling various services and building your application. There is also a
README file with additional pertinent information.
Your application's back-end has been successfully created. To get started, click the appropriate SDK download button below to grab the latest SDKs. Unzip the file into a local folder to create a workspace for building your application. To get started run the index.html file in a browser, example code is included to help you get started calling the services and building your application. There is also a README file with additional pertinent information.
You can return to the Admin Console at any time to explore the other available features of your DSP. This includes, but is not limited to, managing users, creating and assigning roles, creating and managing schemata, and wiring up external services!
Download the SDK!
The DreamFactory Javascript SDK is a convenient wrapper for the entire DreamFactory Services Platform™ REST API. The SDK automatically contains your REST endpoint(s) and API key so you can get coding immediately. When a new service is provisioned from the Admin Console, the Javascript SDK will automatically be rebuilt to include the new service the next time your application is loaded.
Access to the SDK is provided by the Javascript df
object. The df
object is a property added to the DOM's
window
object. To see this in action, using Chrome's Developer Tools, open the console and start typing
window.df.apis
. This will display all the methods and properties of the df
object.
Your DSP's REST API accepts (and subsequently returns) data in either JSON or XML format. The details of switching input/output formats is explained in the API / SDK tab. That tab brings up the known and/or available services with detailed documentation, example code, and a live testing system. There are also examples for iOS and Android at http://www.dreamfactory.com/developers/examples .
Use of the SDK on the local file system requires your list of allowed hosts to be open (i.e. "*"). This setting is available in the CORS Access section under the Config tab.


The charming adventures of an analyst and his solver.
A look at trends in Hedonometer.org's happiness data using Julia, JuMP, and GLPK.
Hedonometer.org popped onto my radar a couple weeks ago. It's a nifty project, attempting to convert samples of words found in the Twitter Gardenhose feed into a time series of happiness.
While I'm not a computational social scientist, I must say the data does have a nice intuitive quality to it. There are obvious trends in happiness associated with major holidays, days of the week, and seasons. It seems like the sort of data that could be decomposed into trends based on those various components. The Hedonometer group has, of course, done extensive analyses of their own data which you can find on their papers page.
This post examines another approach. It follows the structure of Robert Vanderbei's excellent "Local Warming" project to separate out the Hedonometer averages into daily, seasonal, solar, and day-of-the-week trends. We'll be using Julia v0.2 with JuMP and GLPK for linear optimization, Gadfly for graphing, and a few other libraries. If you haven't installed Julia, first do that. Then add all the extra requirements in the Julia shell as shown below.
Pkg.add("IJulia")
Pkg.add("DataFrames")
Pkg.add("Datetime")
Pkg.add("Gadfly")
Pkg.add("JuMP")
Pkg.add("GLPKMathProgInterface")
Pkg.add("HypothesisTests")
Alternatively, if you'd rather use an IJulia notebook for reading this post, you can download it here.
Hedonometer.org introduced an API recently, but as this post was started before that happened, we'll work with the CSV data. For the time being, you can still download it here. The ambitious reader is referred to the exercise at the end of the page regarding their API.
We use the readtable function from the DataFrames package to import the data. One oddity that this function doesn't seem to deal with yet is that the header has three columns while the rest of the CSV has two. So we skip the header and add column names in later. x1 is the date column and x2 is the happiness value, in Hedonometer units.
A brief inspection shows that we are not missing any data, so we don't have to merge in a list of dates with NAs.
using DataFrames
tmpData = readtable("dailyFull.csv", header=false, skipstart=1)
head(tmpData)
6x2 DataFrame: x1 x2 [1,] "2008-09-10" 6.01006 [2,] "2008-09-11" 6.00417 [3,] "2008-09-12" 6.01726 [4,] "2008-09-13" 6.02889 [5,] "2008-09-14" 6.03155 [6,] "2008-09-15" 5.98864
Note that the dates are strings. We vectorize the date constructor and use it to build a proper data frame.
import Datetime: date
@vectorize_1arg String date
data = DataFrame(date=date(tmpData[:, :x1]), value=tmpData[:, :x2])
head(data)
6x2 DataFrame: date value [1,] 2008-09-10 6.01006 [2,] 2008-09-11 6.00417 [3,] 2008-09-12 6.01726 [4,] 2008-09-13 6.02889 [5,] 2008-09-14 6.03155 [6,] 2008-09-15 5.98864
Let's start by making a plot of the raw data. This should look similar to the graph on their website, only not as pretty.
using Gadfly
set_default_plot_size(24cm, 12cm)
colorGradient = Scale.lab_gradient(color("crimson"), color("blue"), color("greenyellow"))
plot(
data,
x = "date",
y = "value",
color = "value",
Guide.xlabel("Date"),
Guide.ylabel("Hedonomoeter.org Units"),
Guide.title("Average Happiness on Twitter"),
Scale.ContinuousColorScale(colorGradient)
)
The data looks right, so we're off to a good start. Now we have to think about what sort of components we believe are important factors to this index. We'll start with the same ones as in the Vanderbei model:
We'll add to this weekly trends, as zooming into the data shows we tend to be happier on the weekends than on work days. In the next section we'll build a model to separate out the effects of these trends on the Hedonometer.org index.
Vanderbei's model analyzes daily temperature data for a particular location using least absolute deviations (LAD). This is similar to the well-known least squares approach, but while the latter penalizes the model quadratically more for bigger errors, the former does not. In mathematical notation, the least squares model takes in a known $m \times n$ matrix $A$ and $m \times 1$ vector $y$ of observed data, then searches for a vector $x$ such that $Ax = \hat{y}$ and $\sum_i \left\lVert y_i - \hat{y}_i \right\rVert_2^2$ is minimized.
The LAD model is similar in that it takes in the same data, but instead of minimizing the sum of the squared $L^2$ norms, it minimizes the sum of the $L^1$ norms. Thus we penalize our model using simply the absolute values of its errors instead of their squares. This makes the LAD model more robust, that is, less sensitive to outliers in our input data.
Using a robust model with this data set makes sense becuase it clearly contains a lot of outliers. While some of them, such as December 25th, may be recurrent, we're going to ignore that detail for now. After all, not every day is Christmas.
We formulate our model below using JuMP with GLPK as the solver. The code works by defining a set of variables called coefficientVars that will converge to optimal values for $x$. For each observation we compute a row of the $A$ matrix that has the following components:
We then add a linear variable representing the residual, or error, of the fitted model for each observation. Constraints enforce that these variables always take the absolute values of those errors.
Minimizing the sum of those residuals gives us a set of eight coefficients for the model. We return these and a function that predicts the happiness level for an offset from the first data record. (Note that the first record appears to be from Wednesday, September 10, 2008.)
using JuMP
using GLPKMathProgInterface
function buildPredictor(d)
m = Model(solver=GLPKSolverLP())
# Define a linear variable for each of our regression coefficients.
# Note that by default, JuMP variables are unrestricted in sign.
@defVar(m, coefficientVars[1:8])
# Residuals are the absolute values of the error comparing our
# observed and fitted values.
residuals = Array(Variable, nrow(d))
# This builds rows for determining fitted values. The first value is
# 1 since it is multiplied by our our trend line's offset. The other
# values correpond to the trends described above. Sinusoidal elements
# have two variables with sine and cosine terms.
function buildRowConstants(a1)
[
1, # Offset
a1, # Daily trend
cos(2pi*a1/365.25), # Seasonal variation
sin(2pi*a1/365.25), #
cos(2pi*a1/(10.66*365.25)), # Solar cycle variation
sin(2pi*a1/(10.66*365.25)), #
cos(2pi*a1/7), # Weekly variation
sin(2pi*a1/7) #
]
end
# This builds a linear expression as the dot product of a row's
# constants and the coefficient variables.
buildRowExpression(a1) = dot(buildRowConstants(a1), coefficientVars)
# Add a variable representing the absolute value of each error.
@defVar(m, residuals[1:nrow(d)] >= 0)
for a1 = 1:nrow(d)
fitted = buildRowExpression(a1)
# Linear way of representing residual >= |observed - fitted|
@addConstraint(m, residuals[a1] >= fitted - d[a1,:value])
@addConstraint(m, residuals[a1] >= d[a1,:value] - fitted)
end
# Minimize the total sum of these residuals.
@setObjective(m, Min, sum(residuals))
solve(m)
# Return the model coefficients and a function that predicts happiness
# for a given day, by index from the start of the data set.
coefficients = getValue(coefficientVars)[:]
# And we would like our model to work over vectors.
predictHappiness(a1) = dot(buildRowConstants(a1), coefficients)
return coefficients, predictHappiness
end
coefficients, predictor = buildPredictor(data)
coefficients
8-element Array{Float64,1}: 5.92575 9.88082e-5 -0.00144717 -0.000969706 0.125979 -0.0184735 -0.00210457 -0.0135585
The optimal values for $x$ are output above. The second value is the change in happiness per day. We can see from this that there does seem to be a small positive trend.
We can call our predictor function to obtain the fitted happiness level for any day number starting from September 10, 2008.
predictor(1000)
6.011059321544039
We now have a set of coefficients and a predictive model. That's nice, but we'd like to have some sense of a reasonable range on our model's coefficients. For instance, how certain are we that our daily trend is really even positive? To deal with these uncertanties, we use a method called bootstrapping.
Bootstrapping involves building fake observed data based on our fitted model and its associated errors. We then fit the model to our fake data to determine new coefficients. If we repeat this enough times, we may be able to generate decent confidence intervals around our model coefficients.
First step: compute the errors between the observed and fitted data. We'll construct a new data frame that contains everything we need to construct fake data.
# Compute fitted data corresponding to our observations and their associated errors.
fitted = DataFrame(
date = data[:, :date],
observed = data[:, :value],
fitted = map(x -> predictor(x), 1:nrow(data)),
error = map(x -> predictor(x) - data[x, :value], 1:nrow(data))
)
describe(fitted)
date Length 2120 Type Date{ISOCalendar} NAs 0 NA% 0.0% Unique 2119 observed Min 5.881438 1st Qu. 5.97705775 Median 6.0043168 Mean 6.0123386060849064 3rd Qu. 6.041977749999999 Max 6.369217 NAs 0 NA% 0.0% fitted Min 5.952580232891465 1st Qu. 5.979249825128504 Median 6.004645249686684 Mean 6.009750200106241 3rd Qu. 6.0430211895222055 Max 6.074612040265234 NAs 0 NA% 0.0% error Min -0.32526218693298237 1st Qu. -0.012411496422179535 Median 0.0 Mean -0.002588405978665409 3rd Qu. 0.011507686354999658 Max 0.13504237942082575 NAs 0 NA% 0.0%
Note that the median for our errors is exactly zero. This is a good sign.
Now we build a function that creates a fake input data set using the fitted values with randomonly selected errors. That is, for each observation, we add a randomly selected error with replacement to its corresponding fitted value. Once we've done that for every observation, we have a complete fake data set.
function buildFakeData(fitted)
errorIndexes = rand(1:nrow(fitted), nrow(fitted))
DataFrame(
date = fitted[:, :date],
value = map(x -> fitted[x, :fitted] - fitted[errorIndexes[x], :error], 1:nrow(fitted))
)
end
# Plot some fake data to see if it looks similar.
plot(
buildFakeData(fitted),
x = "date",
y = "value",
color = "value",
Guide.xlabel("Date"),
Guide.ylabel("Hedonomoeter.org Units"),
Guide.title("Average Happiness on Twitter (Fake)"),
Scale.ContinuousColorScale(colorGradient)
)
Visually, the plot of an arbitrary fake data set looks a lot like our original data, but not exactly.
Now we generate 1,000 fake data sets and run them through our optimization function above. This generates 1,000 sets of model coefficients and then computes $2\sigma$ confidence intervals around them.
The following code took just under an hour on my machine. If you're intent on running it yourself, you may want to get some coffee or a beer in the meantime.
using HypothesisTests
coeffData = [buildPredictor(buildFakeData(fitted))[1] for i = 1:1000]
confidence2Sigma = Array((Float64, Float64), length(coefficients))
for i = 1:length(coefficients)
sample::Array{Float64}
sample = [x[i] for x in coeffData]
confidence2Sigma[i] = ci(OneSampleTTest(sample))
end
confidence2Sigma
8-element Array{(Float64,Float64),1}: (5.925619370480642,5.926250874471983) (9.833341267642056e-5,9.891872644950681e-5) (-0.0014847851702785182,-0.0014095257149699567) (-0.001010722418756523,-0.0009395274541342083) (0.12556589790988862,0.12604831225962906) (-0.018511935408237066,-0.018348960520221922) (-0.002133611417832446,-0.002058047348615918) (-0.01359168539278434,-0.01351903662344336)
From the above output we can tell that we do seem to be trending happier, with a daily trend of 9.88082e-5 in Hedonometer units and a 95% confidence interval on that trend of 9.833341267642056e-5, and 9.891872644950681e-5. Cool!
Now we take a quick look at our model output. First, we plot the fitted happiness values for the same time period as the observed data. We can see that this resembles the same general trend minus the outliers. The width of the curve is due to weekly variation.
plot(
fitted,
x = "date",
y = "fitted",
Guide.xlabel("Date"),
Guide.ylabel("Hedonomoeter.org Units"),
Guide.title("Expected Happiness on Twitter"),
Geom.line
)
Now we take a look at what a typical week looks like in terms of its effect on our happiness. As September 10, 2008 was a Wednesday, we index Sunday starting at 5. The resulting graph highlights what I think we all already know about the work week.
daily(a1) = coefficients[6]*cos(2pi*a1/7) + coefficients[7]*sin(2pi*a1/7)
plot(
x = ["Sun", "Mon", "Tues", "Wed", "Thurs", "Fri", "Sat"],
y = map(daily, [5, 6, 7, 1, 2, 3, 4]),
Guide.xlabel("Day of the Week"),
Guide.ylabel("Hedonomoeter.org Units"),
Guide.title("Variation in Happiness on Twitter Due to Day of the Week"),
Geom.line,
Geom.point
)
That's it for this analysis. We've learned that, for the being at least, we seem to be trending happier. Which is pretty cool. Of course, 6 years is not a huge amount of time for this sort of data. Let's try doing the analysis again in 10 or 20.
Will anyone still be using Twitter then?...
The particularly ambitious reader may find the following exercises interesting.


Hey guys, I've been thinking about blowing the dust off of my old cd's and reinstalling NWN. Just wondering if any of the old timers were still around. I am hoping to be able to find a decent server or two to play on, and will probably run through the OC just to get my bearings back. It's amazing how there are just some games that you will always go back to. _________________ If everyone is thinking alike, then someone is failing to think.
Great thread to look at for you.  Glad to see I am not alone in coming back to nwn. _________________ If everyone is thinking alike, then someone is failing to think.


Let me introduce myself in brief: I am a hardcore open source bigot. I used to say evangelist but I realised the term is too weak.


Please make sure that the input is admissible:
- Age is a number between 0 and 100 years
- Respiratory rate is a number between 10 and 60 breaths per minute
- Lactate dehydrogenase (LDH) level is a number between 50 and 4000 units per liter (U/L)
- C-reactive protein (CRP) level is a number between 1 and 500 milligrams per liter (mg/L)
- Serum albumin is a number between 10 and 60 grams per liter (g/L)
- Serum urea is a number between 1 and 80 millimoles per liter (mmol/L)


Bob Ross
was
the man
That painting above has nothing to do with Bob Ross by the way. It's by Matisse. And I really like it.
One of my favorite Colombian musicians of all time is Joe Arroyo. Try to listen and not move. Go on...
And since we need a link to our email address, holla atcha boy.


Today...
We will learn more about medical imaging...
And...
We will learn how to use 3D Slicer!
This is lecture 2.


| OLD line(s): 31
| NEW line(s): 32
<xsd:include schemaLocation="../../CorporateIncomeTax/Common/IRS4562FY/IRS4562FY.xsd"/> |
| OLD line(s): 93
| NEW line(s): 95,97
<xsd:include schemaLocation="../../CorporateIncomeTax/Common/Dependencies/50YearADSDeductionStatement.xsd" /> <xsd:include schemaLocation="../../CorporateIncomeTax/Common/Dependencies/50YearGDSDeductionStatement.xsd" /> <xsd:include schemaLocation="../../CorporateIncomeTax/Common/Dependencies/AmortizationElectionStatement.xsd" /> |
| OLD line(s): 111
| NEW line(s): 116
<xsd:include schemaLocation="../../CorporateIncomeTax/Common/Dependencies/Section168f1PropertyExplanationStatement.xsd" /> |
| OLD line(s): 268,273
<!-- Includes for 4562 (supporting info, schedules, computations, and such) --> <xsd:include schemaLocation="../../CorporateIncomeTax/Common/IRS4562/50YearADSDeductionStatement.xsd"/> <xsd:include schemaLocation="../../CorporateIncomeTax/Common/IRS4562/50YearGDSDeductionStatement.xsd"/> <xsd:include schemaLocation="../../CorporateIncomeTax/Common/IRS4562/AmortizationElectionStatement.xsd"/> <xsd:include schemaLocation="../../CorporateIncomeTax/Common/IRS4562/Section168f1PropertyExplanationStatement.xsd"/> | NEW line(s): 272
|
| OLD line(s): 573
| NEW line(s): 573,579
<!-- Form 4562FY --> <xsd:element ref="IRS4562FY" minOccurs="0" maxOccurs="unbounded"> <xsd:annotation> <xsd:documentation>IRS Form 4562FY - Depreciation and Amortization (Including Information on Listed Property)</xsd:documentation> </xsd:annotation> </xsd:element> |
| OLD line(s): 1820
<!-- Form 4562 statement --> | NEW line(s): 1826
<!-- Form 4562 and Form 4562FY Statement --> |
| OLD line(s): 1827
<!-- Form 4562 statement --> | NEW line(s): 1833
<!-- Form 4562 and Form 4562FY Statement --> |
| OLD line(s): 1834
<!-- Form 4562 statement --> | NEW line(s): 1840
<!-- Form 4562 and Form 4562FY Statement --> |
| OLD line(s): 1841
<!-- Form 4562 statement --> | NEW line(s): 1847
<!-- Form 4562 and Form 4562FY Statement --> |


{% extends 'projects/project_base.html' %} {% load humanize static escape_extras source_extras project_permissions %} {% block title %}Project {{ project.name }} : Files{% endblock %} {% block project_header_right %}
{% endblock %} {% block project_content %}Are you sure you want to unlink [[ SOURCE_TYPE_NAME_LOOKUP[unlinkSourceType] ]] [[ unlinkSourcePath ]] from this project?
All files from this source will be removed from this project, but will still be available in the [[ SOURCE_TYPE_NAME_LOOKUP[unlinkSourceType] ]] repository.


Dump all your CSS/SCSS files in the 'css' folder. Call them all with a single call to 'application.css'
Dump all your JS/Coffee files in the 'js' folder. Call them all with a single call to 'application.js'
Put your layout in the 'layouts' folder. Call 'render_page' anywhere in it & Rush will bring the correct page.
No more repeating HTML ever. Put it in a partial and call it from anywhere with render_partial.


Log in.