{% extends "/base.html" %} {% block content%} What is important to you? Overall Happiness GDP Social Support Life Expectancy Freedom Generosity Corruption Submit form Explore What Matters To You! Explore Select those values You see as important to Happiness. **Corruption is seen as a negative statistic and will be subtracted from the overall total. Scores are normalized for this purpose. {% endblock content %}
Select those values You see as important to Happiness. **Corruption is seen as a negative statistic and will be subtracted from the overall total. Scores are normalized for this purpose.


Setup¶
pip package installation¶
To install last stable version of the package use pip package manager. Recommended way of installing Python packages is to createVirtual Environment and install everything inside it, we encourage you to do so. Navigate to the directory where you’d like to install package and type:
pip install pyinterpolate
If there are no errors you may use package.
Problem with libspatialindex_c.so¶
There is a chance that your system doesn’t have libspatialindex_c.so file. You must install it from a terminal.
Linux installation:
sudo apt install libspatialindex-dev
MacOS installation:
brew install spatialindex
Manual installation of the newest version¶
If you are interested in the newest version of package from main or dev branch then you should install it manually. To do so download package from the repository, move to the package directory and type:
pip install ./
As with pip installation we recommend use of Virtual Environments.
I’d like to use Jupyter Notebooks, what should I do?¶
If you use conda environment you may use Jupyter Notebooks and you are not obligated to create Virtual Environment - you’ll be working in conda environment. To run pyinterpolate with Jupyter Notebooks you must create environment for it and install everything from the conda environment.
Step 1:¶
conda create -n [NAME OF YOUR ENV]
Step 2:¶
conda activate [NAME OF YOUR ENV]
Step 3:¶
conda install -c conda-forge python=3.7 pip notebook
Step 4:¶
pip install pyinterpolate
More questions?¶
Don’t be scared to ask - create request in Issues. We will do our best to simlify installation or help with it.


Mobility Layout Components
index page
converted from composite Components
panel popup and panel confirmation
contentStack
tabSet
Base tests
tabset tests in Templates
accordion
contentStack menu to be replaced 1.4
accordion content menu to be redone with stack menu 1.4
page layout
splitPane
maps


WELCOME TO GIBBERS GANG OFFICIAL WEBSITE!
Gibbers is group of music talented individuals with the below attributes/qualities ;
Welcome to Gibbers Gang official websites
Gibbers is group of music talented individuals with the below attributes/qualities ;


Phospholipases that hydrolyze the acyl group attached to the 2-position of PHOSPHOGLYCERIDES. (mesh:D054467)
| Synonym | Reference | Specificity |
| cytosolic PLA2 | Exact | |
| cytosolic phospholipase A2 | Exact | |
| cytosolic phospholipase A(2) | Exact | |
| cPLA2 | Exact | |
| phospholipase A2 | Exact |
| Identifier | Name | Relation |
| HGNC:9038 | PLA2G5 | isa |
| HGNC:9029 | PLA2G10 | isa |
| FPLX:PLA2G4 | PLA2G4 | isa |
| HGNC:9040 | PLA2G7 | isa |
| HGNC:17934 | PLA2G3 | isa |
| HGNC:9030 | PLA2G1B | isa |
| HGNC:9039 | PLA2G6 | isa |
| FPLX:PLA2G2 | PLA2G2 | isa |
| HGNC:18554 | PLA2G12A | isa |
| HGNC:18555 | PLA2G12B | isa |


Applying Expert System Technology to Code Reuse with Pyke
PyCon 2008, Chicago
| Author: | Bruce Frederiksen |
|---|---|
| Date: | Fri, 14 Mar 2008 |
| Web: | pyke.sourceforge.net |
| Copyright: | © 2008, Bruce Frederiksen |
Abstract
This paper explores a new approach to code reuse using a backward-chaining rule-based system, similar to prolog, to generate a function call graph before the functions are called. This is compared with current solutions which build the call graph as the functions are called.
This approach is introduced through an open source project called Pyke (Python Knowledge Engine).
Finally, the initial results show that the utility of this approach far exceeds expectations; leading to something more akin to automatic programming rather than adaptable libraries. A call for help is given to explore the capabilities of this approach across different domains.
The Thinking that Led to Pyke
The Need for Code Reuse
At one of my contracting jobs, they had many clients running essentially the same program, but each client needed minor code modifications. Their objective was to maximize code reuse.
What is Code Reuse?
The first question is what does "code reuse" mean? And the answer that seems most logical is function reuse. Where code modifications are required, a new function can be created incorporating those modifications.
Then the remaining task is to bring the proper collection of functions together for each client.
This gets more complicated as several versions of many functions will be produced for various clients that are all available for reuse by the next client. So it's not simply the case that there will be one standard default version of each function, and then several one-off customized versions that each only apply to a single client.
The result of this function combination exercise is a function call graph.
Example 1
Let us imagine that we start out with two functions for client1:
And then client2 comes along.
Let us first suppose that we need a new version of function A, but can reuse function B1:
This is easy in any programming language and leads naturally to the idea that the functions to reuse are the lower-level ones, which can be placed into libraries.
But now let us suppose the opposite; that we need a new version of function B, but can reuse function A1:
This is where we need help.
Current Solutions
The current solutions are all run-time solutions that trap the call from function A1 to some function B and figure out which function B to use when the call is made. For example:
- O-O Dynamic Binding
- Zope Adapters
- Generic Functions
Current Solution Limitations
These solutions are all limited for the same reason. Let's look at another example to see why.
Example 2
Real world programs have many more than two functions, but we can start to see the limitations of the current solutions by looking at a three function example.
We start with one client and three functions.
When client2 was added, it could only share function A1 and had to have a new B (B2) that needs a new function with a different call interface than C, so we'll call it D1.
Then along comes client3. This time things are looking up, because all we need is a new version of function D:
Now let's see what happens when we want to call the program for client3. We know we need to start with function A1, since there is only version of function A:
But at this point we have two choices for function B. All we know for client3 is that we're supposed to use function D2, so we're left to guess about function B. So we try the first one, function B1:
It's not until function B1 tries to call some function C that we discover a problem.
This is where the current solutions break down.
Certainly for this example, it is easy to imagine a developer telling the binding system: oh yea and client3 is going to have to use function B2 as well. But more realistic call graphs are much more complicated than this; so the developer would have to specify which functions to use going back many levels.
And then when there is a change in these upper level shared functions later on, it will affect the call graphs for many clients.
So the current solutions don't scale well.
Continuing on with our example; what we need to do at this point is back up and try the other B function:
After doing this, we discover the solution for the final call graph:
What's Needed
By looking at this example, we discover two things about how to solve this problem:
Do function selection prior to calling any of the functions.
We can't wait until one function calls another to figure out what to do, because we may change our minds!
Use a standard backward-chaining rule-based algorithm.
The process of first trying function B1, then backing up and trying function B2 is exactly the process used in backward-chaining rule-based systems like prolog. They call it backtracking.
Applying Backward-Chaining to Code Reuse
The next question is how do we use a backward-chaining system to produce function call graphs?
Let's examine, conceptually, what a set of backward-chaining rules would look like to find a solution to this problem. Then we can determine how to turn this into a function call graph.
The following diagram shows goals as dotted line boxes around the rules that prove that goal. In this example, some goals only have one rule and some have two.
We also see how rules link to other goals. For example, rule Use B1 and rule Use B2 both prove the same goal: Find B. But Use B1 links to the Find C goal, while Use B2 links to Find D.
Now we can follow how these rules would be run by the knowledge engine:
- The whole process is kicked off by asking the knowledge engine for a solution to Find A.
- There is only one rule for Find A: Use A1, so the knowledge engine tries this rule.
- Use A1 needs a solution to Find B.
- The knowledge engine tries the first rule for Find B: Use B1.
- Use B1 needs a solution to Find C.
- The knowledge engine tries the only rule for Find C: Use C1, which fails for client3!
The situation now looks like:
Continuing on:
- Since there are no other rules for Find C, the Find C goal fails.
- Which means that the Use B1 rule fails.
- So the knowledge engine tries the next rule for Find B: Use B2.
- Use B2 needs a solution for Find D.
- The knowledge engine tries the first rule for Find D: Use D1, which fails for client3.
- The knowledge engine tries the next rule for Find D: Use D2, which succeeds for client3!
- The Find D goal succeeds.
- The Find B goal succeeds.
- And the Find A goal succeeds.
When we achieve final success, we have the following situation:
What remains is to translate this into a function call graph.
It becomes obvious that we want to attach our python functions directly to the backward-chaining rules:
Pyke
Pyke KRB Syntax
How does all of this look in Pyke?
Pyke has its own language for rules, which it compiles into python source modules and then imports. This gives a performance boost by circumventing nearly all of the inference engine interpretation logic. It also makes it very easy to embed short python code snippets directly within the rules to help out with the inferencing. This keeps the inference mechanism simpler as it does not have to deal with things that are already easy in a procedural language (like arithmetic and simple list manipulation).
The Pyke rule source files are called knowledge rule bases and have a .krb suffix.
We'll continue with the previous example here.
First, let's look at the rules before we attach the python functions. Here's three of the rules:
use_B2 use find_B($client) when check_function($client, B, 2) find_D($client) use_D1 use find_D($client) when check_function($client, D, 1) use_D2 use find_D($client) when check_function($client, D, 2)
Note that Pyke uses a $ to indicate pattern variables (anonymous pattern variables start with $_).
The check_function goal checks to see what version of the indicated function should be used for this client. If this is the incorrect version, it fails. If there is no indication for this function, it succeeds to allow guessing.
Attaching Python Functions to Backward-Chaining Rules
Here are the last two rules with the python code added. The rules have the python function attached to them so that the function can be returned from the goal as an additional parameter. Because this parameter does not affect the inferencing process, it is a hidden parameter.
These examples just show one line of python code, but you may have as many lines as you want:
use_D1 use find_D($client) when check_function($client, D, 1) with print "D1" use_D2 use find_D($client) when check_function($client, D, 2) with print "D2"
Pyke calls the function call graphs plans. This terms applies to both the final top-level call graph, as well as intermediate call graphs.
Calling Subordinate Plans
Now we do the same thing to add python code to the use_B2 rule:
use_B2 use find_B($client) when check_function($client, B, 2) find_D($client) with print "B2"
We have code for the B2 function, but how does it call the plan returned from the find_D goal?
The most common way is:
use_B2 use find_B($client) when check_function($client, B, 2) find_D($client) $$() with print "B2"
In general, there may be many goals in the when clause that produce plans. Each would have an indented line of python code under it with $$ indicating the subordinate function. These indented lines are combined with the lines in the with clause to form the complete python function for this rule (with the differences in indenting levels corrected).
But in this case, this would mean that print "Dx" would be executed before print "B2", which seems backwards.
To call the subordinate plan within the with clause, there is an alternate mechanism:
use_B2 use find_B($client) when check_function($client, B, 2) find_D($client) as $d with print "B2" $d()
The as $d clause stores the plan function in pattern variable $d rather than adding a call to it to the with clause. Then you can decide in the with clause whether to call it, when to call it, how many times to call it, etc.
Note that pattern variables in general can be used within the python code. These are replaced by their final bound values (as constants) after the top-level goal has been proven. Thus, the rules can also be used to determine and set constant values within the plan functions to further customize the code. This is the reason that the code for the attached python functions is placed directly in the .krb file rather than in a separate python module.
Some Final Points about Plans
Function parameters are specified at the end of the use clause with an optional taking clause:
use_B2 use find_B($client) taking (a, b = None) ...
A completed plan appears as a normal python function.
Plans may be pickled and reused.
- If you add functools.partial to copy_reg.
You don't need to import all of Pyke to unpickle and run a plan.
- Only one small Pyke module is needed.
Other Capabilities
Pyke also supports forward-chaining rules:
fc_rule_name foreach fact_base_name.fact_name(pattern...) ... assert fact_base_name.fact_name(pattern...) ...
- Pyke runs all of the forward-chaining rules whose foreach clause succeeds prior to running any backward-chaining rules. Thus, forward-chaining rules can not call backward-chaining rules and vice versa. But backward-chaining rules can examine facts asserted by forward-chaining rules.
There are different kinds of knowledge bases:
- Fact Bases:
- simply store facts.
- Rule Bases:
- store both forward-chaining and backward-chaining rules.
- can use rule base inheritance to inherit, and build upon, the rules from
another rule base.
- But only single inheritance.
- Thus each rule base has a unique root rule base.
- All rule bases that share the same root form a rule base category.
- allow selection of which rule base(s) to use through rule base
activation.
- But only one rule base per rule base category may be active at one time.
- Extensibility. You can write your own knowledge bases. These might:
- look up facts in a database
- ask users questions
- probe hardware/software settings
- Fact Bases:
Initial Results
After writing Pyke's younger brother, it occurred to me that backward-chaining could be used to automatically figure out how to join database tables together and generate SQL statements.
And if the backward-chaining rules could see which substitution variables are needed by an HTML templating system, it could automatically generate the SQL to get these data and build the code to update the template.
It seemed that it would no longer be necessary to include anything that looks like code in the HTML templates. The graphic designers could just add simple attributes to their tags and the backward-chaining system would figure out the rest. This would mean that the programmers don't need to modify the HTML templates, and the graphic designers could maintain full ownership of the HTML.
I had a WSGI front-end that would simply assert the data passed to it as facts.
The forward-chaining rules took these starting facts, parsed the cookie information, form information, browser information, and url, determined whether the user was logged in, figured out which client the request was for, established all of this as additional facts and activated the appropriate rule base for this client.
Then the WSGI front-end simply asked for a proof of the process() goal and executed the resulting plan function which returned the final HTTP status codes and HTML document.
For a page retrieval (vs. form action), the process goal used two sub-goals:
- A format_retrieval goal that read the HTML template, and built a plan to render the template, given the needed data. This goal also returned a simple descriptor of this needed data as part of its inferencing.
- A retrieve goal then took that descriptor of the needed data, built the necessary SQL statements, and cooked them into a plan to execute those statements and return the needed data as a simple dictionary.
Then the two sub plans were combined in the reverse order, to first retrieve the data and then populate the template, for the final plan that went back to the WSGI front-end.
The Pyke examples/sqlgen and examples/web_framework are simplified examples that you can look at.
Now, as it turned out, the company had been running without a president for quite awhile, and had finally hired a new president.
So just as I finished the SQL generation logic to handle unique data (vs. multi-row data) and was preparing to show some demonstrations; our new president, coming from a java background and apparently never having heard of python, decided to cancel the project.
End of contract!
Code Reuse through Automatic Programming
The fundamental lesson learned was that this technique ends up being far more capable than what I had first imagined.
More than producing adaptable libraries capable of using B1 or B2 at some point in their call graphs, this approach leads to something more akin to the back-end of a compiler -- except that the compiler front-end does not target a textual language that needs to be written and parsed; but is rather a simple observer of already known facts:
Show me your schema, and I'll build your SQL statements.
Show me your HTML templates, and I'll build the code to populate them for you.
This seems to change the whole concept of code reuse; elevating it from the realm of static libraries, to the realm of dynamic automatic programming.
Going Forward
Thinking that others might find this useful, I've re-implemented the underlying knowledge engine from scratch, with numerous improvements gained from the experience of the first attempt, and made it open source.
With the backward-chaining rule base system, many applications are possible:
- Complicated decision making applications.
- Compiler back-ends.
- The .krb compiler uses Pyke.
- Automatic SQL statement generation.
- Automatic HTML generation/template processing.
- The control module for a web framework tool.
- Incorporate new custom functions into a large set of standard functions, which may change the selection or configuration of standard functions in other parts of the program.
- Automatically re-distribute the modules of a system over different programs and computers to meet a wide range of performance and capacity goals.
- Diagnosis systems.
- E.g., Automated customer service systems.
- Program or library customization for specific uses.
- Instantiate, configure, and interconnect networks of objects to meet a specific need or situation.
Up to this point, I've been flying solo. For this project to move forward to fully explore its capabilities, I'm going to need help!
I'd like to see several early adopters run with this and try it out in different domains. Pyke is in alpha status now and is ready to start to lean on.
This paper is licensed under a Creative Commons Attribution 3.0 Unported License.


