Analysis with Programming

One of the problems often dealt in Statistics is minimization of the objective function. And contrary to the linear models, there is no analytical solution for models that are nonlinear on the parameters such as logistic regression, neural networks, and nonlinear regression models (like Michaelis-Menten model). In this situation, we have to use mathematical programming or optimization. And one popular optimization algorithm is the gradient descent , which we're going to illustrate here. To start with, let's consider a simple function with closed-form solution given by \begin{equation} f(\beta) \triangleq \beta^4 - 3\beta^3 + 2. \end{equation} We want to minimize this function with respect to $\beta$. The quick solution to this, as what calculus taught us, is to compute for the first derivative of the function, that is \begin{equation} \frac{\text{d}f(\beta)}{\text{d}\beta}=4\beta^3-9\beta^2. \end{equation} Setting this to 0 to obtain the stationary point gives us \begin{align}...

In my previous article , we talked about implementations of linear regression models in R, Python and SAS. On the theoretical sides, however, I briefly mentioned the estimation procedure for the parameter $\boldsymbol{\beta}$. So to help us understand how software does the estimation procedure, we'll look at the mathematics behind it. We will also perform the estimation manually in R and in Python, that means we're not going to use any special packages, this will help us appreciate the theory. Linear Least Squares Consider the linear regression model, \[ y_i=f_i(\mathbf{x}|\boldsymbol{\beta})+\varepsilon_i,\quad\mathbf{x}_i=\left[ \begin{array}{cccc} 1&x_{11}&\cdots&x_{1p} \end{array}\right],\quad\boldsymbol{\beta}=\left[\begin{array}{c}\beta_0\\\beta_1\\\vdots\\\beta_p\end{array}\right], \] where $y_i$ is the response or the dependent variable at the $i$th case, $i=1,\cdots, N$. The $f_i(\mathbf{x}|\boldsymbol{\beta})$ is the deterministic part of the model tha...

Consider the linear regression model, $$ y_i=f_i(\boldsymbol{x}|\boldsymbol{\beta})+\varepsilon_i, $$ where $y_i$ is the response or the dependent variable at the $i$th case, $i=1,\cdots, N$ and the predictor or the independent variable is the $\boldsymbol{x}$ term defined in the mean function $f_i(\boldsymbol{x}|\boldsymbol{\beta})$. For simplicity, consider the following simple linear regression (SLR) model, $$ y_i=\beta_0+\beta_1x_i+\varepsilon_i. $$ To obtain the (best) estimate of $\beta_0$ and $\beta_1$, we solve for the least residual sum of squares (RSS) given by, $$ S=\sum_{i=1}^{n}\varepsilon_i^2=\sum_{i=1}^{n}(y_i-\beta_0-\beta_1x_i)^2. $$ Now suppose we want to fit the model to the following data, Average Heights and Weights for American Women , where weight is the response and height is the predictor. The data is available in R by default.

A family of pdfs or pmfs $\{g(t|\theta):\theta\in\Theta\}$ for a univariate random variable $T$ with real-valued parameter $\theta$ has a monotone likelihood ratio (MLR) if, for every $\theta_2>\theta_1$, $g(t|\theta_2)/g(t|\theta_1)$ is a monotone (nonincreasing or nondecreasing) function of $t$ on $\{t:g(t|\theta_1)>0\;\text{or}\;g(t|\theta_2)>0\}$. Note that $c/0$ is defined as $\infty$ if $0< c$. Consider testing $H_0:\theta\leq \theta_0$ versus $H_1:\theta>\theta_0$. Suppose that $T$ is a sufficient statistic for $\theta$ and the family of pdfs or pmfs $\{g(t|\theta):\theta\in\Theta\}$ of $T$ has an MLR. Then for any $t_0$, the test that rejects $H_0$ if and only if $T >t_0$ is a UMP level $\alpha$ test, where $\alpha=P_{\theta_0}(T >t_0)$. Example 1 To better understand the theorem, consider a single observation, $X$, from $\mathrm{n}(\theta,1)$, and test the following hypotheses: $$ H_0:\theta\leq \theta_0\quad\mathrm{versus}\quad H_1:\theta>\theta...