| Opcode/Instruction | Op/En | 64-Bit Mode | Compat/Leg Mode | Description |
| 66 0F 2A /r CVTPI2PD xmm, mm/m64* | RM | Valid | Valid | Convert two packed signed doubleword integers from mm/mem64 to two packed double-precision floating-point values in xmm. |
Notes: *Operation is different for different operand sets; see the Description section.
| Op/En | Operand 1 | Operand 2 | Operand 3 | Operand 4 |
| RM | ModRM:reg (w) | ModRM:r/m (r) | NA | NA |
Converts two packed signed doubleword integers in the source operand (second operand) to two packed doubleprecision floating-point values in the destination operand (first operand).
The source operand can be an MMX technology register or a 64-bit memory location. The destination operand is an XMM register. In addition, depending on the operand configuration:
In 64-bit mode, use of the REX.R prefix permits this instruction to access additional registers (XMM8-XMM15).
DEST[63:0] ← Convert_Integer_To_Double_Precision_Floating_Point(SRC[31:0]); DEST[127:64] ← Convert_Integer_To_Double_Precision_Floating_Point(SRC[63:32]);
| CVTPI2PD: | __m128d _mm_cvtpi32_pd(__m64 a) |
None.
See Table 22-6, “Exception Conditions for Legacy SIMD/MMX Instructions with XMM and without FP Exception,” in the Intel® 64 and IA-32 Architectures Software Developer’s Manual, Volume 3B.


On 9 December 2014, the Ministry of Health of Brazil reported a case of West Nile Virus (WNV) in the state of Piauí (PI). This is the first detection of a human case of WNV infection in Brazil.
The case has been under investigation since August 2014, when the patient had the first symptoms. The infection was laboratory-confirmed on 28 November 2014. The case, who was admitted to hospital in Teresina (PI), has been discharged and will undergo rehabilitation and physical therapy to recover.
Four other people presented symptoms; however, laboratory tests ruled out WNV infection. In addition to those with symptoms, tests were performed on 18 other individuals in the area and all results were negative.


procedure LoadFromResourceID(Instance: THandle; ResID: Integer);
Loads an image from a resource. Instance is the handle of the module that contains the resource. ResID is the resource ID for the bitmap
[See also]: LoadFromFile, LoadFromStream, LoadFromResourceName


string geoip_time_zone_by_country_and_region ( string $country_code [, string $region_code ] )
The geoip_time_zone_by_country_and_region() function will return the time zone corresponding to a country and region code combo. In the United States, the region code corresponds to the two-letter abbreviation of each state. In Canada, the region code corresponds to the two-letter province or territory code as attributed by Canada Post. For the rest of the world, GeoIP uses FIPS 10-4 codes to represent regions. You can check » http://www.maxmind.com/app/fips10_4 for a detailed list of FIPS 10-4 codes. This function is always available if using GeoIP Library version 1.4.1 or newer. The data is taken directly from the GeoIP Library and not from any database.
Returns the time zone on success, or FALSE if the country and region code combo cannot be found.


Dynamic linear model tutorial
1 Dynamic linear models
This text gives an introduction to using state space based dynamic regression analysis of time series. We are especially interested in extracting trends in climatic observations. The computer examples are given using my DLM toolbox for Matlab.
1.1 Introduction
Statistical analysis of time series data is usually faced with the fact that we have only one realization of a process whose properties we might not fully understand. In analysis of correlation structures, for example, we need to assume that some distributional properties of the process generating the variability stay unchanged in time. In linear trend analysis, we assume that there is an underlying change in the background that stays approximately constant over time. Dynamic regression with state space approach tries to avoid some of the problems. By explicitly allowing for variability in the regression coefficients we let the system properties change in time. Furthermore, the use of unobservable state variables allows direct modelling of the processes that are driving the observed variability, such as seasonality or external forcing, and we can explicitly allow for some modelling error.
Dynamic regression can, in very general terms, be formulated using state space representation of the of the observations and the state of the system. With a sequential definition of the processes, having conditional dependence only on the previous time step, the classical Kalman filter formulas can be used to estimate the states given the observations. When the operators involved in the definition of the system are linear we have so called dynamic linear model, DLM.
A basic model for many climatic time series consists of four elements: slowly varying background level, seasonal component, external forcing of known processes modelled by proxy variables, and stochastic noise. The noise components might contain an autoregressive correlation structure to account for correlated model residuals. However, autocorrelated errors can be problematic as they might be caused by a long memory of the process as well as by some systematic features not included in the model.
This is a short tutorial on using dynamic linear models in time series analysis. It is based on my attempts to analyse some climatic time series. For that purposes I wrote the Matlab code described here. The examples deal with univariate time series, but the code can handle multivariate data, to some extent. Here we are mostly interested in extracting the components, and using these to infer about some interesting features of the system, but not to produce predictions about the behaviour of the system in the future, although understanding the system is a first step to be able to make predictions.
The use of DLMs in time series analysis is well documented in statistical literature, see the references at the end.
1.2 State space description
General dynamic linear model can be written with a help of observation equation and model equation as
\begin{eqnarray*} \label{dlmDef1} y_{t} &=& F_{t}x_{t}+v_{t},\quad v_{t}\sim N(0,V_{t}), \\ \label{dlmDef2} x_{t} &=& G_{t}x_{t-1}+w_{t},\quad w_{t}\sim N(0,W_{t}). \end{eqnarray*}Above \(y_{t}\) are the \(p\) observations at time \(t\), with \(t=1,\dots,n\). Vector \(x_{t}\) of length \(m\) contains the unobserved states of the system that are assumed to evolve in time according to a linear system operator \(G_{t}\) (a \(m\times m\) matrix). In time series settings \(x_t\) will have elements corresponding to various components of the time series process, like trend, seasonality, etc. We observe a linear combination of the states with noise and matrix \(F_{t}\) (\(p\times m\)) is the observation operator that transforms the model states into observations. Both observation end system equations can have additive Gaussian errors with covariance matrices \(V_{t}\) and \(W_{t}\).
This formulation is quite general and flexible and it allows handling of many time series analysis problems in a single framework. Moreover, an unified computational tool can be used, i.e. a single DLM computer code can be used for various purposes. Below we give examples of different analyses. As we are dealing with linear models, we assume that the operators \(G_{t}\) and \(F_{t}\) are linear, i.e. matrices. However, they can change with the time index \(t\). The state space framework can be extended to non linear model and non Gaussian errors, and to spatial-temporal analyses as well, see, e.g., Cressie.
1.3 DLM as hierarchical statistical model
The DLM formulation can be seen as a special case of a general hierarchical statistical model with three levels: data, process and parameters (see e.g. Cressie). In terms of statistical distributions we have, first, the observations uncertainty \(p(y_{t}|x_{t},\theta)\) described by the observation equation. Second, the process uncertainty of the unknown states \(x_{t}\) and their evolution given by the process equations as \(p(x_{t}|\theta)\), and third, the uncertainty related to model parameters \(p(\theta)\). These conditional formulations allows both efficient description of the system and computational tools to estimate its components. It also combines different statistical approaches, as we can have full prior probabilities for the unknowns (Bayesian approach), estimate them by maximum likelihood and plug them back (frequentistic approach), or even fix the model parameters by expert knowledge (a trivial non statistical approach).
1.4 State and parameter estimation
Recall that \(y_{t}\) are the observations and \(x_{t}\) the hidden system states for \(t=1,\dots,n\). In addition, we have (static) model parameter \(\theta\) that contains auxiliary parameters needed in defining the model and observation errors \(W_{t}\) and \(V_{t}\) and the system matrices \(G_{t}\) and \(F_{t}\). For dynamic linear models we have efficient and well founded computational tools for all the relevant statistical distributions.
- \(p(x_{t+1}|x_{t},y_{1:t},\theta)\) by Kalman filter
- \(p(x_{t}|y_{1:t},\theta)\) by Kalman filter
- \(p(x_t|y_{1:n},\theta)\) by Kalman smoother
- \(p(x_{1:n}|y_{1:n},\theta)\) by simulation smoother
- \(p(y_{1:t}|\theta)\) by Kalman filter likelihood
- \(p(x_{1:n},\theta|y_{1:n})\) by MCMC
- \(p(x_{1:n}|y_{1:n})\) by MCMC
1.5 Recursive Kalman formulas
Here are the relevant parts of the recursive formulas for Kalman filter and smoother to estimate the marginal distributions of DLM states given the observations. We assume that the initial distributions at \(t=1\) are available. First, we perform Kalman filter forward recursion for the predicted states \(p(x_{t+1}|x_{t},y_{1:t},\theta) = N(\hat{x}_{t+1},\hat C_{t+1})\), \(t=1,2,\dots,n-1\)
\begin{align*} v_t &= y_t -F_t\hat{x}_t &\text{prediction error}\\ C^y_t &= F_t\hat C_tF^T + V_t &\text{prediction error covariance}\\ K_t &= G_{t}\hat C_tF^T C_t^{y-1} &\text{Kalman gain}\\ \hat{x}_{t+1} &= G\hat{x}_t + K_tv_t&\text{next state prior mean}\\ \hat C_{t+1} &= G_t\hat C_t(G_t-K_tF_t)^T +W&\text{next state prior covariance} \end{align*}Then, we apply Kalman smoother backward recursion to obtain the smoothed states \(p(x_t|y_{1:n},\theta) = N(\tilde{x}_t,\tilde C_t)\), for \(t=n,n-1,\dots,2,1\)
\begin{align*} L &= G_t - K_tF_t\\ r &= F_t^T C^{y-1}_t v_t + L^T r\\ N &= F_t^T C^{y-1}_t F_t + L^T N L\\ \tilde{x}_t &= \hat{x}_t + \hat C_tr &\text{smoothed state mean}\\ \tilde C_t &= \hat C_t - \hat C_tN\hat C_t &\text{smoothed state covariance} \end{align*}
Both these recursion are implemented in dlmsmo.m
Matlab code.
1.6 Simulation smoother
Kalman smoother algorithm provides the marginal distributions \(p(x_{t}|y_{1:n},\theta)\) for each \(t\). These are all Gaussian. However, for studying trends and other dynamic features in the system, we are interested in the joint distribution spanning the whole time range \(p(x_{1:n}|y_{1:n},\theta)\). Note that we are still conditioning on the unknown parameter vector \(\theta\) and will account for it later. This very high dimensional joint distribution is not easily accessible directly. As in many cases, instead of analytic expressions, it is more important to be able to draw realizations from the distribution.
The system equations provide a direct way to recursively produce realizations of both the states \(x_{1:n}\) and the observations \(y_{1:n}\). However, the generated states will be independent of the original observations. It can be shown (Durbin & Koopman, Section 4.9) that the residual process of generated vs. smoothed state is independent of the \(\tilde{x}_{1:n}\) and \(\tilde{y}_{1:n}\). This means that if we add these residuals over the smoothed state \(x_{1:n}\), we get a new realization that is conditional on \(y_{1:n}\), the original observations. So to produce \({x}^*_{1:n} \sim p(x_{1:n}|y_{1:n},\theta)\), we:
- Sample from the system equations to get \(\tilde{x}_{1:n}\) and \(\tilde{y}_{1:n}\).
- Smooth \(\tilde{y}_{1:n}\) to get \(\tilde{\tilde{x}}_{1:n}\).
- Add the residuals to the original smoothed state, \({x}^*_{1:n} = \tilde{x}_{1:n} - \tilde{\tilde{x}}_{1:n} + {x}_{1:n}\).
This simulation smoother will be used in trend studies and as a part of more general simulation algorithm that will sample from the joint posterior distribution \(p(x_{1:n},\theta|y_{1:n})\), and by marginalization argument also from \(p(x_{1:n}|y_{1:n})\) where the uncertainty in \(\theta\) has been integrated out. Simulation smoother is implemented in dlmsmo
and dlmsmosam
functions. As a note, the simulation smoother can be seen as a special case of the randomize-then-optimize (RTO) method in Bardsley et al. (2104).
1.7 Analysing trends
In general terms, trend is some change in the distributional properties, such as in the mean, of the process generating the observations. We are typically interested in slowly varying changes in the background level, i.e. in the mean process after known sources of variability, such as seasonality, have been accounted for. A common way to explore trends is to fit some kind of a smoother, such as a moving average, over the time series. Many smoothing methods do not provide statistical ways to estimate the smoothness parameters or asses the uncertainty related to the level of smoothing.
In DLM trend analysis, a slowly varying (relative to the time scale we are interested in) background level of the system is modelled as a random walk process with variance parameters that determine the time wise smoothness of the level. These variance parameters must be estimated. In an optimal case, the data provides information also on the smoothness of the trend component, but sometimes we need to use subject level prior information to decide the time scale of the changes we want to extract. DLM models offer means of providing qualitative prior information in the form of the model equations and quantitative information by prior distributions on variance parameters, see (Gamermann).
Let \(x_{\mathrm{level},t}\) be the model state that defines the background level of the process. For statistical analysis, we need to estimate the whole state, as either \(p(x_{1:n}|y_{1:n},\hat{\theta})\), where we plug in some estimates of the auxiliary parameters \(\theta\), (maximum likelihood approach) or by \(p(x_{1:n}|y_{1:n}) = \int p(x_{1:n},\theta|y_{1:n})\,d\theta\) where the uncertainty of auxiliary parameters \(\theta\) are integrated out (Bayesian approach). The latter is typically done by Markov chain Monte Carlo (MCMC) simulation.
1.8 Examples
1.8.1 Spline smoothing
A simple level and trend model can be used as a basis for many trend related studies. In terms of local level \(\mu_{t}\) and local trend \(\alpha_{t}\) it can be written as
\begin{align*} y_t &= \mu_{t} + \epsilon_\mathrm{obs}, &\epsilon_\mathrm{obs}\sim N(0,\sigma^2_\mathrm{obs}), &\text{ observations}\\ \mu_{t} &= \mu_{t-1} + \alpha_{t-1} + \epsilon_\mathrm{level}, & \epsilon_\mathrm{level}\sim N(0,\sigma^2_\mathrm{level}), &\text{ local level}\\ \alpha_{t} &=\alpha_{t-1} +\epsilon_\mathrm{trend}, &\epsilon_\mathrm{trend}\sim N(0,\sigma^2_\mathrm{trend}), &\text{ local trend}, \end{align*}which in state space representation transfers into
\begin{equation*} G=\begin{bmatrix} 1 & 1 \\ 0 & 1 \\ \end{bmatrix},\quad F= \begin{bmatrix}1 & 0 \\\end{bmatrix}, \quad x_t = \begin{bmatrix} \mu_{t}&\alpha_{t}\end{bmatrix}^T, \end{equation*}with three parameters for the error variances
\begin{equation*} \theta = \begin{bmatrix} \sigma^2_\mathrm{obs}&\sigma^2_\mathrm{level} &\sigma^2_\mathrm{trend} \end{bmatrix}^T. \end{equation*}It is interesting to note, that the local level and trend model with \(\sigma^2_\mathrm{level} = 0\) and \(\lambda = \sigma^2_\mathrm{trend}/\sigma^2_\mathrm{obs} > 0\) is equivalent to cubic spline smoothing (Durbin & Koopman, 2012) with smoothing parameter \(\lambda\).
Below, we apply this to classical Nile flow data set used in numerous time series text books. We use Matlab with DLM toolbox.
% Nile flow data y = [1120 1160 963 1210 1160 1160 813 1230 1370 1140 995 935 ... 1110 994 1020 960 1180<PHONE_NUMBER> 1100 1210 1150 1250 ... 1260 1220 1030 1100 774 840 874 694 940 833 701<PHONE_NUMBER> 1050 969 831 726 456<PHONE_NUMBER> 1100 832 764 ... 821 768 845 864 862 698 845<PHONE_NUMBER> 759 781 865 ... 845 944 984<PHONE_NUMBER> 771 676 649 846 812 742 801 ... 1040 860 874 848 890 744<PHONE_NUMBER> 918 986 797 923 ...<PHONE_NUMBER><PHONE_NUMBER> 912 746 919 718 714 740]'; t = (1871:1970)'; % time
options = struct('order',1);
dlmplotfit(dlmfit(y,122,[0 1.65],[],[],[],options),t);
Above we used known observation uncertainty standard deviation (\(\sigma_{\mathrm{obs}}=122\)) and known level and trend variability standard deviations ([0 1.65]). In later examples these values are estimated from the data. The Figure 1 below shows the data together with estimated level component \(\mu_{1:n} = x_{\mathrm{level},1:n}\) and its 95% probability envelope, formed from individual probability limits.
Figure 1: DLM smoother fitted on Nile flow time series
1.8.2 Classical linear regression
Linear regression without dynamical evolution of the regression components is a special case of DLM model, although using DLM code to fit linear regression model is not the most efficient way to achieve this. A non dynamic model is obtained by setting all model error variances to zero. Below, this is done for the Nile example, The 95% probability limits for the level coincides with classical confidence limits of the linear regression.
dlmplotfit(dlmfit(y,122,[0 0],[],[],[],options),t);
Figure 2: Linear trend with DLM
1.8.3 Seasonal component
For DLM time series models, we have two main ways to model seasonality. If ns is the number of seasons, then the full seasonal model has ns-1 components. The number is one smaller than the number of seasons, as the model typically has already a component for the mean level. Model for 4 seasons would be obtained as
[G,F] = dlmgensys(struct('fullseas',1,'ns',4,'order',-1)) G = -1 -1 -1 1 0 0 0 1 0 F = 1 0 0
Above order=-1
suppresses the default mean level component, so the function returns only the seasonal parts.
An alternative description of seasonality uses trigonometric functions and is achieved by using two model states for each harmonic component. Monthly data with annual and semiannual cycles would be defined by
[G,F] = dlmgensys(struct('trig',2,'ns',12,'order',-1)) G = 0.86603 0.5 0 0 -0.5 0.86603 0 0 0 0 0.5 0.86603 0 0 -0.86603 0.5 F = 1 0 1 0
For illustration we use monthly temperatures at Kaisaniemi, Helsinki. Fit seasonal component with one harmonic function, but allow some variability in the amplitude, with \(\sigma_\mathrm{trend} = 0.005\), \(\sigma_\mathrm{seas} = 0.4\).
load kaisaniemi.mat; % time and temp options = struct('order',1,'trig',1); w0 = [0, 0.005, 0.4, 0.4]; % mode error stds s = 2; % error in the observation operator dlm = dlmfit(temp,s,w0,[],[],[],options); dlmplotfit(dlm,time); grid; datetick; xlabel('time'); ylabel('temperature [C]') title('Monthly temperatures at Kaisaniemi, Helsinki')
Figure 3: Seasonal model on Kaisaniemi temperatures 2004-2013
See also the ozone time series example.
1.8.4 Autoregressive process
Autoregressive processes have serial dependence between the observations. For example, AR(3) process is defined by three coefficients \([\rho_{1},\rho_{2},\rho_{3}]\) and an independent innovation term \(\epsilon\) as
\begin{equation*} y_{t} = \rho_{1}y_{t-1} + \rho_{2}y_{t-2} + \rho_{3}y_{t-3} + \epsilon, \qquad \epsilon \sim N(0,\sigma^{2}_\mathrm{eps}) \end{equation*}For including autoregressive component in state space formulation we need to include state variables that "remember" the previous values. This can be achieved by suitable evolution operator \(G\). For example, AR(3) process with coefficients \([\rho_{1},\rho_{2},\rho_{3}]\), has three extra states with
\begin{equation*} G = \begin{bmatrix} \rho_1 & 1 & 0 \\ \rho_2 & 0 & 1 \\ \rho_3 & 0 & 0 \\ \end{bmatrix},\quad F = \begin{bmatrix} 1 & 0 & 0 \\ \end{bmatrix}. \end{equation*}A pure AR(3) process would then be obtained by setting the "observation error" \(\sigma_\mathrm{obs}\) to zero and the first model error component equal to the innovation variance \(\sigma^{2}_\mathrm{eps}\).
In the Matlab toolbox, you can set autoregressive coefficients by using option arphi
.
options = struct('arphi',[0.1 0.2 0.3],'order',-1); [G,F] = dlmgensys(options) G = 0.1 1 0 0.2 0 1 0.3 0 0 F = 1 0 0
This representation allows the use of DLM code to to fit AR process and estimate its parameters, see below for an example. Note also that the basic AR formulation assumes that the observations are obtained without error. In DLM formulation it is easy to have a model that has both observation error and an AR component.
1.9 Estimating the parameters
In examples above, the variance parameters defining the model error covariance matrix \(W_t\) were assumed to be known. In practice we need some estimation methodology. Basically there are three alternatives. One. Use subject level knowledge with trial and error to fix the parameters without any algorithmic tuning. Two. Use the likelihood function with a numerical optimization routine to find maximum likelihood estimates of the parameters and plug the estimates back to the equations and re-fit the DLM model. Three. Use MCMC to sample from the posterior distribution of the parameters to estimate the parameters or to integrate out their uncertainty.
To estimate the free parameters \(\theta\) in the model formulation we need the marginal likelihood function \(p(y_{1:n}|\theta)\). By the assumed Markov properties of the system, this is obtained sequentially as a by product of the Kalman filter recursion,
\begin{equation*} -2\log(p(y_{1:n}|\theta)) = \mathrm{constant}+ \sum_{t=1}^n\left[(y_t-F_t\hat{x}_{t})^TC_{t}^{y-1}(y_t-F_t\hat{x}_{t}) + \log(|C_{t}^y|)\right]. \end{equation*}On the right hand side, the parameter \(\theta\) will appear in the model predictions \(\hat{x}_{t}\) as they depend on the model formulation, on the matrix \(G_t\), as well as on the model error \(W_t\). For the same reason we need the determinant of the model prediction covariance matrix \(C_{t}^y\). A good thing is that the above likelihood can be calculated along the DLM filter and smoother recursions without much extra effort.
The dlmfit
Matlab function has code for both maximum likelihood and MCMC calculations of parameters in the diagonal of the model error matrix \(W\), autoregressive coefficients in \(G\), as well as for estimating one extra multiplicative factor for the observation error matrix \(V\). Other parameter estimation tasks can be programmed, with the help of the likelihood returned by the dlmsmo
function, also. See the AR example below.
1.9.1 Nile data again.
Next, we use the Nile data and estimate the smoothness parameter \(\sigma_\mathrm{trend}\). By default the code displays the optimization iterations performed by fminsearch
. The option winds
contains the indexes of the diagonal of \(W\) that are to be optimized. For example, you can set winds = [0 1 2 2 2 2]
to have one parameter for the second element and one common parameter for the elements 3 to 6. The input and output for the optimizer will be in square root scale, i.e. as standard deviations, but the optimization is performed over logarithms of those to ensure positivity of the results.
options = struct('order',1,'opt',1,'winds',[0 1]); dlm = dlmfit(y,122,[0 1.65],[],[],[],options); Iteration Func-count min f(x) Procedure 0 1 1102.4 1 2 1102.4 initial simplex 2 4 1102.39 expand 3 6 1102.39 contract inside 4 8 1102.39 contract inside 5 10 1102.39 contract inside 6 12 1102.39 contract inside 7 14 1102.39 contract inside 8 16 1102.39 contract inside 9 18 1102.39 contract inside 10 20 1102.39 contract inside 11 22 1102.39 contract inside Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04 wopt = 1.579 sqrt(diag(dlm.W))' ans = 0 1.579
1.9.2 AR parameters
In this example we use Kalman filter likelihood to fit parameters in the process evolution matrix \(G\). Take the residuals from the last DLM run and fit AR(1) model. For numerical stability you have to set some small non zero observation variance (this might change in the next version). This functionality in now build into dlmfit
, also. The same idea can be used to estimate any parameter in model definition using the Kalman filter likelihood.
Below, the messy looking cost function ofun
takes all the elements of previous DLM fit stored in variable dlm
except matrix G
for which the first column is replaced with given phi
. Then it runs dlmsmo
for DLM calculations and returns the -2log(likelihood), which is minimized with respect to phi
using fminsearch
.
y = dlm.resid2; % data are residual from the previous fit phi = 0.3; % initial guess for AR parameter opts = struct('order',-1,'seas',0,'arphi',phi); dlm = dlmfit(y,0.00001,[1],[],[],[],opts); % initial DLM fit % define object function for minimizing the -2*log(likelihood) ofun = @(phi)getfield(dlmsmo(dlm.y,dlm.F,dlm.V,dlm.x0,assifun(dlm.G,phi,1:size(dlm.G,1),1),dlm.W,dlm.C0,dlm.XX,0),'lik'); phiopt = fminsearch(ofun,phi) phiopt = 0.62045
2 Toolbox functions
2.1 Function dlmsmo
out = dlmsmo(y,F,V,x0,G,W,C0,X, sample);
The main workhorse for the DLM calculations. Given the definition of DLM model as input, return Kalman filter and Kalman smoother estimates of the states, an optional sample from the state and observations, and statistics related to the fit.
Input:
y |
Observations, \(n\times p\) matrix |
F |
Observation operator, \(p\times m-q\) matrix |
V |
Observation error, \(n\times p\) matrix of standard deviations |
x0 |
initial value for the state, vector of length \(m\) |
G |
Model operator, \(m-q \times m-q\) |
W |
Model error covariance, \(m \times m\) |
C0 |
Uncertainty of x0 , \(m \times m\) |
X |
Optional matrix of regression covariates, \(n\times q\) |
sample |
1 if also a sample from the state is needed |
Above \(m\) is the dimension of the state vector, \(n\) is the number of observation times, \(p\) is the dimension of the observation at each time and \(q\) is the number of regression covariates in the model.
If covariates are given in matrix X
, the observation operator and the model operator are appended with suitable entries for X
, hence their input dimensions are of size \(m-q\). However, the initial values x0
and C0
and model error covariance matrix W
should already contain the elements corresponding to columns of X
. An empty X
as input means no covariates.
Output is a structure that contains the element below, and many more (to be added).
x |
smoothed state, \(m\times n\) |
xstd |
std of x , \(n\times m\) |
xf |
predicted states |
lik |
-2log(likelihood) |
The function dlmsmo
is called by dlmfit
that contains code for constructing the system matrices and estimating for some of the model parameters. The ploting functions described below can be used to plot the model output.
2.2 Function dlmfit
out = dlmfit(y,s,w,x0,C0,X,options);
The function dlmfit
uses dlmsmo
to fit univariate DLM time series model to observations. It generated the suitable model structure according to options given. Optionally it performs optimization or MCMC sampling analysis for parameters in the diagonal of model error matrix \(W\).
The inputs are as the following.
y
- The observations to be fitted.
s
- The observation uncertainties as standard deviations, scalar or matrix of same size as
y
. w
- Square roots in the diagonal of the model error matrix \(W\), up to the last non zero element.
x0
- The initial values for the model states, if not given, the code tries to use some reasonable first guess (typically 0).
C0
- Covariance matrix for the initial model state uncertainty. Defaults to some large values in the diagonal.
X
- Covariate proxy time series. Defaults to empty matrix
[]
. options
- System matrix generation and optimization options as a structure.
This function uses dlmgensys
to generate the system matrices. So, all the options defined for it can be used. In addition, the options
structure can have the following options related to the estimation of variance parameters
opt |
if 1, perform optimization |
mcmc |
if 1, do mcmc sampling in the parameters |
nsimu |
number of mcmc simulations |
winds |
indexes in diag(W) to be estimated |
fitv |
if 1, estimate observation error, also |
As an output the function returns the result of dlmsmo
after optional estimation of model parameters. Even if no iterative optimization is performed, the function dlmsmo
is run at least twice to iterate the initial values x0
and C0
. Also, some additional elements are added to the output structure. (These are to be documented later, see the example and source.)
Currently, the optimization and MCMC simulation can both have target from the diagonal of the model error matrix \(W\) and one scale parameter for the observation error \(V\). For some models, there are predefined optimization targets. In general you should use winds
to define which indexes in diag(\(W\)) are to be estimated by maximum likelihood or sampled by MCMC.
To use winds
, you need to know the order of the model components in \(W\): trend, seasonal, AR parameters, covariates. For example if the model has trend of order 1, and two harmonic seasonal components (trig=2
) then having one parameter for the local trend variability and one common parameter for the seasonal variability would be achieved by setting options.winds = [0 1 2 2 2 2]
.
If you need to estimate something else than the diagonal of \(W\), then you can build the target function by using the lik
output from dlmsmo
. See the AR example.
2.3 Function dlmgensys
[G,F] = dlmgensys(options);
This function generates system evolution matrix and corresponding observation operator matrix for several common dynamic linear models used in time series analyses. The input variable options
is a structure having elements that define the components. The following elements are possible at the moment:
order |
order of polynomial trend |
ns |
number of seasons |
fullseas |
if 1, use full seasonal model |
trig |
number of harmonic components |
arphi |
AR coefficients |
Example
[G,F] = dlmgensys(struct('order',1,'fullseas',1,'ns',4)) G = 1 1 0 0 0 0 1 0 0 0 0 0 -1 -1 -1 0 0 1 0 0 0 0 0 1 0 F = 1 0 1 0 0
2.4 Function dlmsmosam
sample = dlmsmosam(dlm,nsam);
Sample of states from DLM fit. Input: dlm
fit object, nsam
number of samples to draw. Output: matrix of size nstate*ntime*nsam. Uses the simulation smoother, i.e. samples from joint distribution of states at all times given the full set of observations.
2.5 Plotting functions
dlmplotfit(dlm, time, yscale, yind);
Plots the observations with the estimated level component of the model.
Input: dlm
is output from dlmfit
, time
is optional x axis for the plots, yscale
optional scaling factor for y
, and yind
the column of y
to be plotted.
dlmplotdiag(dlm);
dlmplotdiag(dlm);
Plot residual diagnostics of DLM fit performed by dlmsmo
or dlmfit
.
Figure 4: Nile flow time series residual diagnostics
dlmplotmcmc(dlm);
Plot the MCMC chain used in estimation of model parameter in function dlmfit
.
3 Ozone trend analysis
To illustrate the DLM toolbox, we analyse univariate time series with local level and trend, seasonality and some proxy time series as regressors. The observations are remote sensing satellite data for stratospheric ozone. This example reproduces some of the calculations for article: M. Laine, N. Latva-Pukkila and E. Kyrölä: Analyzing time-varying trends in stratospheric ozone time series using state the space approach, Atmospheric Chemistry and Physics 14(18), 2014, doi: 10.5194/acp-14-9707-2014.
The model can be written as
\begin{equation} y_t = \mu_t + \gamma_t + \beta_{1,t} X_{1,t} + \beta_{2,t} X_{2,t} + \beta_{3,t} X_{3,t} + \epsilon_t \end{equation}Where \(\mu_{t}\) is the signal level, \(\gamma_{t}\) is the seasonal component, and \(\beta_{i,t}\) and \(X_{i,t}\) are regression coefficients and values of the proxy time series.
The DLM matrices involved here are
\begin{equation*} x_t = \begin{bmatrix} \mu_t & \alpha_t & \psi_{t,1} & \psi^*_{t,1} & \psi_{t,2} & \psi^*_{t,2} &\beta_1& \beta_2 &\beta_3 \end{bmatrix}^T, \end{equation*} \begin{equation*} G = \begin{bmatrix} 1 & 1 & 0 & 0 & 0 & 0& 0 & 0& 0 \\ 0 & 1 & 0 & 0 & 0 & 0& 0 & 0& 0 \\ 0 & 0 &\cos\left(\frac\pi6\right)&\sin\left(\frac\pi6\right)&0&0&0&0&0\\ 0 & 0 &-\sin\left(\frac\pi6\right)&\cos\left(\frac\pi6\right)&0&0&0&0&0\\ 0 & 0 &0&0&\cos\left(\frac\pi3\right)&\sin\left(\frac\pi3\right)&0&0&0\\ 0 & 0 &0&0&-\sin\left(\frac\pi3\right)&\cos\left(\frac\pi3\right)&0&0&0\\ 0 & 0 & 0 & 0 &0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 &0 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 &0 & 0 & 0 & 0 & 1\\ \end{bmatrix}, \end{equation*} \begin{equation*} F_t = \begin{bmatrix} 1 & 0 & 1 & 0 & 1 & 0& \mathrm{solar}(t) &\mathrm{qbo1}(t)&\mathrm{qbo2}(t) \\ \end{bmatrix}, \end{equation*} \begin{equation*} \theta = \begin{bmatrix} \sigma^2_\mathrm{level} &\sigma^2_\mathrm{trend}&\sigma^2_\mathrm{seas} \end{bmatrix}^T. \end{equation*}
In the toolbox you define the model as options=struct('trig',2,'order',1);
. The elements defining the proxy variable are added by dlmfit
.
Function dlmgensys
can be used to construct the matrices \(G\) and \(F\) if needed.
options = struct('trig',2,'order',1); [G, F] = dlmgensys(options) G = 1 1 0 0 0 0 0 1 0 0 0 0 0 0 0.86603 0.5 0 0 0 0 -0.5 0.86603 0 0 0 0 0 0 0.5 0.86603 0 0 0 0 -0.86603 0.5 F = 1 0 1 0 1 0
First load the example data set. It consists of stratospheric ozone observations from two satellite instruments SAGE II and GOMOS. The original satellite observations has been combined and processed to produce averaged zonal data for different altitudes. Here we use monthly averages between 45 km to 55 km and 40°N to 50°N for years 1984 to 2012. The data set given here is for demonstration purposes only.
load ozonedata.mat time = data(:,1); % time in years from 1984 to 2012 y = data(:,2); % ozone density [1/cm^3] s = data(:,3); % uncertainty standard deviation for y X = data(:,4:6); % proxies
We scale \(y\) for numerical stability.
ys = stdnan(y); yy = y./ys; ss = s./ys;
Specify prior means for some components of \(W\), the model error matrix. The w0
variable below has the square roots of the diagonal of \(W\) and it is given as an argument to dlmfit
.
ym = meannan(yy); % mean observations wtrend = abs(ym)*0.005/12; % trend std wseas = 0.01; % seasonal component std w0 = [0 wtrend wseas wseas wseas wseas];
Calculate the DLM smoother solution and do MCMC simulation over some components in the matrix \(W\). As we are doing MCMC analysis we need to specify Bayesian prior distributions for the unknown parameters. Now we have two of them, the trend standard deviation \(\sigma_\mathrm{trend}\) and the seasonal standard deviation \(\sigma_\mathrm{seas}\). By default, the estimation is done in logarithmic scale with Gaussian priors, i.e. log-normal priors on the original scale. The prior means were already given above. Prior log-normal scale parameter is given with option varcv
below. Value 1 means prior standard deviation that is approximately 100% on the original scale of the parameter.
options = struct('trig',2,'mcmc',1,'nsimu',2000); options.winds = [0 1 2 2 2 2]; options.varcv = [1 1]; dlm = dlmfit(yy,ss,w0,[],[],X,options); Sampling these parameters: name start [min,max] N(mu,s^2) w2: -5.36761 [-Inf,Inf] N(-5.36761,1^2) w3: -4.60517 [-Inf,Inf] N(-4.60517,1^2) woptv = 1x0 empty double row vector woptw = 0.0011264 0.22004 woptg = 1x0 empty double row vector
Next, plot the fitted observations and the level component.
dlmplotfit(dlm, time, ys) title(label);xlabel('time');ylabel('average O3 density [cm^{-3}]')
Figure 5: Stratospheric ozone trends
Some model diagnostic from the prediction residuals are plotted by dlmplotdiag
.
dlmplotdiag(dlm, time, ys)
Figure 6: Stratospheric ozone trends diagnostics
The residuals do not look optimal, there is some cyclic behaviour that might be caused by non-modelled seasonal term. We try again with 3 harmonics instead of two.
options = struct('trig',3,'mcmc',1,'nsimu',2000); options.winds = [0 1 2 2 2 2 2 2]; options.varcv = [1 1]; w0 = [0 wtrend wseas wseas wseas wseas wseas wseas]; dlm = dlmfit(yy,ss,w0,[],[],X,options); Sampling these parameters: name start [min,max] N(mu,s^2) w2: -5.36761 [-Inf,Inf] N(-5.36761,1^2) w3: -4.60517 [-Inf,Inf] N(-4.60517,1^2) woptv = 1x0 empty double row vector woptw = 0.0042701 0.096897 woptg = 1x0 empty double row vector
Still not happy with the result. Next thing to try is to use an autoregressive component. This usually can explain some non-modelled feature, on the cost of possibly explaining some unknown systematic effects with correlated noise.
The next version has 3 harmonics and an AR(1) component. We need some extra options for this, arphi
has the initial value for autoregressive coefficient, fitar
tells to try to estimate this value, also, and gcv
has the relative prior standard deviation (the log-normal scale parameter).
options = struct('trig',3,'mcmc',1,'nsimu',2000); options.arphi = 0.4; options.winds = [0 1 2 2 2 2 2 2 3]; options.varcv = [1 1 1]; options.fitar = 1; options.gcv = [1]; war = 0.4; w0 = [0 wtrend wseas wseas wseas wseas wseas wseas war]; dlm = dlmfit(yy,ss,w0,[],[],X,options); Sampling these parameters: name start [min,max] N(mu,s^2) w2: -5.36761 [-Inf,Inf] N(-5.36761,1^2) w3: -4.60517 [-Inf,Inf] N(-4.60517,1^2) w9: -0.916291 [-Inf,Inf] N(-0.916291,1^2) g1: 0.4 [0,1] N(0.4,1^2) woptv = 1x0 empty double row vector woptw = 0.00099304 0.013988 0.37175 woptg = 0.34958
Plot the fit and the residual diagnostics.
figure(1); clf dlmplotfit(dlm, time, ys) title(label);xlabel('time');ylabel('average O3 density [cm^{-3}]') figure(2); clf dlmplotdiag(dlm, time, ys)
Figure 7: Stratospheric ozone trends
Figure 8: Stratospheric ozone trends diagnostics
Expect for one relatively negative autocorrelation at lag 13, the results look satisfactory. Next, for trend analysis, we produce sample from the model states using dlmsmosam
. It accounts the posterior uncertainty in \(W\) and in \(G\) using the MCMC chain in dlm.chain
.
Also, to illustrate, we draw dome sample realizations of the level component over the fit plot.
nsam = 200; % number of sampled to draw from the posterior dlm_sample = dlmsmosam(dlm,nsam); figure(1) hold on for i=1:5:nsam plot(time,ys*squeeze(dlm_sample(1,:,i)),'-') end hold off
Figure 9: Simulated samples of the level component
Next, Figure 10 shows prior and posterior probability distributions for the unknown parameters, three from the diagonal of model error matrix \(W\) and one on \(G\) for autoregressive coefficient. The estimated posterior distributions look relatively well identified with respect to the priors, which mean that there is some information about the values of the parameters in the observations.
mcmcplot(dlm.chain,[],dlm.res,'denspanel',2);
Figure 10: Parameter prior and posterior from MCMC chain
Next, we sample trend statistics form the DLM sample. Calculate 10 year (+-5) running trend for each model state sample in dlm_sample
, collect the results and plot 95% probability region of our uncertainty about the decadal trend.
nyear = 10; tsamp = ys*squeeze(dlm_sample(1,:,:)); % sample of levels ysm = mean(tsamp(:)); % their mean ysf = 1/(ysm*nyear)*100; % scale factor to get % change / 10 year t10 = mean((tsamp(nyear*12+1:end,:)-tsamp(1:end-nyear*12,:))')*ysf; % mean trend s10 = std((tsamp(nyear*12+1:end,:)-tsamp(1:end-nyear*12,:))')*ysf; % std in the sample time10 = time(fix(nyear/2)*12+1:end); time10 = time10(1:length(t10)); % time axis for plot
confband(time10,t10,s10);grid; xlim([time(1),time(end)]); % match axis to other plots title('10 year trend'); ylabel('% change / year')
Figure 11: Ten year trend in ozone
As also seen from the fit plot with the observations and the estimated level, the negative trend has started to increase toward positive somewhere after the year 1995, this increase has leveled out around 2000 and the current trend is slightly positive, but not necessarily in statistical significant way, i.e. zero or negative slope (although larger than in the beginning of the period) would also be consistent with the observations and with our model.
4 References
- A.C. Harvey.
Forecasting, structural time series and the Kalman filter. Cambridge University Press, 1990.
- T.J. Durbin and S.J. Koopman.
Time Series Analysis by State Space Methods. Oxford University Press, second edition, 2012.
- G. Petris, S. Petrone, and P. Campagnoli.
Dynamic Linear Models with R. Springer, 2009.
- S. Särkkä.
Bayesian Filtering and Smoothing. Cambridge University Press, 2013.
- D. Gamerman.
Markov chain Monte Carlo – Stochastic simulation for Bayesian inference. Chapman & Hall, second edition, 2006.
- N. Cressie and C.K. Wikle.
Statistics for Spatio-Temporal Data. Wiley, 2011.