More on Likelihood Ratio Test, the following problem is originally from Casella and Berger (2001), exercise 8.12. Problem For samples of size $n=1,4,16,64,100$ from a normal population with mean $\mu$ and known variance $\sigma^2$, plot the power function of the following LRTs (Likelihood Ratio Tests). Take $\alpha = .05$. $H_0:\mu\leq 0$ versus $H_1:\mu>0$ $H_0:\mu=0$ versus $H_1:\mu\neq 0$ Solution The LRT statistic is given by $$ \lambda(\mathbf{x})=\frac{\displaystyle\sup_{\mu\leq 0}\mathcal{L}(\mu|\mathbf{x})}{\displaystyle\sup_{-\infty<\mu<\infty}\mathcal{L}(\mu|\mathbf{x})}, \;\text{since }\sigma^2\text{ is known}. $$ The denominator can be expanded as follows: $$ \begin{aligned} \sup_{-\infty<\mu<\infty}\mathcal{L}(\mu|\mathbf{x})&=\sup_{-\infty<\mu<\infty}\prod_{i=1}^{n}\frac{1}{\sqrt{2\pi}\sigma}\exp\left[-\frac{(x_i-\mu)^2}{2\sigma^2}\right]\\ &=\sup_{-\infty<\mu<\infty}\frac{1}{(2\pi\sigma^2)^{1/n}}\exp\left[-\displaystyle\sum...

In this post, I would like to share a simple problem about sampling analysis. And I will demonstrate how to solve this using Python and R. The first two problems are originally from Sampling: Design and Analysis book by Sharon Lohr. Problems Let $N=6$ and $n=3$. For purposes of studying sampling distributions, assume that all population values are known. $y_1 = 98$ $y_2 = 102$ $y_3=154$ $y_4 = 133$ $y_5 = 190$ $y_6=175$ We are interested in $\bar{y}_U$, the population mean. Consider eight possible samples chosen. Sample No. Sample, $\mathcal{S}$ $P(\mathcal{S})$ 1 $\{1,3,5\}$ $1/8$ 2 $\{1,3,6\}$ $1/8$ 3 $\{1,4,5\}$ $1/8$ 4 $\{1,4,6\}$ $1/8$ 5 $\{2,3,5\}$ $1/8$ 6 $\{2,3,6\}$ $1/8$ 7 $\{2,4,5\}$ $1/8$ 8 $\{2,4,6\}$ $1/8$

As promised from my recent article , here's my tutorial on how to layout and design an infographic in R. This article will serve as a template for more infographic design that I plan to share on future posts. Hence, we will go through the following sections: Layout - mainly handles by grid package. Design - style of the elements in the layout. Texts - use extrafont package for custom fonts; Shapes (lines and point characters) - use grid , although this package has been removed from CRAN (as of February 26, 2015), the compressed file of the source code of the package is still available. But if I am not mistaken, by default this package is included in R. You might check it first before installing. Plots - several choices for plotting data in R: base plot, lattice , or ggplot2 package. The Infographic We aim to obtain the following layout and design in the final output of our code:

The Philippine government has launched Open Data Philippines (data.gov.ph) last year, January 16, 2014. Accordingly, the data.gov.ph aims to make national government data searchable, accessible, and useful, with the help of the different agencies of government, and with the participation of the public. This website consolidates the data sets of different government agencies, allowing users to find specific information from a rich and continuously growing collection of public data sets. Data.gov.ph provides information on how to access these datasets and tools, such infographics and other applications, to make the information easy to understand. Users may not only view the datasets, but also share and download them as spreadsheets and other formats, for their own use. The primary goal of data.gov.ph is to foster a citizenry empowered to make informed decisions, and to promote efficiency and transparency in government. For more, check out the video:

In imaging, we deal with multivariate data, like in array form with several spectral bands. And trying to come up with interpretation across correlations of its dimensions is very challenging, if not impossible. For example let's recall the number of spectral bands of AVIRIS data we used in the previous post . There are 152 bands, so in total there are 152$\cdot$152 = 23104 correlations of pairs of random variables. How will you be able to interpret that huge number of correlations? To engage on this, it might be better if we group these variables into two and study the relationship between these sets of variables. Such statistical procedure can be done using the canonical correlation analysis (CCA). An example of this on health sciences (from Reference 2) is variables related to exercise and health. On one hand you have variables associated with exercise, observations such as the climbing rate on a stair stepper, how fast you can run, the amount of weight lifted on bench press,...

Ever wonder what's the mathematics behind face recognition on most gadgets like digital camera and smartphones? Well for most part it has something to do with statistics. One statistical tool that is capable of doing such feature is the Principal Component Analysis (PCA). In this post, however, we will not do (sorry to disappoint you) face recognition as we reserve this for future post while I'm still doing research on it. Instead, we go through its basic concept and use it for data reduction on spectral bands of the image using R. Let's view it mathematically Consider a line $L$ in a parametric form described as a set of all vectors $k\cdot\mathbf{u}+\mathbf{v}$ parameterized by $k\in \mathbb{R}$, where $\mathbf{v}$ is a vector orthogonal to a normalized vector $\mathbf{u}$. Below is the graphical equivalent of the statement:

Authors: Arnold R. Salvacion                                                                        arsalvacion@gmail.com Data Analysis and Visualization using R (blog)                                          Al-Ahmadgaid B. Asaad (maintainer) alstated@gmail.com Agricultural Land Use Evaluation System (ALUES) is an R package that evaluates land suitability for different crop production. The package is based on the Food and Agric...

Enough with the theory we recently published, let's take a break and have fun on the application of Statistics used in Data Mining and Machine Learning, the k -Means Clustering. k-means clustering is a method of vector quantization , originally from signal processing, that is popular for cluster analysis in data mining. k -means clustering aims to partition n observations into k clusters in which each observation belongs to the cluster with the nearest mean, serving as a prototype of the cluster. (Wikipedia, Ref 1.) We will apply this method to an image, wherein we group the pixels into k different clusters. Below is the image that we are going to use, Colorful Bird From Wall321 We will utilize the following packages for input and output: jpeg - Read and write JPEG images; and, ggplot2 - An implementation of the Grammar of Graphics.

Currently, I am taking Statistics for Image Analysis on my masteral, and have been exploring this topic in R. One package that has the capability in this field is the EBImage from Bioconductor , which will be showcased in this post. Installation For those using Ubuntu, you may likely to encounter this error: It has something to do with the tiff.h C header file , but it's not that serious since mytechscribblings has an effective solution for this, do check that out. Importing Data To import a raw image, consider the following codes:

There will be upcoming meet ups for R User Group Philippines and Python Philippines (PythonPH) Community. Below are the details: R Meetup topic: R for SAS users, and planning of RUG activities   venue: 9/F Sun Life Centre, 5th Avenue corner Rizal Drive,             Bonifacio Global City, 1634, Taguig date: Thursday, June 19, 2014          7:00 pm outline: Introducing R to SAS users; common SAS functions used at PPD - c/o Mark Javellosa; group discussion on equivalent packages in R; and, Sharing of experiences of actual SAS converts. Question? Ask here .

Warning: Twitter have redesigned the interface of their developers' page, thus the screenshots below are now useless. But this has nothing to do with the codes, so you can still use it. It's been two weeks for search and rescue operations of the Malaysia Airlines Flight MH370, after it vanished from the radar on March 8, 2014. And wherever they are, we hope and pray for them. Photo from VENUS - Wall of Hope & Prayers for MH370 In this post, we are going to do text data mining on Twitter tweets containing #PrayForMH370 from March 8, to March 20, 2014 using Twitter API. First, we need to have an authentication on the Twitter API, to obtain the data. In the proceeding tutorial, the idea and codes for Twitter authentication were based from Julianhi's amazing blog , and I am going to replicate his code to save a copy of it.

To the R users based in the Philippines , there will be upcoming meetup, here are the details: topics: Markov Switching Models in R            by Ohly Santos            How to use the optim function in R            by Joe Brillantes venue: 9/F Sun Life Centre, 5th Avenue corner Rizal Drive,             Bonifacio Global City, 1634, Taguig date:    Friday March 28, 2014            7:00 PM to 9:00 PM (Philippine Standard Time) deadline for RSVP:            Sunday March 23, 2014            11:59 PM (Philippine Standard Time) important:            Please fill up the form before the deadline for RSVP. Your admi...

There is one package that I've been longing. A package that will give me the power to manipulate and do any 3D stuffs in R. I tried persp and wireframe , but I find them difficult to use especially on complicated mathematical functions, like doing parametric plots. And I am just frustrated about that, since I envy the 3D graphics of Wolfram Mathematica a lot, which I exploited for about a year, before I embrace R. However, that has come to an end after Joseph Rickert introduced the plot3D (authored by Karline Soetaert ) package in his post . And for the moment, we will be playing with the surf3D function. Here is the first one, the Mollusc Shell surface plot: With parametric equations: $$  \begin{eqnarray} x(u,v)&=&\left(1.16^v\right)(1 + \cos(u))\cos(v);\nonumber\\ y(u,v)&=&\left(-1.16^v\right)(1 + \cos(u))\sin(v);\nonumber\\ z(u,v)&=&\left(-2\times 1.16^v\right)(1 + \sin(u));\nonumber \end{eqnarray} $$where $u\in[0, 2\pi],\,v\in[-15, 6]$. B...

One great package in R is the animation made by Yihui Xie . And just for fun , we are going to explore that. Our aim is to create simple animated 2D and 3D plots. Here is the first one, 2D of course The code, It's a piece of cake right? The function we used for wrapping the plot is saveGIF , this function basically collects all the plots made and use these as frames of the GIF file. In other words, the above plot was generated/looped 100 times through the curve function, and in every iteration we increased the limits of the x axis; hence rolling all the generated plots, animates x-axis towards positive values. What about 3-dimensional? Speechless,

A friend of mine asked me to code the following in R: Generate samples of size 10 from Normal distribution with $\mu$ = 3 and $\sigma^2$ = 5; Compute the $\bar{x}$ and $\bar{x}\mp z_{\alpha/2}\displaystyle\frac{\sigma}{\sqrt{n}}$ using the 95% confidence level; Repeat the process 100 times; then Compute the percentage of the confidence intervals containing the true mean.  So here is what I got, Staying with the default values, one would obtain The output is a list of Matrix and Decision , wherein the first column of the first list ( Matrix ) is the computed $\bar{x}$; the second and third columns are the lower and upper limits of the confidence interval, respectively; and the fourth column, is an array of ones -- if true mean is contained in the interval and zeros -- true mean not contained. Now how fast it would be if I were to code this in Python?

I've been thinking that it might be better to explore the Box-Jenkins ARIMA (Autoregressive Integrated Moving-Average) three-iterative modelling on Shiny. So here is what I got, this app is intended for ARIMA(2, 2, 2) subclass family only. The app has six tabs, and these are: Historical Plot; Identification; Estimation; Diagnostic; Forecast; and Data The first tab is where the time plot of the simulated time series, the series can be simulated from different subclass family of ARIMA(2, 2, 2). The order is assigned using the controls in the side panel. The values of the parameters are set in the text field right below the plot. So for example, the ARIMA(1, 1, 1) has two text fields for AR (Autoregressive) and MA (Moving-Average) parameters as shown below. The default parameter for all models is 0.3. The next tab is Identification, this is the first stage of Box-Jenkins iterative modelling. Here, the model is identified using the correlograms, Autocorrelation Fun...