## Introduction Easily communicate among browser tabs(supports cross-origin). ### Features 1. Safely enables [cross-origin](https://developer.mozilla.org/en-US/docs/Web/HTTP/Access_control_CORS) communication among different browser tabs. Uses `PostMessage` API for communication. 2. Easy to hook custom callback at various levels. Eg: executing a custom method in Child's tab on receiving a message from Parent tab. 3. Option to provide `data-tab-opener="name"` attribute on the target link/btn(which opens up a new tab), so that it remains disable until `Child` tab initiates a handshake and is received by the `Parent` tab 4. Fully fledged API to get information regarding the tabs(Parent and Child tabs) and other communication related methods. 5. Exports in a UMD format i.e. library works everywhere. ### Browser Support Tested in Chrome, Firefox and IE(versions >= 9). ### Example ##### Create an instance / reference before using. **Opener(Parent) tab** ```js var config = { onHandshakeCallback: function () { ... }, onPollingCallback: function () { ... } } var parent = new AcrossTabs.Parent(config); ``` **New(Child tab)** ```js var config = { onReady: onReady, onInitialize: onInitialize, handshakeExpiryLimit: 4000, // msec onParentDisconnect: onParentDisconnect, onParentCommunication: onParentCommunication }; var child = new AcrossTabs.Child(); ``` That's it! The communication setup would be set up once the new tab is opened up.


Sign in using your site profile
Need a Profile? Register Now.
{* #userInformationForm *}
*EMAIL
{* traditionalSignIn_emailAddress *}
*PASSWORD
Forgot your password?
{* traditionalSignIn_password *}
{* traditionalSignIn_signInButton *}
{* /userInformationForm *}


How much space would the filing cabinets of the Stasi and the NSA use up, if the NSA would print out their 5 Zettabytes?
Made by OpenDataCity. CC-BY 3.0.
The German President, Joachim Gauck, concluded in an interview with the ZDF on 30.6.2013, that the NSA was not to be compared with the Stasi:
We know for example, that it is not like it was with the Stasi and the KGB – that there exist big filing cabinets in which all the content of our conversations are written down and nicely filed. This is not the case.
Wir wissen zum Beispiel, dass es nicht so ist wie bei der Stasi und dem KGB, dass es dicke Aktenbände gibt, in denen unsere Gesprächsinhalte alle aufgeschrieben und schön abgeheftet sind. Das ist es nicht.
This statement is completely correct. At the NSA, conversation contents are not written down nor filed - but digitally recorded, saved and can be searched and found within seconds.
In contrast to the Stasi, the NSA can count on new technologies and can therefore collect information in gigantic quantities. To get the picture, we compared the data volume in this little app:
According to a report by the NPR, the data center of the NSA in Utah will be capable of saving 5 Zettabytes (5 billion Terabyte). Assuming that a filing cabinet with 60 files (30.000 pages of paper) uses up 0,4 m², which would correspond to 120 MB of data, the printed out Utah data center would use up 17 million square kilometers. Thereby the NSA can capture 1 billion times more data than the Stasi!
This page is uses Google Analytics. The NSA may already know that you're right here. contact


this file is located at /www/index.html
Platform:
Version:
Device Name:
UUID:
X:
Y:
Z:
Watch Accelerometer
Get Location
Call 411
Send SMS
Custom Alert
Beep
Vibrate
Offline Storage: Put a couple values
Offline Storage: Retrieve Values
Exit


Final project
- Form groups of 3 or 4.
- Agree on a small project to work on.
- Use this project to each individually accomplish all the goals in the following checklist
Checklist
- Use version control.
- Collaborate using github.
- Write a well documented function.
- Write a test for your function.
- Use outputs of code directly in LaTeX.
- Open a PR.
- Review a PR.
Some potential ideas
- Choose a theorem/result from your own research or perhaps from a course you're teaching on;
- Investigate a data set;
- Write about Newton's method for finding roots of continuous functions (it doesn't really matter what you choose :))


{% extends "base.html" %} {% block content %}
Discuss subjects you are worried or excited about
{{ section.description }}
All topics →


| Developer Spotlight |
Jean-Christophe Hugly Jean-Christophe Hugly, otherwise known as J-C, is one of the initial "Project JXTA" developers from Sun. Most of J-C's work so far has been with the peerGroup, the bootstrap, and their relationship with configuration and applications. J-C's professional life before Project JXTA was spent dealing with operating systems, specifically ChorusOS, a distributed real-time OS. He joined Sun as a result of the Chorus Systems aquisition in 1997. To J-C, the most exciting thing to make of the JXTA technology is a self organizing and self healing computing substrate; a neutral primal soup that offers just enough constraints to enable a multitude of higher-level organizations to form, cohabitate, and compete for adoption. You can reach J-C at<EMAIL_ADDRESS>or on the<EMAIL_ADDRESS>mailing list.
Read the archived Developer Spotlight biographies. |


This is so embarrasing and I have no idea where that page might have gone. Are you sure you typed in the correct URL?
We seem to be having some technical difficulties. Hang tight while our Awesome Engineers fix the Bug.
You may want to return to the
Home Page.
Contact Page is available at
Contact Page.


Set items: Documentation
In
Documentation
tab we can add information about the
set item
, and a name for it.


Bank of America, Raymond James to pay US$12 million for college fund overcharges
A U.S. regulator said units of Bank of America Corp and Raymond James Financial Inc will pay roughly 12 million in restitution to customers who incurred excess fees on their investments in so-called 529 plans, which are used to encourage saving for educational expenses.
REUTERS: Bank of America Corp and Raymond James Financial Inc will pay roughly US$12 million in restitution to customers who incurred excessive fees on investments meant to help their children afford college, a U.S. regulator said on Wednesday.
The Financial Industry Regulatory Authority said Bank of America's Merrill Lynch unit will pay at least US$4 million and two Raymond James units will pay US$8.03 million, after failing to ensure that financial advisers properly took fees into account when recommending investments in so-called 529 savings plans.
Such plans are tax-advantaged municipal securities designed to encourage saving for education, and sponsored by U.S. states, state agencies and educational institutions.
FINRA said the supervisory failures caused Merrill and Raymond James customers, especially those with young children many years away from college, to incur higher fees than they should have.
According to the regulator, many customers who invested in Class C shares, which carry no upfront sales charges but higher annual fees, would have been better off investing in Class A shares that carry upfront sales charges but lower annual fees.
The failures allegedly occurred from 2011 and 2015 at Merrill, and from January 2008 to March 2017 at Raymond James.
None of the units admitted or denied wrongdoing, and FINRA said they received credit for "extraordinary cooperation" with the regulator.
(Reporting by Jonathan Stempel in New York; Editing by Chizu Nomiyama and Marguerita Choy)


What this post is about: Data cleansing in practice with R
Data analysis, in practice, consists typically of some different steps which can be subsumed as “preparing data” and “model data” (not considering communication here):
(Inspired by this)
Often, the first major part – “prepare” – is the most time consuming. This can be lamented since many analysts prefer the cool modeling aspects (since I want to show my math!). In practice, one rather has to get his (her) hands dirt…
In this post, I want to put together some kind of checklist of frequent steps in data preparation. More precisely, I would like to detail some typical steps in “cleansing” your data. Such steps include:
- identify missings
- identify outliers
- check for overall plausibility and errors (e.g, typos)
- identify highly correlated variables
- identify variables with (nearly) no variance
- identify variables with strange names or values
- check variable classes (eg. characters vs factors)
- remove/transform some variables (maybe your model does not like categorial variables)
- rename some variables or values (especially interesting if large number)
- check some overall pattern (statistical/ numerical summaries)
- center/scale variables
Don’t get lost in big projects
Before we get in some details, let’s consider some overall guidelines. I have noticed that some projects keep growing like weed, and I find myself bewildered in some jungle… The difficulties then arise not because data or models are difficult, but due to the sheer volume of the analysis.
All in a function
Put analytical steps which belong together in one function. For example, build one function for data cleansing, give as input the raw data frame and let it spit out the processed data frame after all your cleansing steps:
cleansed_data <- cleanse_data(raw_dirty_data,
step_01 = TRUE,
step_02 = TRUE,
step_03 = TRUE)
Although functions are a bit more difficult to debug, at the end of the day it is much easier. Normally or often the steps will be run many times (for different reasons), so it is much easier if all is under one roof.
That said, pragmtic programming suggests to start easy, and to refactor frequently. So better start with a simple solution that works than to have a enormous code that chokes. Get the code running, then improve on it.
Data set for practice
The OKCupid Data set (sanitized version, no names!) is quite nice. You can download it here.
In the following, I will assume that these data are available.
data(profiles, package = "okcupiddata")
data <- profiles
So, the data set is quite huge: 59946, 22 (rows/cols)
Let’s have a brief look at the data.
library(dplyr)
glimpse(data)
## Observations: 59,946
## Variables: 22
## $ age <int> 22, 35, 38, 23, 29, 29, 32, 31, 24, 37, 35, 28, 24...
## $ body_type <chr> "a little extra", "average", "thin", "thin", "athl...
## $ diet <chr> "strictly anything", "mostly other", "anything", "...
## $ drinks <chr> "socially", "often", "socially", "socially", "soci...
## $ drugs <chr> "never", "sometimes", NA, NA, "never", NA, "never"...
## $ education <chr> "working on college/university", "working on space...
## $ ethnicity <chr> "asian, white", "white", NA, "white", "asian, blac...
## $ height <int> 75, 70, 68, 71, 66, 67, 65, 65, 67, 65, 70, 72, 72...
## $ income <int> NA, 80000, NA, 20000, NA, NA, NA, NA, NA, NA, NA, ...
## $ job <chr> "transportation", "hospitality / travel", NA, "stu...
## $ last_online <dttm> 2012-06-28 20:30:00, 2012-06-29 21:41:00, 2012-06...
## $ location <chr> "south san francisco, california", "oakland, calif...
## $ offspring <chr> "doesn't have kids, but might want them", "doesn't...
## $ orientation <chr> "straight", "straight", "straight", "straight", "s...
## $ pets <chr> "likes dogs and likes cats", "likes dogs and likes...
## $ religion <chr> "agnosticism and very serious about it", "agnostic...
## $ sex <chr> "m", "m", "m", "m", "m", "m", "f", "f", "f", "m", ...
## $ sign <chr> "gemini", "cancer", "pisces but it doesn't matter"...
## $ smokes <chr> "sometimes", "no", "no", "no", "no", "no", NA, "no...
## $ speaks <chr> "english", "english (fluently), spanish (poorly), ...
## $ status <chr> "single", "single", "available", "single", "single...
## $ essay0 <chr> "about me: i would love to think that i was som...
What’s in a name?
With the code getting longer, it is easy to get confused about naming: data_v2_no_missings_collapsed
is the right data matrix to proceed, wasn’t it? Or rather data_dat_edit_noNA_v3
? Probably a helpful (though partial) solution is to prevent typing a lot. “Don’t repeat yourself” - if stuff is put inside a function, the objects will not clutter your environment, and you don’t have to deal with them all the time. Also, you will reduce code, and the number of objects if stuff is put inside functions and loops.
That raises the question who to name the data set? At least two points are worth thinking. First, the “root” name should it be the name of the project such as “OKCupid” or “nycflights13”? Or just something like “data”? Personally, I prefer “data” as it is short and everybody will know what’s going on. Second question: Should some data manipulations should be visible in the name of the object? I, personally, like this and do that frequently, eg., carat_mean
, so I will remember that the mean is in there. However, for bigger projects I have made the experience that I lose track on which is the most recent version. So better do not put data manipulations in the name of the data frame. However, what one can do is use attributes, eg.:
data_backup <- removing_missings_somehow(data) # pseudo funtion here for illustration
data <- data_backup
attr(data, "NA_status") <- "no_NA"
Checklist
Clearly, there are different ways to get to Rome; I present just one, which has proved helpful for me.
Identify missings
A key point is here to address many columns in one go, otherwise it gets laborious with large data sets. Here’s one way:
library(knitr)
data %>%
summarise_all(funs(sum(is.na(.)))) %>% kable
| age | body_type | diet | drinks | drugs | education | ethnicity | height | income | job | last_online | location | offspring | orientation | pets | religion | sex | sign | smokes | speaks | status | essay0 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 5296 | 24395 | 2985 | 14080 | 6628 | 5680 | 3 | 48442 | 8198 | 0 | 0 | 35561 | 0 | 19921 | 20226 | 0 | 11056 | 5512 | 50 | 0 | 5485 |
The function kable
prints a html table (package knitr
).
There seem to be quite a bit missings. Maybe better plot it.
library(ggplot2)
library(tidyr)
data %>%
summarise_all(funs(sum(is.na(.)))) %>%
gather %>%
ggplot(aes(x = reorder(key, value), y = value)) + geom_bar(stat = "identity") +
coord_flip() +
xlab("variable") +
ylab("Absolute number of missings")
With this large number of missings, we probably will not find an easy solution. Skipping cases will hurt, but imputating may also not be appropriate.  1. i am a workaholic. 2. i love to cook regardless of whether i am at work. 3. i love to drink and
## 3 i'm not ashamed of much, but writing public text on an online dating site makes me pleasantly uncomfortable. i'll try to be as earnest as po
## 4 i work in a library and go to school. . .
## 5 hey how's it going? currently vague on the profile i know, more to come soon. looking to meet new folks outside of my circle of friends. i'm
## 6 i'm an australian living in san francisco, but don't hold that against me. i spend most of my days trying to build cool stuff for my company
OK, that are length(cols_with_too_many_NA)
columns. We would not want to exclude them because they are too many.
Identify outliers
Obviously, that’s a story for numeric variable only. So let’s have a look at them first.
data %>%
select_if(is.numeric) %>% names
## [1] "age" "height" "income"
Histograms are a natural and easy way to spot them and to learn something about the distribution.
data %>%
select_if(is.numeric) %>%
gather %>%
ggplot(aes(x = value)) + facet_wrap(~ key, scales = "free", nrow = 3) +
geom_histogram()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 48445 rows containing non-finite values (stat_bin).
Especially income
may be problematic.
Find out more on gather
eg., here.
Box plots (or violin plots/bean plots) may also be a way:
data %>%
select_if(is.numeric) %>%
gather %>%
ggplot(aes(x = 1, y = value)) + facet_wrap(~ key, scales = "free") +
geom_violin() +
ylab("Value") +
xlab("Variable")
## Warning: Removed 48445 rows containing non-finite values (stat_ydensity).
Identify variables with unique values
Similar to outliers, for categorical variable we can look whether some values are seldom, e.g, only 0.1% of the times. What “often” or “barely” is, depends on … you!
library(purrr)
data %>%
select_if(negate(is.numeric)) %>%
select(-matches("essay")) %>%
select(-last_online) %>%
gather %>%
ggplot(aes(x = value)) + geom_bar() +
facet_wrap(~ key, scales = "free", ncol = 3)
Pooh, that takes ages to plot. And does not look too pretty. Maybe better don’t plot, since we are not after exploration, but just want to know if something is going wrong.
Maybe it is better if we do the following:
for each non-numeric variable do divide most frequent category by least frequent category
This gives an indication whether some categories are quite frequent in relation to others.
data %>%
select_if(is.character) %>%
summarise_each(funs(max(table(.)/(min(table(.)))))) %>%
arrange %>%
kable
| body_type | diet | drinks | drugs | education | ethnicity | job | location | offspring | orientation | pets | religion | sex | sign | smokes | speaks | status | essay0 |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 74 | 1507.727 | 129.7516 | 92.00976 | 2178.091 | 32831 | 37.20098 | 31064 | 360 | 18.65052 | 336.6818 | 209.5385 | 1.485632 | 43.46341 | 29.65946 | 21828 | 5569.7 | 14 |
Plausibility check
Plausibility check can includes checking orders of magnitude, looking for implausible values (negative body weight), among others. A good starter is to differentiate between numeric and non-numeric variables.
Numeric
data %>%
select_if(is.numeric) %>%
map(summary)
## $age
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 18.00 26.00 30.00 32.34 37.00 110.00
##
## $height
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.0 66.0 68.0 68.3 71.0 95.0 3
##
## $income
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 20000 20000 50000<PHONE_NUMBER>00 1000000 48442
The function map
comes from package purrr
; it maps each selected column to the function summary
.
For instance, the max. age of 110 appears somewhat high… And the min. height of 1 should rather be excluded before further operations start. For income
similar reasoning applies.
Non-numeric
Let’s do not look at all these essay
variables here, because “strange” values are not so straight forward to identify compared to more normal categorical variables (with less distinct values).
data %>%
select(-matches("essay")) %>%
select_if(is.character) %>%
mutate_all(factor) %>%
map(summary)
Output truncated, TL;DR
We need to convert to factor
because for character
variables, no nice summary
is available of the shelf.
Highly correlated variables
data %>%
select_if(is.numeric) %>%
cor
## age height income
## age 1 NA NA
## height NA 1 NA
## income NA NA 1
Constants/ near zero variance variables
Applies to numeric variables obviously.
data %>%
select_if(is.numeric) %>%
na.omit %>%
summarise_all(c("sd", "IQR"))
## age_sd height_sd income_sd age_IQR height_IQR income_IQR
## 1 9.746844 4.003725 201433.5 12 5 80000
Rename many variables
ncol_data <- ncol(data)
names_data_new <- paste("V",1:ncol_data, sep = "")
dummy <- data
names(dummy) <- names_data_new
Recode values
A quite frequent use case is to get “strange” value or variables names in order, e.g., “variable 2” or “I rather not say” (including blanks or some other non-normal stuff).
One approach is to use dplyr::recode
.
dplyr::distinct(data, body_type)
## body_type
## 1 a little extra
## 2 average
## 3 thin
## 4 athletic
## 5 fit
## 6 <NA>
## 7 skinny
## 8 curvy
## 9 full https://sebastiansauer.github.io/images/2017-02-13d
## 10 jacked
## 11 rather not say
## 12 used up
## 13 overweight
dummy <- dplyr::recode(data$body_type, `a little extra` = "1")
unique(dummy)
## [1] "1" "average" "thin" "athletic"
## [5] "fit" NA "skinny" "curvy"
## [9] "full https://sebastiansauer.github.io/images/2017-02-13d" "jacked" "rather not say" "used up"
## [13] "overweight"
A second approach is to use base::levels
for factors.
dummy <- data$body_type
dummy <- factor(dummy)
levels(dummy)
## [1] "a little extra" "athletic" "average" "curvy"
## [5] "fit" "full https://sebastiansauer.github.io/images/2017-02-13d" "jacked" "overweight"
## [9] "rather not say" "skinny" "thin" "used up"
levels(dummy) <- 1:12
levels(dummy)
## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12"
I do not use plyr:mapvalues
because plyr
can interfere with dplyr
(and dplyr
seems more helpful to me).
Center/ scale variables
As often, several approaches. One is:
data %>%
select_if(is.numeric) %>%
scale() %>%
head
## age height income
## [1,] -1.0938888 1.67836042 NA
## [2,] 0.2813681 0.42673429 -0.1211069
## [3,] 0.5987351 -0.07391616 NA
## [4,] -0.9880998 0.67705952 -0.4189719
## [5,] -0.3533659 -0.57456662 NA
## [6,] -0.3533659 -0.32424139 NA
base::scale
performs a z-transformation of the matrix (like object) given as input.
But wait; one issue is that the data frame now consists of the numeric variables only. We have to manually knit it together with the rest of the party. Not so convenient. Rather, try this:
data %>%
select(-matches("essay")) %>%
mutate_if(is.numeric, scale) %>%
glimpse
## Observations: 59,946
## Variables: 21
## $ age <dbl> -1.09388884, 0.28136809, 0.59873508, -0.98809985, ...
## $ body_type <chr> "a little extra", "average", "thin", "thin", "athl...
## $ diet <chr> "strictly anything", "mostly other", "anything", "...
## $ drinks <chr> "socially", "often", "socially", "socially", "soci...
## $ drugs <chr> "never", "sometimes", NA, NA, "never", NA, "never"...
## $ education <chr> "working on college/university", "working on space...
## $ ethnicity <chr> "asian, white", "white", NA, "white", "asian, blac...
## $ height <dbl> 1.67836042, 0.42673429, -0.07391616, 0.67705952, -...
## $ income <dbl> NA, -0.1211069, NA, -0.4189719, NA, NA, NA, NA, NA...
## $ job <chr> "transportation", "hospitality / travel", NA, "stu...
## $ last_online <dttm> 2012-06-28 20:30:00, 2012-06-29 21:41:00, 2012-06...
## $ location <chr> "south san francisco, california", "oakland, calif...
## $ offspring <chr> "doesn't have kids, but might want them", "doesn't...
## $ orientation <chr> "straight", "straight", "straight", "straight", "s...
## $ pets <chr> "likes dogs and likes cats", "likes dogs and likes...
## $ religion <chr> "agnosticism and very serious about it", "agnostic...
## $ sex <chr> "m", "m", "m", "m", "m", "m", "f", "f", "f", "m", ...
## $ sign <chr> "gemini", "cancer", "pisces but it doesn't matter"...
## $ smokes <chr> "sometimes", "no", "no", "no", "no", "no", NA, "no...
## $ speaks <chr> "english", "english (fluently), spanish (poorly), ...
## $ status <chr> "single", "single", "available", "single", "single...
If we only want to center (or something similar), we could do
data %>%
mutate_if(is.numeric, funs(. - mean(.))) %>%
head
## age body_type diet drinks drugs
## 1 -10.34029 a little extra strictly anything socially never
## 2 2.65971 average mostly other often sometimes
## 3 5.65971 thin anything socially <NA>
## 4 -9.34029 thin vegetarian socially <NA>
## 5 -3.34029 athletic <NA> socially never
## 6 -3.34029 average mostly anything socially <NA>
## education ethnicity height income
## 1 working on college/university asian, white NA NA
## 2 working on space camp white NA NA
## 3 graduated from masters program <NA> NA NA
## 4 working on college/university white NA NA
## 5 graduated from college/university asian, black, other NA NA
## 6 graduated from college/university white NA NA
## job last_online
## 1 transportation 2012-06-28 20:30:00
## 2 hospitality / travel 2012-06-29 21:41:00
## 3 <NA> 2012-06-27 09:10:00
## 4 student 2012-06-28 14:22:00
## 5 artistic / musical / writer 2012-06-27 21:26:00
## 6 computer / hardware / software 2012-06-29 19:18:00
## location offspring
## 1 south san francisco, california doesn't have kids, but might want them
## 2 oakland, california doesn't have kids, but might want them
## 3 san francisco, california <NA>
## 4 berkeley, california doesn't want kids
## 5 san francisco, california <NA>
## 6 san francisco, california doesn't have kids, but might want them
## orientation pets
## 1 straight likes dogs and likes cats
## 2 straight likes dogs and likes cats
## 3 straight has cats
## 4 straight likes cats
## 5 straight likes dogs and likes cats
## 6 straight likes cats
## religion sex
## 1 agnosticism and very serious about it m
## 2 agnosticism but not too serious about it m
## 3 <NA> m
## 4 <NA> m
## 5 <NA> m
## 6 atheism m
## sign smokes
## 1 gemini sometimes
## 2 cancer no
## 3 pisces but it doesn't matter no
## 4 pisces no
## 5 aquarius no
## 6 taurus no
## speaks status
## 1 english single
## 2 english (fluently), spanish (poorly), french (poorly) single
## 3 english, french, c++ available
## 4 english, german (poorly) single
## 5 english single
## 6 english (fluently), chinese (okay) single
## essay0
## 1 about me: i would love to think that i was some some kind of intellectual: either the dumbest smart guy, or the smartest dumb guy. can't
## 2 i am a chef: this is what that means. 1. i am a workaholic. 2. i love to cook regardless of whether i am at work. 3. i love to drink and
## 3 i'm not ashamed of much, but writing public text on an online dating site makes me pleasantly uncomfortable. i'll try to be as earnest as po
## 4 i work in a library and go to school. . .
## 5 hey how's it going? currently vague on the profile i know, more to come soon. looking to meet new folks outside of my circle of friends. i'm
## 6 i'm an australian living in san francisco, but don't hold that against me. i spend most of my days trying to build cool stuff for my company
That’s it; happy analyzing!


{{wordsC.error|translate}}
{{wordsC.text}}
Wallet seed not available.
Your Wallet Seed:
Write it down and keep it somewhere safe.
Write it down and keep it somewhere safe.
{{wordsC.mnemonicWordsJoined}}
WARNING: This seed was created with a passphrase. To recover this wallet both the mnemonic and passphrase are needed.
Once you have written your wallet seed down, it is recommended to delete it from this device.
Failed to export


Portfolio
Game Day Trivia
Log in to play Colorado Trivia
Compete with friends to get the top score on the leaderboard
Includes trivia of Denver Broncos and Colorado
Game Day Trivia Git Hub RepoIcarus
Create an account to access your Comic Character Vault
An application that will store your favorite comic book characters
Create and store blog posts
Icarus Git Hub RepoFriendFinder
Take a short survey and find a friend with your interests
An application that will match your responses with
Save previous moods and connect it to your local weather
FriendFinder Git Hub RepoMoodList
Can't decide what music to listen to?
An application that would give you a playlist based on your mood
Save previous moods and connect it to your local weather
MoodList Git Hub RepoNifty-Gifty
Love Gifs of your favorite TV shows?
Enter your tv show and get your Gifs
Nifty-Gifty Git Hub Repo
All Aboard!
Train Schedule
Enter your departing train and your final destination
The remaining input fields and it will log your train data on the page.
All Aboard! Git Hub Repo


Once upload completes, you will be sent to a new page. You may then navigate away while the rest of the ingest completes; you will receive an email to notify you when the ingest has completed.


last updated: 17 February 2016
This page has no lang attribute on the html element
This page has no title element
This page has no h1
This page has a div which is being styled to look like an H1


Search
Clear Search
Address
Use the following abbreviations - Rd for Road, Dr for Drive, Pl for Place and Ave for Avenue
Keywords
Explore
Notable Architecture
Nabokov Trail
Manhattan Project
Farmhouses
Major Award Winners
Natural Environment
Notable Women
Public Site
Village History
Sports History
Technology & Inventors
Musicians
Artists
Writers


An experiment in tech innovation and collaboration
Create an ecosystem of innovation and collaboration between tech entrepreneurs and British Columbia's public sector.
How?
BCDevExchange is an experiment; an idea that needs to be proven. We will test, validate and refine the idea with tech entrepreneurs and early adopters in the public sector.
The BCDevExchange experiment aims to make it easier for the B.C. public and tech sectors to innovate and collaborate on digital products and services. Meet the BCDevExchange early adopters who are testing this new model of innovation and collaboration. View Profiles
Create an ecosystem of innovation and collaboration between tech entrepreneurs and British Columbia's public sector.
How?
Kick-started by the Province of B.C., the BCDevExchange is in its early stages. Contact us for how your public sector organizations can get involved, and for the following guides:
Email us to find out more or sign up and join our community.
BCDevExchange is an experiment; an idea that needs to be proven. We will test, validate and refine the idea with tech entrepreneurs and early adopters in the public sector.
The Government of B.C. is the kick-starter of the BCDevExchange experiment. We will work closely with the tech industry and early adopters in the public sector to explore the idea.
Tech entrepreneurs on the hunt for new business opportunities and innovation champions in the public sector are invited to help us test the idea. Learn how you can participate.
BCDevExchange is on GitHub, openly sharing the code and concepts that drive this website. Fork a repo and share your ideas with us!
As part of the BCDevExchange experiment, the Province of British Columbia joined GitHub. Check out the Province’s open sourced code.
The BCDevExchange experiment is about innovation and collaboration. One way to innovate is to share public sector data, code and APIs with tech entrepreneurs. Below is a list of data and APIs you can use right now for your business.
Discover data under permissive licenses such as the Open Government License - British Columbia.
BC Laws is an electronic library providing free public access to the laws of British Columbia made available for use under the permissive Queen's Printer License.
The TransLink Open API provides access to data associated with the regional bus network, bus service operations and major road network allowing developers to build applications enriched with TransLink data.
Check out what's happening in our featured projects, find opportunities, and how you can contribute!
The Labs showcase and share in-progress ideas and work from BCDevExchange early adopter public sector organizations. As safe places to experiment, the Labs contain work under design and construction that may change, disappear, or be buggy.
Visit the Labs
We are keeping you informed about the BCDevExchange on DataBC’s blog. Learn about the vision for the BCDevExchange, events, and opportunities to contribute.
Check out the blog here
Join the BCDevExchange:
You are invited to signup as a member of the BCDevExchange using your GitHub account. Sign up here.
Email:
Please contact us via email to stay informed of our progress or participate in the BCDevExchange experiment.
Twitter:
Follow @Data_BC and the hashtag #BCDev on Twitter for announcements.
The BCDevExchange experiment aims to make it easier for the public sector to innovate and to collaborate when developing digital products and services.
Please contact us for more information.


?
ÁÁ?ÁM
i
1
II
?PHH
K*wwnr?i
ii
m
ÀSet
n
3
NO.
3^
DET.
BY
MILDRED
E.
MATHIAS
*0rTf,
^/ãea
.
**£*/c
PLa;
:
DE
CUBA
**
oJ^uxaaS(L
J>
y^ULU
JltUru^U
(gaiUi/n
(
<&Ç)
À0?
ffirnant/o.
@uc£l*
'_____naif
BOTANICAL
QABDEN
Ó'"01501286
dt>
j.9
£s
01501286
Legend -
Level of confidence that token is an accurately-transcribed word
extremely low
very low
low
undetermined
medium
high
very high


SHAMS ALHASHIMIAH
energy generation project under the name of Shams Al-Hashimiah, which produces electrical energy with a capacity of 5 megawatts and at a cost of 5.6 million dinars. It contributed to reducing the university’s annual energy bill from 2.5 million dinars to zero, as the university recovered the cost of the project within a year and six months. The project covers double the university's need of electricity and meets the university's requirements and expansion, including for the next 25 years. The project was inaugurated with a ceremony attended by Prince El Hassan bin Talal .
Zaatari camp (Jordan)
November 13, 2017 - The solar power plant, which is the largest solar power plant in the world in refugee camps, was inaugurated today in the Zaatari camp in Jordan, with funding from the German government through the German Development Bank (KFW). ). The plant, valued at 15 million euros, will provide renewable and clean energy to 80,000 Syrian refugees both inside the camp and to urban areas around the camp. The solar power plant, with a maximum capacity of 12.9 megawatts, will allow UNHCR to increase the hours of electricity provision to refugee shelters from 8 to 14 hours. These additional hours will alleviate the difficult living conditions for families in the camp and improve their safety and security, in addition to the possibility of storing food and allowing children to do their homework.
HOME
A consumer who used to pay about $133.31 per month for the cost of electricity consumption could save $79.9 per month as a result of installing a solar energy system for his home, thus making his monthly bill only $53.33. If the consumer follows a rationalization system, he will be able to store a surplus that is calculated at the end of the year and the cost of solar energy for homes will be calculated based on the fees of the regular meter only.
survey
rating