text:
## Introduction

Welcome to “The Annotated GPT-2”.

One of the most brilliant and well-explained articles I have ever read is The Annotated Transformer. It introduced Attention like no other post ever written. The simple idea was to present an “annotated” version of the paper Attention is all you need along with code.

Something I have come to realize with my little experience in Machine Learning, when you write things in code, the implementation and the secrets become clearer. It is not magic anymore.

There is nothing magic about magic. The magician merely understands something simple which doesn’t appear to be simple or natural to the untrained audience. Once you learn how to hold a card while making your hand look empty, you only need practice before you, too, can “do magic.”

– Jeffrey Friedl in the book Mastering Regular Expressions

The GPT-2 might seem like magic at first with all it’s glitter and beauty too, but hopefully I would have uncovered that magic for you and revealed all the tricks by the time you finish reading this post. That is my goal. To make it as simple as possible for the keen to understand how the GPT-2 model works underneath.

Note: Pretty much the entirety of the code has been copied, inspired and referenced from Hugging Face’s implementation of the GPT-2, keeping merely the essentials for simplicity. If you want to train the GPT-2 model on parallel GPUs, save checkpoints while fine-tuning, run inference tasks on multiple CPUs and much more, I would recommend using the Hugging Face API. A simple tutorial on how to do so was recently released by Hugging Face and can be found here.

In this post, I am not trying to reinvent the wheel, but merely bringing together a list of prexisting excellent resources to make it easier for the reader to grasp GPT-2. I leave it up to the reader to further build upon these foundations in any area they choose.

You can’t build a great building on a weak foundation. You must have a solid foundation if you’re going to have a strong superstructure.

– Gordon B. Hinckley

## Prerequisites

This post assumes that the reader has a solid understanding of Attention and Transformers. The GPT-2 utilizes a 12-layer Decoder Only Transformer architecture. If you want a refresher or understand Attention and Transformers, here is an excellent list of resources to aid your understanding regarding:

1. The illustrated Transformer by Jay Alammar
2. The Annotated Transformer by Harvard NLP
3. Introduction to the Transformer by Rachel Thomas and Jeremy Howard

If you’re just beginning your journey into NLP or you’re an expert, I would definitely recommend the fast.ai NLP course taught by Rachel Thomas and Jeremy Howard. The course starts with the basics including Sentiment Classification using Naive Bayes and Logistic Regression, moves on to RNNs and also talks about Transfer Learning, ULMFiT, Seq2Seq translation and Transformers amongst other things. It is an excellent resource put together by the fast.ai team free of cost.

Another amazing resource on GPT-2 itself, is The Illustrated GPT-2 by Jay Alammar. This post starts with a basic introduction to Language Models and explains the GPT-2 model step-by-step in a very easy to understand manner. I would highly recommend the reader to give this post a read.

The Annotated Transformer by Harvard NLP implements the complete Transformer architecture using PyTorch and is great way to understand Attention in depth.

Let’s then build upon these excellent existing resources and implement GPT-2 in code.

## Abstract

Natural language processing tasks, such as question answering, machine translation, reading comprehension, and summarization, are typically approached with supervised learning on taskspecific datasets. We demonstrate that language models begin to learn these tasks without any explicit supervision when trained on a new dataset of millions of webpages called WebText. Our largest model, GPT-2, is a 1.5B parameter Transformer that achieves state of the art results on 7 out of 8 tested language modeling datasets in a zero-shot setting but still underfits WebText. Samples from the model reflect these improvements and contain coherent paragraphs of text. These findings suggest a promising path towards building language processing systems which learn to perform tasks from their naturally occurring demonstrations.

A Zero-shot setting is one where you do not finetune the language model and directly run inference on the target dataset. For example, pretrain a LM on WebText and directly try and predict the next words of Amazon Movie reviews dataset.

## Model Architecture (GPT-2)

We use a Transformer (Vaswani et al., 2017) based architecture for our LMs. The model largely follows the details of the OpenAI GPT model (Radford et al., 2018) with a few modifications. Layer normalization (Ba et al., 2016) was moved to the input of each sub-block, similar to a pre-activation residual network (He et al., 2016) and an additional layer normalization was added after the final self-attention block. We scale the weights of residual layers at initialization by a factor of 1/√N where N is the number of residual layers. The vocabulary is expanded to 50,257 words. We also increase the context size from 512 to 1024 tokens and a larger batchsize of 512 is used.

This is the entirety of model explanation inside the GPT-2 research paper. This warrants a need for us to look at the architecture inside the GPT model.

## Model Specifications (GPT)

Our model largely follows the original transformer work. We trained a 12-layer decoder-only transformer with masked self-attention heads (768 dimensional states and 12 attention heads). For the position-wise feed-forward networks, we used 3072 dimensional inner states. We used the Adam optimization scheme with a max learning rate of 2.5e-4. The learning rate was increased linearly from zero over the first 2000 updates and annealed to 0 using a cosine schedule. We train for 100 epochs on minibatches of 64 randomly sampled, contiguous sequences of 512 tokens. Since layernorm is used extensively throughout the model, a simple weight initialization of N(0, 0.02) was sufficient. We used a bytepair encoding (BPE) vocabulary with 40,000 merges and residual, embedding, and attention dropouts with a rate of 0.1 for regularization. We also employed a modified version of L2 regularization proposed in, with w = 0.01 on all non bias or gain weights. For the activation function, we used the Gaussian Error Linear Unit (GELU).

As can be seen from the GPT Architecture, to implement it, we will first need to implement Masked Self Attention and Feed Forward layer.

## Imports

import torch
import copy
import torch.nn as nn
import torch.nn.functional as F
from torch.nn.modules import ModuleList
from torch.nn.modules.normalization import LayerNorm
import numpy as np
import os
from tqdm import tqdm_notebook, trange
import logging
logging.basicConfig(level = logging.INFO)
logger = logging.getLogger()

## Transformer Decoder inside GPT-2

To re-use the terminology used to describe the Transformer, the attention is a function of a query (Q) and set of key (K) and value (V) pairs. To handle longer sequences, we modify the multi-head self-attention of the Transformer to reduce memory usage by limiting the dot products between Q and K in:

class Conv1D(nn.Module):
def __init__(self, nx, nf):
super().__init__()
self.nf = nf
w = torch.empty(nx, nf)
nn.init.normal_(w, std=0.02)
self.weight = nn.Parameter(w)
self.bias = nn.Parameter(torch.zeros(nf))

def forward(self, x):
size_out = x.size()[:-1] + (self.nf,)
x = torch.addmm(self.bias, x.view(-1, x.size(-1)), self.weight)
x = x.view(*size_out)
return x

### CONV1D Layer Explained

The CONV1D layer can be thought of as a LINEAR layer itself. Essentially, it is casting an initial tensor x (having the final dimension of x.size(-1)) being passed to it to have a final dimension of size self.nf.

Here’s an example output of the same:

d_model = 768
conv1d  = Conv1D(d_model, d_model*3)
x       = torch.rand(1,4,d_model) #represents a sequence of batch_size=1, seq_len=4 and embedding_sz=768, something like "Hello how are you"
x       = conv1d(x)
x.shape

>> torch.Size([1, 4, 2304])

As can be seen in the example above, the final dimension of tensor returned by CONV1D is 3 times the initial size. We do this to be able to cast the input to query, key and value matrices.

It is possible then to retrieve the query, key and value matrices like so:

query, key, value = x.split(d_model, dim=-1)

query.shape, key.shape, value.shape
>> (torch.Size([1, 4, 768]), torch.Size([1, 4, 768]), torch.Size([1, 4, 768]))

Another way to cast the input to Q, K and V matrices would have to been to have separate Wq, Wk and Wv matrices. I have explained this under the EXTRA section of this post at the bottom. I find this other approach more intuitive and relatable, but we use the CONV1D layer in this post, because we reuse the CONV1D pretrained weights from Hugging Face.

### FEEDFORWARD Layer Explained

class FeedForward(nn.Module):
def __init__(self, dropout, d_model=768, nx=768*4):
super().__init__()
self.c_fc    = Conv1D(d_model, nx)
self.c_proj  = Conv1D(nx, d_model)
self.act     = F.gelu
self.dropout = nn.Dropout(dropout)

def forward(self, x):
return self.dropout(self.c_proj(self.act(self.c_fc(x))))

Something, that’s just so well explained in Jay Alammar’s post - also referenced above, is how the inputs are passed through ATTENTION layer first and then on to FEEDFORWARD layer. The Feedforward network, is a normal neural network that accepts the outputs from the ATTENTION layer (768), casts them to nx (768*4) dimension, adds an activation function self.act (GELU), casts them back to d_model (768) and adds dropout (0.1).

This is also mentioned in the GPT research paper referenced below.

For the position-wise feed-forward networks, we used 3072 dimensional inner states

### ATTENTION Layer Explained

The below extract is from the paper Attention is all you need.

#### Scaled Dot-Product Attention

We call our particular attention “Scaled Dot-Product Attention”. The input consists of queries and keys of dimension dk, and values of dimension dv. We compute the dot products of the query with all keys, divide each by √dk, and apply a softmax function to obtain the weights on the values.

In practice, we compute the attention function on a set of queries simultaneously, packed together into a matrix Q. The keys and values are also packed together into matrices K and V . We compute the matrix of outputs as:

The two most commonly used attention functions are additive attention, and dot-product (multiplicative) attention. Dot-product attention is identical to our algorithm, except for the scaling factor of 1/√dk. Additive attention computes the compatibility function using a feed-forward network with a single hidden layer. While the two are similar in theoretical complexity, dot-product attention is much faster and more space-efficient in practice, since it can be implemented using highly optimized matrix multiplication code. While for small values of dk the two mechanisms perform similarly, additive attention outperforms dot product attention without scaling for larger values of dk. We suspect that for large values of dk, the dot products grow large in magnitude, pushing the softmax function into regions where it has extremely small gradients. To counteract this effect, we scale the dot products by 1/√dk.

To implement the The Attention layer in code, we first utilize the CONV1D layer and get the q, k and v matrices as explained before.

Once we have the q, k and v matrices, we can perform attention using the function _attn. This function replicates the formula mentioned above inside Attention Dot Product.

class Attention(nn.Module):
super().__init__()
self.d_model = d_model
self.c_attn  = Conv1D(d_model, d_model*3)
self.scale   = scale
self.softmax = nn.Softmax(dim=-1)
self.register_buffer("bias", torch.tril(torch.ones(n_ctx, n_ctx)).view(1, 1, n_ctx, n_ctx))
self.dropout = nn.Dropout(0.1)
self.c_proj  = Conv1D(d_model, d_model)

"return shape [batch, head, sequence, features]"
x = x.view(*new_shape)
return x.permute(0, 2, 1, 3)

def _attn(self, q, k, v, attn_mask=None):
scores  = torch.matmul(q, k.transpose(-2, -1))
if self.scale: scores = scores/math.sqrt(v.size(-1))
nd, ns  = scores.size(-2), scores.size(-1)
scores  = self.softmax(scores)
scores  = self.dropout(scores)
outputs = torch.matmul(scores, v)
return outputs

x         = x.permute(0, 2, 1, 3).contiguous()
new_shape = x.size()[:-2] + (x.size(-2)*x.size(-1),)
return x.view(*new_shape)

def forward(self, x):
x        = self.c_attn(x) #new x shape - [1,3,2304]
q, k, v  = x.split(self.d_model, dim=2)
out      = self._attn(q, k, v)
out      = self.c_proj(out)
return out

Another way to implement Attention is explained in the Extras section at the bottom of this blog. I find it to be more intuitive and easy to compare with the research paper. It utilizes Linear layers instead of CONV1D to cast inputs to Q, K and V matrices. The reason why we haven’t used it is because we use the pretrained weights for CONV1D layer from Hugging Face.

The below extract is from the paper Attention is all you need.

Instead of performing a single attention function with dmodel-dimensional keys, values and queries, we found it beneficial to linearly project the queries, keys and values h times with different, learned linear projections to dk, dk and dv dimensions, respectively. On each of these projected versions of queries, keys and values we then perform the attention function in parallel, yielding dv-dimensional output values. These are concatenated and once again projected, resulting in the final values, as depicted in Figure below.

Multi-head attention allows the model to jointly attend to information from different representation subspaces at different positions. With a single attention head, averaging inhibits this.

In this work we employ h = 8 parallel attention layers, or heads. For each of these we use dk = dv = dmodel/h = 64. Due to the reduced dimension of each head, the total computational cost is similar to that of single-head attention with full dimensionality.

Not to be confused by this, in essence all that’s being done is to add another dimension to the Q, K and V matrices. That is, if the matrices were before of size [1, 4, 768] which represents [bs, seq_len, d_model], these matrices are projected to dimension [1, 12, 4, 64] which represents [bs, n_head, seq_len, d_model//n_head]. GPT-2 utizes 12 parallel heads. We split the Q, K, V matrices inside split_heads function. Finally, once we get an output from applying parallel attentions we concatenate it inside merge_heads back to matrices of dimension [bs, seq_len, d_model].

## GPT-2 Model Architecture in Code

So far, we have implemented Multi Head Attention and FeedForward layers. The two layers form the building blocks of the Transformer Decoder block, shown in the picture above. The GPT-2 consists of 12 of these Transformer Blocks.

This has been shown in Jay Alammar’s post like so:

### Transformer Decoder Block Explained

class TransformerBlock(nn.Module):
super(TransformerBlock, self).__init__()
self.feedforward = FeedForward(dropout=0.1, d_model=768, nx=768*4)
self.ln_1        = LayerNorm(d_model)
self.ln_2        = LayerNorm(d_model)

def forward(self, x):
x = x + self.attn(self.ln_1(x))
x = x + self.feedforward(self.ln_2(x))
return x

The Transformer Block consists of Attention and FeedForward Layers. As referenced from the GPT-2 Architecture Model Specification,

Layer normalization (Ba et al., 2016) was moved to the input of each sub-block Here are the sub-blocks are Attention and FeedForward.

Thus, inside a Transformer Decoder Block, essentially we first pass the inputs to a LayerNorm followed by the first sub-block Attention. Next, we pass the outputs of this sub-block to LayerNorm again and finally to FeedForward layer.

## The GPT-2 Architecture Explained

As referenced from the GPT paper,

We trained a 12-layer decoder-only transformer with masked self-attention heads (768 dimensional states and 12 attention heads).

Thus, the complete GPT-2 architecture is the TransformerBlock copied over 12 times.

def _get_clones(module, n):
return ModuleList([copy.deepcopy(module) for i in range(n)])

class GPT2(nn.Module):
def __init__(self, nlayers=12, n_ctx=1024, d_model=768, vcb_sz=50257):
super(GPT2, self).__init__()
self.nlayers = nlayers
self.h       = _get_clones(block, 12)
self.wte     = nn.Embedding(vcb_sz, d_model)
self.wpe     = nn.Embedding(n_ctx, d_model)
self.drop    = nn.Dropout(0.1)
self.ln_f    = LayerNorm(d_model)
self.out     = nn.Linear(d_model, vcb_sz, bias=False)
self.loss_fn = nn.CrossEntropyLoss()
self.init_weights()

def init_weights(self):
self.out.weight = self.wte.weight
self.apply(self._init_weights)

def _init_weights(self, module):
if isinstance(module, (nn.Linear, nn.Embedding, Conv1D)):
module.weight.data.normal_(mean=0.0, std=0.02)
if isinstance(module, (nn.Linear, Conv1D)) and module.bias is not None:
module.bias.data.zero_()
elif isinstance(module, nn.LayerNorm):
module.bias.data.zero_()
module.weight.data.fill_(1.0)

def forward(self, src, labels=None, pos_ids=None):
if pos_ids is None: pos_ids = torch.arange(0, src.size(-1)).unsqueeze(0)
inp = self.drop((self.wte(src)+self.wpe(pos_ids)))
for i in range(self.nlayers): inp = self.h[i](inp)
inp     = self.ln_f(inp)
logits  = self.out(inp)
outputs = (logits,) + (inp,)

if labels is not None:
shift_logits = logits[..., :-1, :].contiguous()
shift_labels = labels[..., 1:].contiguous()
loss = self.loss_fn(shift_logits.view(-1, shift_logits.size(-1)), shift_labels.view(-1))
outputs = (loss,) + outputs
return outputs
return logits

Something I have not mentioned yet is Positional Encoding and Token Embeddings. Since, we cannot pass words such as “hey” or “hello” directly to the model, we first Tokenize our inputs. Next, we use Embeddings to represent the tokens as numbers. This post by Jay Alammar again explains Embeddings very well.

Also, since unlike the RNNs where the input words are passed sequentially, Transformers take input matrices in parallel thus losing the sense of position for the words being input. To make up for the loss, before handling the Token Embeddings to the model, we add Positional Encoding - a signal that indicates the order of the words in the sequence. Since, as mentioned before, the context size of GPT-2 is 1024, the positional encodings are of dimensions [1024, 768].

Thus, the inputs to the GPT-2 architecture is the sum of Token Embeddings and Positional Encodings passed through a Dropout, to add regularization. Once, we have the input matrix, we pass this through each of the 12 Layers of the GPT-2 architecure, where each layer is a Transformer Decoder Block that consists of two sublayers - Attention and FeedForward Network.

#### Language Modeling or Classification

When using GPT-2 as a language model, we pass the inputs to a final LayerNorm and through a Linear layer with a final dimension of size [768, vocab_sz] (50257) and get an output of size [1, 4, 50257]. This output represents the next word logits and we can very easily now pass this through a Softmax layer and take argmax to get the positional of the word inside the vocabulary with the highest probability.

For classification task, we can pass the outputs received from the GPT-2 architecture through a Linear layer with a dimension of size [768, n] to get probabilities for each category (where n represents number of categories), pass it through a softmax, get the highest predicted category and use CrossEntropyLoss to train the architecture to do classification.

And that’s really all the magic behind GPT-2. It’s a Decoder only Transformer Based architecture that takes inputs parallely with Positional Encodings unlike RNNs, passes them through each of it’s 12 Transformer Decoder layers (which consist of Multi head Attention and FeedForward Network) to return the final output.

Let’s see this model in action in a language model task.

## Sample text generation using Hugging Face Pretrained Weights

First, let’s initialize the model with the Pretrained Weights already provided by Hugging Face.

model = GPT2()
# load pretrained_weights from hugging face
# download file https://s3.amazonaws.com/models.huggingface.co/bert/gpt2-pytorch_model.bin to .

model_dict = model.state_dict() #currently with random initialization

old_keys = []
new_keys = []
for key in state_dict.keys():
if "mlp" in key: #The hugging face state dict references the feedforward network as mlp, need to replace to feedforward be able to reuse these weights
new_key = key.replace("mlp", "feedforward")
new_keys.append(new_key)
old_keys.append(key)

for old_key, new_key in zip(old_keys, new_keys):
state_dict[new_key]=state_dict.pop(old_key)

pretrained_dict = {k: v for k, v in state_dict.items() if k in model_dict}

model_dict.update(pretrained_dict)
model.eval() #model in inference mode as it's now initialized with pretrained weights

Let’s now generate text. We will utilize Hugging Face’s pretrained Tokenizer to convert words to input embeddings.

from transformers import GPT2Tokenizer
tokenizer = GPT2Tokenizer.from_pretrained("gpt2")
context   = torch.tensor([tokenizer.encode("The planet earth")])

def generate(context, ntok=20):
for _ in range(ntok):
out = model(context)
logits = out[:, -1, :]
indices_to_remove = logits < torch.topk(logits, 10)[0][..., -1, None]
logits[indices_to_remove] = np.NINF
next_tok = torch.multinomial(F.softmax(logits, dim=-1), num_samples=1).squeeze(1)
context = torch.cat([context, next_tok.unsqueeze(-1)], dim=-1)
return context

out = generate(context, ntok=20)
tokenizer.decode(out[0])

>> 'The planet earth is the source of all of all the light," says the study that the government will'

## Extras

Another way to implement Attention as shown in the NLP Course by fast.ai referenced from here, that I find to be more intuitive is as below:

class Attention_FASTAI(nn.Module):
super().__init__()
self.softmax  = nn.Softmax(dim=-1)
self.scale    = scale
self.atn_drop = nn.Dropout(0.1)
bias=bias) for o in range(3)]

x = layer(x)

def _attn(self, q, k, v, attn_mask=None):
scores  = torch.matmul(q, k.transpose(-2, -1))
if self.scale: scores = scores/math.sqrt(v.size(-1))
attn_prob  = self.atn_drop(self.softmax(scores))
attn_vec   = attn_prob @ v
return attn_vec

x         = x.permute(0, 2, 1, 3).contiguous()
return x.view(bs, seq_len, -1)

def forward(self, q, k, v, mask=None):
bs, seq_len = q.size(0), q.size(1)
wq, wk, wv  = map(lambda o:self.split_heads(*o, bs),
zip((q,k,v), (self.wq, self.wk, self.wv)))
attn_vec    = self._attn(wq, wk, wv)

The key difference between the implementation above and the one we have used is that this implementation does not use CONV1D. Instead, we first pass the input x to self.wq, self.wk and self.wv Linear Layers to get wq, wk and wv matrices and then perform attention as before.


text:
## SWAG Group Title how o u solve (3x^2y)^3 one year ago one year ago

1. Mertsj Group Title

$(3x^2y)^3=3^3(x^2)^3(y^1)^3$

2. Mertsj Group Title

Does that help?

3. Mertsj Group Title

What is 3 to the third?

4. campbell_st Group Title

the power of 3 is operating on each part of the expression and can be written in expanded form as $3^3 \times (x^2)^3 \times (y^1)^3$ to simplify the variables use the index law for power of a power which means multiply the powers $(x^a)^b = x^{a \times b}$

5. SWAG Group Title

thank u guys

format: plain


text:
## De Rham Cohomology

De Rham Cohomology is a very cool sounding term in advanced math. This blog post is a short introduction on how it is defined.

Definition:
A differential form $\omega$ on a manifold $M$ is said to be closed if $d\omega=0$, and exact if $\omega=d\tau$ for some $\tau$ of degree one less.

Corollary:
Since $d^2=0$, every exact form is closed.

Definition:
Let $Z^k(M)$ be the vector space of all closed $k$-forms on $M$.

Let $B^k(M)$ be the vector space of all exact $k$-forms on $M$.

Since every exact form is closed, hence $B^k(M)\subseteq Z^k(M)$.

The de Rham cohomology of $M$ in degree $k$ is defined as the quotient vector space $\displaystyle H^k(M):=Z^k(M)/B^k(M).$

The quotient vector space construction induces an equivalence relation on $Z^k(M)$:

$w'\sim w$ in $Z^k(M)$ iff $w'-w\in B^k(M)$ iff $w'=w+d\tau$ for some exact form $d\tau$.

The equivalence class of a closed form $\omega$ is called its cohomology class and denoted by $[\omega]$.


text:
# Event¶

The Event type represents events that can be applied to contracts. The Event type is built-in. The definition is given below.

type Event {
agent : Agent,
timestamp : DateTime
}

## ContractInstance¶

The built-in abstract type ContractInstance represents a reference to an instantiated contract. Values of this type cannot be constructed directly in CSL syntax, so whenever such a value is needed it must be supplied to the runtime system directly – either as an instantiation argument to a template or as an argument to a report.


text:
# Preprocessing survey data for clustering

I want to find 4-10 clusters in survey data with 100 questions answered by 2000 individuals using a technique such as K-means or Gaussian Mixture Models. There is no response variable so the clustering technique needs to be unsupervised. The outcome of the clustering highly depends on the way the data is represented and preprocessed, causing me to identify multiple potential issues.

The responses to the questions can be either ordinal with M possible values or categorical with N possible values. M will be either 3 or 10, whereas N can be anything in the interval [2,8]. All questions have been answered by all users.

Below are the issues that I would like to get input on:

1. How should find the similarity between responses when I have these different kinds of variables? I would assume that the best way to input the ordinal information would be as a single variable/dimension, whereas a categorical variable will be represented using N variables/dimensions. Should I scale the variables in a particular way to ensure equal weight to each question or is there some similarity measure I could use that fits this kind of data?

2. The questions may be subgrouped categories such as "Sociodemographic Information", "Communication Assessment", "Product Assessment" and "Brand Assessment". Suppose that "Communication Assessment" contains 80% of the questions, but we would still like to get a clustering result that uses equal amount of information from each category. How should I preprocess the data to get this result?

3. The surveys are sent out to a representative sample of the population. However, the responses are biased such that males are overrepresented and females are underrepresented. The difference is statistically significant. How should I include this information in the clustering process? Should I oversample users to get a representative dataset across gender, age, job type? The number of duplicate users created this way would probably be 10-15% of the total.

You should feel free to address other issues with clustering survey data that I have not identified.

A few suggestions:

Wrt (#1), you have a mixture of scales and information. Therefore, K-means (which assumes continuously distributed information only) is not an appropriate technique. Depending on the mixture of scale types you may want to consider traditional PCA (for continuous data only) vs an approach that integrates an unsupervised mixture of scales such as correspondence analysis or latent class PCA.

Latent class models are another approach to clustering but are supervised and presume the existence of a response or target variable, which you may not have available.

In terms of your concern with (#2) "equal weighting" and pre-processing of the various features in your data (e.g., communication), you can create a separate dimension for each one either using judgement or based on a principal components analysis to collapse and combine the features appropriate for each.

To obtain results that are weighted and representative, I wouldn't oversample. Consider developing a weight factor instead. The general rule of thumb is not to weight the algorithms while factoring or clustering your data, but to weight the final results of those algorithms. There are lots of ways to weight your data, one approach is to use an IPF (iterative proportional fitting) model based on the marginals of each factor used.

• Thanks! Sadly, I don't have any response or target variable. I'm looking into correspondence analysis - thanks for the suggestion! However, I cannot find anything on latent class PCA, but only on latent class models? Could you explain the arguments behind that general rule of thumb? – pir Nov 25 '15 at 21:37
• Latent class PCA is covered in these publications (under Discrete Factor Models) from the vendors of Latent Gold, an LC software tool. I'm not pushing LG but having used it for years, it's the best developed and supported statistical routine for this class of models that I'm aware of. And, at ~\$1k per license, it's not expensive. statisticalinnovations.com/publicationsarticles/#dfm – Mike Hunter Nov 25 '15 at 21:40
• Sadly, I'm still a student and do not have that kind of money. Would it make sense to compute similarities separately for continuous and categorical variables and then combine then when clustering? – pir Nov 26 '15 at 11:58
• I wouldn't do them separately. You can do correspondence analysis in R, I'm sure, although I can't recommend a specific module. Or if you have access to a package such as SPSS or SAS through your school, they both have modules for this. – Mike Hunter Nov 26 '15 at 12:16
• Based on en.wikipedia.org/wiki/Correspondence_analysis it seems that correspondence analysis needs to be applied to only categorical data? I also have continuous variables. – pir Nov 26 '15 at 15:42


text:
# Tag Info

### Why do compounds like SF6 and SF4 exist but SH6 and SH4 don't?

TL;DR Fluorine is electronegative and can support the extra negative charge that is dispersed on the six X atoms in $\ce{SX6}$, whereas hydrogen cannot. First, let's debunk a commonly taught myth, ...
• 65.5k
Accepted

### Why is the bond angle H-P-H smaller than H-N-H?

Starting point: 2s orbitals are lower in energy than 2p orbitals. The $\ce{H-N-H}$ bond angle in ammonia is around 107 degrees. Therefore, the nitrogen atom in ammonia is roughly $\ce{sp^3}$ ...
• 81.6k
Accepted

### Why is [PdCl4]2- square planar whereas [NiCl4]2- is tetrahedral?

The geometry of the complex changes going from $\ce{[NiCl4]^2-}$ to $\ce{[PdCl4]^2-}$. Clearly this cannot be due to any change in the ligand since it is the same in both cases. It is the other factor,...
• 65.5k
Accepted

### What are t2g and eg in CFT?

As stated in both the links Geoff and Philipp have kindly commented (1, 2) they are to do with symmetry labels we chemists like to assign to orbitals. Knowing an orbitals symmetry class can lead to a ...
• 1,973
Accepted

### What is a non-classical carbocation?

Here is a picture of a "classical" carbocation, there is an electron deficient carbon bearing a positive charge. There are many examples of "non-classical" carbocations, but the 2-norbornyl ...
• 81.6k

### When are two orbitals orthogonal?

Unfortunately, the sense in which orbitals are orthogonal is more or less impossible to define rigorously without recourse to functions of some kind. So, I'll give an explanation a shot using some ...
• 16.8k

### When is it true that more nodes equals higher energy?

General case There is indeed a mathematical theorem that deals with the number of nodes an eigenfunction corresponding to a certain eigenvalue can possess. It was laid down by Courant$^{[1, 2]}$ and ...
• 17.2k

### Why do n AOs only form n MOs?

If you have $n$ functions (e.g. AOs) you can make a maximum of $n$ new linearly independent functions (e.g. MOs). If you try to make $n+1$ MOs, then any one of them can be expressed as a linear ...
• 4,759

### Molecular orbital (MO) diagram for N2 and N2^-

I have been taught that the MO diagram is different for molecules with 14 or less electrons than the one used for molecules with 15 or more electrons. This is (partly) wrong because the change in ...
• 18.4k
Accepted

### Canonical MOs vs. Localized MOs: Do both represent reality in the same way?

NOTE: In the below, I'm implicitly discussing a ground-state, closed-shell wavefunction, where all occupied orbitals are doubly occupied. The discussion would be similar for open-shell wavefunctions, ...
• 16.8k

### How can the dipole moment of carbon monoxide be rationalised by molecular orbital theory?

Unfortunately, nothing in the bonding situation in carbon monoxide is easily explained, especially not the dipole moment. According to the electronegativities of the elements, you would expect the ...
• 42.2k

### Physical intuition behind negative values for wave function?

The wavefunction of a particle actually has no physical interpretation to it until an operator is applied to it such as the Hamiltonian operator, or if you square it which gives its probability of ...
• 7,036
Accepted

### Cyclobutadiene - Jahn–Teller effect or not?

Very interesting question, and it kept me up despite daylight saving time cheating me of one hour of sleep last night... A good reference is Albright, Burdett and Whangbo, Orbital Interactions in ...
• 65.5k

### Why is the bond length of CO+ less than that of CO?

Unfortunately, the arguments presented by buckminst and Uncle Al aren't completely right. The MO schemes are correct but the ...
• 17.2k
Accepted

### How do atomic orbitals interfere both constructively AND destructively to give two molecular orbitals?

What does MO formation entail? A very common misconception is that the formation of MOs involve addition, or subtraction, of two physical objects. And the confusion arises because, how can the same ...
• 65.5k
Accepted

### What is actually the difference between valence bond theory and molecular orbital theory?

TL;DR VB theory treats atomic orbitals (including hybridized orbitals) as providing a good mathematical/physical description of the true form of the molecular wavefunction. MO theory uses atomic ...
• 16.8k

### Sp5 hybridization in cyclopropane?

Until this post, $\mathrm{sp^1,~sp^2~\text{and}~sp^3}$ meant for me, that the hybrid orbitals would consist from one s orbital and one, two or three p orbitals. I thought of it being so, as this seems ...
Accepted

### How can the gauche-effect be explained?

The gauche effect is commonly explained with LCAO-based bond orbitals. LCAO is short for linear combination of atomic orbitals and implies that we can take two atomic orbitals $\phi_1, \phi_2$ and ...
• 63.6k

### Destabilization of antibonding MO vs stabilization of bonding MO

TL;DR "The destabilisation of an antibonding MO is equal to the stabilisation of a bonding MO, relative to the constituent AOs." Several key assumptions are inherent in this statement. Born–...
• 65.5k
Accepted

### Why is the bond order of CO+ 3.5?

For a long time it was taught in school and universities that the HOMO of carbon monoxide is anti-bonding. Without more context it was also often taught that the bond order in CO is three, since there ...
• 42.2k
Accepted

### When are two orbitals orthogonal?

One atomic (or molecular) orbital is said to be orthogonal to another atomic (or molecular) orbital if there is no interaction between the electrons in one orbital with the electrons (wavefunction) in ...
• 81.6k

### Delocalization of π-electrons in nitrate ion

Unfortunately the key to understanding the delocalisation of electrons lies in understanding rudimentary molecular orbital theory. It is also necessary to understand that hybridisation is a ...
• 42.2k
Accepted

### What is the usage of orbitals more complex than f orbitals?

Surprisingly, I learned that there are also usages for orbitals g,h,i and even j. Actually, the letter "j" is not used, so it is s, p, d, f, g, h, i, k, l, etc. The higher angular momentum ...
• 38.8k
Accepted

### What does it mean when it is said that an sp³ orbital has 25% s character?

Hybridisation is a purely mathematical concept, which makes it possible to explain experimentally found structures. The most prominent example for this is methane, where you can consider the central ...
• 42.2k
Accepted

### Jahn-Teller Distortions in Square Planar Complexes?

Introduction Thank you for prompting me to look for square planar complexes that are not $\mathrm{d^8}$; I learnt some valuable stuff while researching the answer to this question! Do take this ...
• 63.6k
Accepted

### Can the idea of entropy be extended to orbitals?

No. The reason why a gas particle in a large volume has a large entropy is not because it has a lot of space to move around per se. A better explanation is that for a given energy, there are many ...
• 65.5k
Accepted

### How can the intense color of potassium permanganate be explained with molecular orbital theory?

Let’s take a look at a qualitative MO scheme for a tetrahedric transition metal complex whose ligands have three p-type orbitals each. On the left of figure 1 you have the metal orbitals ($\mathrm{3d}$...
• 63.6k
Accepted

### Why/how is blood red? (colours of hemoglobin)

The $\ce{Fe}$-porphyrin complex is a vastly conjugated $\pi$-system, in which the oxidation state of the $\ce{Fe}$ plays a central role. The wavelength of light absorbed by hemoglobin exactly ...
• 23.2k

### Does a neutral dimercury molecule exist?

There has been recent research on mercury dimers in its ground state. Bond energies, transition energies, band spectrum and other spectroscopic parameters has been calculated. Here are the abstracts ...
• 21.7k
No 6-centre-10-electron bond. The bonding situation in $\ce{IF7}$ is explained by Christe et. al. (see ssavec's answer) as a 6-centre-10-electron bond. In the article itself I found no real ...


text:
# Normal distribution and star density in a Galaxy

1. May 19, 2013

### corndog16

Hey all.

I'm working on a personal programming project where I'm attempting to simulate (to a small degree) a galaxy. And I have come across a decent 2D density map for a spiral galaxy. This map (array actually) defines a 128x128 grid of values between 0 and 255 representing the frequency of stars found in that 'sector' of space.

I want to expand this density map into the 3rd dimension and I figure a good way to do that would be to define the height of the galaxy as being n sectors tall and then use a normal distribution curve to distribute the density value (something between 0 and 255) vertically across those 5 sectors.

My question is, how would I take a number, say 234, and divide it among n 'boxes' or 'sectors' so that it has a normal curve type distribution?
I attempted to Google this, but I found that what was coming up did not describe what I was looking for so I'm thinking that I don't currently have the vocabulary to describe this particular intent.

2. May 20, 2013

### Stephen Tashi

A normal distribution has two parameters, the mean $\mu$ and the standard deviation $\sigma$. Lets say we have 128 boxes. Imagine the boxes as spaces between a ruler that has marks at 0,1,2,..128. Assuming you want the middle of the distribution to be at the middle of the ruler, set $\mu = 64$.

Then you have to decide on a value for $\sigma$. A normal distribution has non-zero values from $- \infty$ to $+\infty$, so you have to decide what part of the distribution you are going to leave out or you have to settle for a distribution that doesn't exactly match a normal distribution by forcing all of the stars to be within the 128 boxes.

For most programming languages, you can find a library with a function that give the values of the cumulative normal distibution. Let's assume that such a function is

$F(\mu,\sigma,x) = \int_{-\infty}^x f(\mu,\sigma,x) dx$ where $f(\mu,\sigma,x)$ is the density of a normal distribution.

The fraction of stars between marks $A$ and $B$ on the ruler is

$F(64,\sigma,B) - F(64,\sigma,A)$.

So the fraction of stars between the 0 and 128 marks is

$M = F(64,\sigma,128) - F(64,\sigma,0)$

The fraction that are outside those marks is $1 - M$.

The fraction of stars in the kth box is $s_k = F(64,\sigma,k) - F(64,\sigma,k-1)$. If you want to "renormalize" so that all the stars fall in the marks between 0 and 128, you can make the fraction of stars in the kth box $\frac{s_k}{1-M}$

For 234 total stars, the number of stars in the kth box will be $(s_k)(234)$. So it won't always be a whole number. To convert to a while number, you'll have to decide how to round things off. (If you want the grand total to remain 234, you'll have to round off in a special way!)

The above calculations don't determine a unique value for $\sigma$. You'll have to try various values to find one that you like if your goal is only to have a pleasing picture. If you are doing a physical simulation of moving stars, we'll have to investigate the problem further.

3. Nov 27, 2013

### davidhayes

Hi,
Any chance you could post your source for the 2D Density Map you have? I'm looking for exactly that
Thanks
Dave


text:
Jukes-Cantor Probabilities

In an earlier post I went through a derivation of the Jukes Cantor model (based on a Poisson distribution of mutation events) and got the following type of relationship.

,

where  is the fraction of divergent sites between two DNA sequences and  is the time the sequences diverged from a common ancestor in units of  where  is the number of generations and  is the per generation per nucleotide mutation rate.

This can be rearranged to

to convert an observed difference between two sequences into a corrected estimate of the divergence between two species.

However, this is only a point estimate and cannot handle cases where the fraction of sites that differ happen to be above 75% (this is possible by chance but becomes a log of a negative number and undefined).

We can rewrite this in a different way, as a binomial probability, and evaluate the probability of a data set over a range of divergence and place intervals around our point estimate.

As above, the probability that two nucleotides are different is

.

Assuming each nucleotide is independent in a comparison of two sequences we can multiply by the number of sites that are different, .

On the other hand two nucleotides can be in the same state if no mutations occurred, which is a probability of  (this comes from the Poisson distribution the model is based on, look back at the earlier post if this isn't clear) or mutations could have occurred along the lineages but they happened to end up as the same nucleotide (an "invisible" mutation). This second outcome has a probability of

.

Combined the probability that we see  sites that are the same in two sequences descended from a common ancestor is

.

Since there are only two possibilities, the sites are the same or different, and the probability of all possible outcomes must be 100% we can rewrite this in a slightly simpler way.

.

Then multiply everything together to get the probability of the data including sites that differ and sites that are the same.

There is also a binomial coefficient that should be multiplied to get the exact probability (because there are different ways to get the same data: the first and second site can be different, or the first and third, or the second and fifth, etc.). However, for a given data set this is a constant and we're going to drop it for the moment.

Say we align two sequences of ten base pairs from the same gene in two different species.

grasshopper: AGCTACAACT
cricket:     AGATACGACT

The sequence differs at two sites. We could rewrite the comparison as a 1 if the base is the same and a 0 if they are different.

DNA: 1101110111

Making a plot of the two parts of the equation we can see that the probability that two sites are the same goes down as the time of divergence increases and the probability that they are different increases.

Importantly, there is a middle ground that can best accommodate both signals from the data, the differences and the similarities, near but lower than the crossover point. Combining these together we get the curve of probability of all the data with a peak just above 0.1.

This can be integrated to find upper and/or lower confidence intervals (which is messy and uses the hypergeometric function); however, it is clear from looking at the curve that a wide range of values are hard to rule out with a data set this small.

If the data set were larger we can narrow down the confidence interval.

So, how does this compare to the point estimate we get from the traditional Jukes-Cantor distance? A difference of 0.2 (two sites out of 10) can be plugged in for ,

,

and gives .

When we take the derivative of

set it equal to zero and solve for  we get

.*

Plugging in  and  gives , exactly the same as with the other method. The original Jukes-Cantor method is a maximum likelihood point estimate (the parameter value of the model that maximizes the likelihood of the data). You can also see it is equivalent in other ways. The  in the denominator makes the log evaluate a negative number if the fraction of sites that vary are greater than 3/4 of the total, which is undefined. (Thinking about this intuitively it is most consistent with positive infinity; however, I suspect the curve gets infinitely flat and almost any large number is essentially just as likely.) However, we can still work with the probabilities which allows us to rule out small values of  as in the example below where the sites that differ are 80% of the total.

* This is only true for real solutions. The more general solution is

,

where  is an element of the integers. This is fun because it combines a base  log with  and , but this isn't really relevant to the discussion (time of divergence doesn't rotate into a second dimension).

Toxorhynchites

We found a total of five Toxorhynchites larvae in a pool containing Culex larvae near Manoa stream (Jan 18, 2018, 21.299882, -157.813210). I haven't seen this genus before; they are possibly T. splendens. These are unusual mosquitoes in that they do not blood feed. The larvae eat aquatic invertebrates including other mosquito larvae and several species were introduced into Hawaiʻi as an attempt at biocontrol of mosquitoes. They are also very large mosquitoes and the genus includes the largest mosquito species.  Thanks to Matt Medeiros and Jessica Mavica for helping me ID them.

Jessica Mavica, Pirbright Institute

Jessica Mavica visited our lab last week and gave a presentation in our department seminar series. She works at the Pirbright Institute (https://www.pirbright.ac.uk/) and took some Culex mosquito larvae back to the institute for research projects there.

Transcriptional effects of a positive feedback circuit in Drosophila melanogaster

A new paper is out at BMC Genomics; it was published a few days before the new year.

Transcriptional effects of a positive feedback circuit in Drosophila melanogaster

This is work that was started many years ago. The tTAV gene expression system was developed into a genetic pest management tool by various groups (e.g., Horn C, Wimmer EA. A transgene-based, embryo-specific lethality system for insect pest management. Nat Biotechnol. 2002;21:64–70; Alphey, Luke. "Expression system for insect pest control." U.S. Patent 9,121,036, issued September 1, 2015). A positive feedback form of the tTAV system was being marketed by Oxitec as "RIDL."  In this system a transgene promotes its own expression by binding to DNA. If left unchecked this causes a runaway feedback loop and is toxic to the organism. However, the protein produced also binds to tetracycline and if tetracycline is added to the insect's diet then the lethal effect is masked (tetracycline bound protein does not bind to DNA and drive gene expression).

In a nutshell, in applications large amounts of insects can be reared with tetracycline. Then released into the wild to mate with wild individuals. The toxic effect is dominant and kills off the offspring of the wild individuals that inherit it (who presumably do not have tetracycline in their diet). If enough individuals can be released relative to the wild population this can cause a reduction in wild population numbers over the following generations. This is similar to classical SIT (Sterile Insect Technique) for population suppression.

We were curious about testing this system in various ways. In the first exploration we added tetracycline to food for fruit flies to see if there was an effect. Tetracycline can disrupt normal mitochondrial function and kill naturally occurring bacteria including endosymbionts like Wolbachia. This was the basis of a student's research project: Müller, H., Krefft, M., Reeves, G., & Reed, F. A. 2010. Trans-generational influence of tetracycline on Drosophila melanogaster (Bachelor Thesis, Fachhochschule Bingen). He found a reduction in egg laying and a delay in development in response to tetracycline. Interestingly there was a transgenerational effect in both the paternal and maternal sides (maternal could be due to mitochondrial or Wolbachia effects, but these do not explain the paternal effect).

Next we engineered the tTAV system in Drosophila melanogaster and performed an EMS (ethyl methanesulfonate) mutation screen to see if any mutations could result in tolerance of the feedback loop and survival off of dietary tetracycline. This was the basis of a grant application: Deutshe Forschungsgemeinschaft (DFG , German National Research Foundation). Die Entstehung von Resistenzen gegen genetisch induzierte Sterilität bei Insekten. (The evolution of resistance to genetically induced sterility in insects.) 2010 RE-3062/2-1. Long story short we generated ~400,000 embryos from flies exposed to EMS and recovered five individuals that had the tTAV insert but could survive without tetracycline.  All five were male. We tried mating them to females to recover offspring---and all five appeared to be sterile which essentially ended that experiment.

In this paper we wanted to determine the genome-wide transcriptional response to the runaway tTAV effect positive feedback loop.

The tTAV insert works as expected in D. melanogaster. The image above shows the curve of lethality as the concentration of two antibiotics, doxycycline and tertracycline, declines. At the upper end of the curve approximately half of the offspring of a heterozygous parent carry the genetic modification, as expected. As the dietary concentration diminishes the curve approaches zero---none of the offspring that carry the tTAV insert survive.

The cause of this lethality is unknown. One can imagine that a positive feedback loop of gene expression and protein production may over-use transcriptional and translational resources and/or overload protein degradation pathways, or disrupt off target gene expression. These hypotheses are not mutually exclusive. Looking at the shift in expression of other genes in the genome in response to tTAV runaway expression could be informative and give us clues as to what process results in lethality.

While we did see thousands of genes across the Drosophila genome that significantly changed expression, there was no clear pattern that emerged. The tTAV system was inserted at different points in the genome and these different lines gave different groups of genes with altered expression. A relative handful of 31 genes with differences in expression were shared across all lines. However, 27 of these 31 genes shifted in expression in different directions, some with increased expression and some with decreased expression across the lines.

This left just four genes with consistent changes in expression when the tTAV system was activated: crok (CG17218), Cyp6a17 (CG10241), olf186-F (CG11430), and Pex23 (CG32226). There is no clear relationship between these genes nor any obvious hypothesis that presents itself in terms of understanding a response to shifts in tTAV activity.

We backup up and looked at gene ontology enrichment. Were there categories of genes involved in certain processes or pathways that tended to be represented in the shift of gene expressions associated with tTAV---again no clear pattern presented itself.

Were there regional changes in expression, perhaps related to different insert sites or nuclear architecture interactions among chromosomes? No, we looked at groups of genes in sliding windows along the genome and there was no clear pattern of clustering by region.

Keep in mind however, there were widespread and significant changes in gene expression. It was not that there was the lack of an effect; just that there is not a clear pattern to the effect.

Gene expression levels between pairwise comparisons of the different insert lines and the control all showed statistically significant correlations. However, the direction of the correlation was not consistent. (This is shown in the figure above.)

I was quite surprised by this result---that the insertion site of a highly expressed transgene has such a profoundly idiosyncratic effect on genome-wide transcription patterns. This was a bit frustrating to be honest but it does suggest some cautions in these kinds of studies and next steps to explore. First of all a a caution, if we had only studied the transcriptional response to a single insert we would have found a significant response, with many genes affected, and may have been misled by the resulting pattern (if for example a number of these genes were in a certain pathway), overgeneralizing the particular signal seen. In a next step, it would be interesting to see if this widespread idiosyncratic effect extends to other genetic modifications or is specific to the tTAV system. For example, what is the genome-wide translational response, if any, to actin5c driven expression of EGFP? Importantly, does any effect seen vary widely by insertion site?

Gene Drive Questions

I received an inquiry about gene drive technology by email and thought that it might be useful to post part of it along with my reply here.

Here is the question:

"Since a few months, gene drive is a very prominent topic. But what I am missing a recent review about the state of the art in building gene drive constructs and testing them. I would like to ask you if you could provide some advice about literature, which would cover an overview about this theme.  At the level where I try to follow the discussions, I have the feeling that things are randomly mixed (like e.g. gene drive = synthetic biology). It would be really great to have an overview about what is possible by now, what has been done, which systems work and which do not work. And perhaps where is gene drive applied."

And my response:

Here are some general reviews:

This is rapidly developing.  A lot of the recent focus is on CRISPR-Cas9 based gene drive and there is some confusion in the media that this is the only form of gene drive.  I've tried to categorize the different forms in Appendix A of this article:

https://arxiv.org/abs/1706.01710

One major problem with CRISPR-Cas9 based gene drive is the rapid accumulation of mutation and naturally occurring resistance in the wild. See:

Work is also progressing on underdominance-type approaches that may get around some of the problems with the CRISPR system.

The Genetic Engineering and Society Center at NCSU has a lot of resources and knowledgeable people in this area.

https://research.ncsu.edu/ges/

Here is a series of recent publications that were recently posted---I haven't had a chance to go through all of them yet.

https://research.ncsu.edu/ges/publications/faculty-publications/?table_filter=JRISpecialIssue

There are also some relevant publications that will come out over the next few months. I can update you as these are published and become publicly available.

To my knowledge there is no direct application as of this moment, but this is likely only a matter of time.

The US military is also very interested in gene drive issues. DARPA is funding some of the research on gene drive and IARPA is working on ways to detect if genetic modifications have occurred.

Let me know if I can help with additional information or if you want more detail about something.

Virus imaging with negative staining

First of all there has been a long lag in posting any updates here about the lab.  There are many things going on that are a higher priority at the moment; however, I also want to try to keep a certain tempo of posts and updates going here if possible.  So, I set aside a few minutes today to make a post.

A couple of undergraduates in the lab have an interest in viruses, specifically bacteriophages. They have worked out methods to isolate and propagate the phages for various experiments. However, we wanted to get a look at the actual viral particle instead of inferring it indirectly from plaques on a plate of bacteria. I talked to Tina Carvalho at our Biological Electron Microscope Facility and she was enthusiastic about imaging the viruses.

Negative staining is a relatively quick and easy method to get viral mugshots. The virus particles are put onto a thin membrane held in a metal grid and electrons are beamed through both the virus and the grid to get an image. The problem is that this doesn't really work that well.  The virus is very transparent to electrons and they will pass right through almost like it wasn't there, i.e., the image contrast defining the virus will be low and look faint.

A stronger image can be made by adding a thin fluid of uranium salt (Uranyl acetate (UO2(CH3COO)2·2H2O)). Uranium is very dense and good at blocking electrons. If all goes well (and Tina pointed out that there can be quite a bit of luck involved with the color of you socks that day possibly playing a role) the fluid will dry a bit and form droplets around the virus particles by adhesion. The electrons pass through the virus, which is making a path of sorts, through the uranium drop, so the outside shell of the viral surface can be better resolved. Here is a figure to try to illustrate.

First we looked at a Vibriophage that Stacy Paulino isolated from seawater from Maunalua Bay. This phage infects a coral pathogen, Vibrio coralliilyticus, and Stacy imaged it below.

The next image is an unknown coliphage (infects Escherichia coli).

Notice the scale in the image indicating 20 nanometers; that is the size of a ribosome. Double stranded DNA is about 2 nm wide of 1/10 of the scale bar. The graininess of the tube is from individual tail tube proteins, single polypeptide molecules, which are typically assembled from ~250 individual peptides.

Maya Shaulsky (in a collaboration with Bob Thomson) is working with T7 (below) and made these coliphage images. Originally we thought we were working with T7 coliphage but we kept getting frustrating results that didn't make sense. When we actually saw the virus (the image above) it confirmed that we were not working with T7 so we ordered actual T7 from a different stock center.

And here is an E. coli bacterium with T7 particles replicating inside the cell!

If you are working with something small I encourage you not to be intimidated with the idea of electron microscopy imaging. It is really something to actually see what you are working with. The protocol to prepare the specimens is not that long or difficult (some people may disagree with this, perhaps we have been lucky) and it is not really that expensive to pay for the training and machine time. It is worth trying at least once.

An alternative GAL4-based genetic sterile insect technique?

Here is a thought that came to me on a walk. Perhaps this could also be applied to mosquitoes to suppress wild populations.

GAL4, UAS, and GAL80 are genetic tools from yeast that are commonly used in Drosophila.

GAL4 is a protein that drives expression of genes with UAS (Upstream Activation Sequence) as a gene expression enhancer. GAL80 inhibits GAL4.

So, place the gene encoding the GAL4 protein under UAS expression. This could create a runaway positive feedback loop (GAL4 drives its own expression of more GAL4) and result in lethality.  However, set this up in Drosophila that also express GAL80 to inhibit the GAL4 from binding to UAS.

If GAL80 were expressed from a Y-chromosome only male offspring would survive to mate with the remaining females. Female-specific lethality is a powerful way to suppress a population.

This could be set up in a stable lab Drosophila stock using compound X chromosomes for X^X Y females and regular X Y males. All the Drosophila stocks, sequences, and insertion sites exist to quickly try this out.

http://flystocks.bio.indiana.edu/Browse/aberration/compound_x_overview.htm

http://flystocks.bio.indiana.edu/Browse/chr_y/special_y.php

https://www.thebestgene.com/

Sea urchin mitochondrial genome evolution

A lab publication is out (link)!

Áki J. Láruson generated the mitochondrial genome of Tripneustes gratilla and compared it to all other existing sea urchin genomes. He found that the Toxopneustidae are a sister family to the Strongylocentrotidae splitting from each other around a time of major climate change (link) and that patterns of DNA substitutions are generally consistent with neutrality for the sea urchin mitochodria. Stay tuned for more Tripneustes gratilla genomics publications!

NIH capping fund distributions to address problems in science in the US

In an article recently posted at NIH a plan to limit individual grant support has been announced.

Implementing Limits on Grant Support to Strengthen the Biomedical Research Workforce

NIH is facing the fact that (1) there is not enough grant funding to go around, (2) the funding is being concentrated into fewer labs (3) at fewer institutions, and as a consequence (4) science is loosing young researchers because the field appears to have no future for people starting out. Without a next generation, science in the US is approach a crash in its workforce. In response NIH plans to spread out research funding among more researchers by placing a cap on individual funds received.

"And we will monitor, on a trans-agency basis, investigators’ Grant Support Index, with the idea that over time and in close consultation with the extramural research community, we will phase in a resetting of expectation for total support provided to any one investigator. We plan to implement a Grant Support Index cap of 21 points, essentially the equivalent of 3 single-PI R01 grants. Over the next few weeks to months, we will meet with NIH Advisory Councils and other stakeholder groups to explore how best to phase in and implement this cap – so that formal assessment of grant support can be used to best inform, on a trans-NIH basis, our funding decisions." link

Why stop at three? I would be overjoyed with just one or two R01 grants to allow me to do my job.

What is science?

We can read the steps of the scientific method and know about conventions such as the peer review process but at its core what is science? I have been thinking about this question more and more recently. In trying to pin it down the definition seems to diffuse into a overly abstract statement that resembles something I dislike as a shortcut to thinking and refer to as "bumper sticker philosophy." It seems odd to phrase it this way but science is, in a way, about what works and what doesn't work, and understanding the reasons why what works and what doesn't work in a very rigorous way. Science, ideally, is a careful balance of keeping an open mind, embracing curiosity, taking details seriously, and reasoned agreement and disagreement with others. Then there is the real world human aspect of how science is or is not practiced. Statements such as "believing" in evolution or "believing" in climate change do not sit well in a scientific perspective. In current American English this implies faith in something... However, pointing this out can be misinterpreted and hijacked by people with non-scientific motivations.

At the risk of sounding like a kōan, I thought two famous historical case examples, one from physics and one from biology, that have interested me for different reasons, might provide more of an illustration of what science is then trying to nail down a precise verbal/text definition. However, this comes at the risk of people not being patient enough to read them. Perhaps therein lies part of the modern problem. If it doesn't fit into a snappy short phrase it never reaches most people. I've linked PDFs of the two examples below.


text:
# BMO and Elasticity: Korn's Inequality; Local Uniqueness in Tension

created by spector on 08 Apr 2020

[BibTeX]

Preprint

Inserted: 8 apr 2020

Year: 2020

ArXiv: 2004.02368 PDF

Abstract:

In this manuscript two $BMO$ estimates are obtained, one for Linear Elasticity and one for Nonlinear Elasticity. It is first shown that the $BMO$-seminorm of the gradient of a vector-valued mapping is bounded above by a constant times the $BMO$-seminorm of the symmetric part of its gradient, that is, a Korn inequality in $BMO$. The uniqueness of equilibrium for a finite deformation whose principal stresses are everywhere nonnegative is then considered. It is shown that when the second variation of the energy, when considered as a function of the strain, is uniformly positive definite at such an equilibrium solution, then there is a $BMO$-neighborhood in strain space where there are no other equilibrium solutions.

Credits | Cookie policy | HTML 5 | CSS 2.1


text:
# 7 Character

A third fundamental type of data is character data (also called string data). In R, character vectors may be used as data for analysis, and also as names of data objects or elements. (R is its own macro language, and we use the same functions to manipulate language elements as we use to manipulate data values for analysis. In R, it is all “data”.)

As data for analysis, character values can signify categories (see Chapter 9). An example might be a variable that classifies people as “Democrat”, “Green”, “Independent”, “Libertarian”, or “Republican” (American political affiliations).

affiliations <- c("Dem", "Dem", "Rep", "Rep", "Ind", "Lib")
table(affiliations)
affiliations
Dem Ind Lib Rep
2   1   1   2

A single character value might also represent multiple categorical variables. The FIPS code “55025” is a combination of state (“55” for Wisconsin) and county (“025” for Dane) codes. And the date “2020-12-23” is a combination of a year, a month, and a day code (see Chapter 8).

You will want to be able to combine character values into a single value, and to separate a single value into parts.

Another aspect of working with character data is that they may represent the raw input from multiple people. If you have ever used social media you will appreciate that people’s views of acceptable capitalization, spelling, and punctuation vary enormously. Cleaning raw data is another important part of working with character data.

## 7.1 Combining Character Values

One basic task when working with character data is to combine elements from two or more vectors. This is useful whenever you need to construct a single variable to represent a value identified by multiple other variables. For example you might have data about calendar dates given as separate month, day, and year variables. To combine these into a single vector, use the paste() function (see help(paste)).

month <- c("Apr", "Dec", "Jan")
day   <- c(3, 13, 23)
year  <- c(2001, 2009, 1997)

date_str <- paste(year, month, day, sep="-")
date_str
[1] "2001-Apr-3"  "2009-Dec-13" "1997-Jan-23"

The paste() operation is vectorized in much the same way that numeric operations are. Notice that the results are character values.

The sep argument specifies a character value to place between the data elements being combined. The default separator is a space. To have nothing added between the elements being combined, we can either specify a null string, sep="" (quotes with NO space between), or we can use the paste0() function.

You might also use this if you were constructing a set of variable names with a common prefix. Notice the recycling in this example.

paste("Q", 1:4, sep="")
[1] "Q1" "Q2" "Q3" "Q4"
paste0("Q", 3, c("a", "b", "c"))
[1] "Q3a" "Q3b" "Q3c"

paste() and paste0() recycles each argument so that it matches the length of the longest argument, and then it concatenates element-wise. In the paste() statement above, the longest argument (1:4) is four elements long, so all others (here, just "Q") are recycled to length four (c("Q", "Q", "Q", "Q")). In the paste0() statement, the longest argument (c("a", "b", "c")) has three elements, so the others ("Q" and 3) are recycled until they are three elements long (c("Q", "Q", "Q") and c(3, 3, 3)). Then, they are concatenated element-by-element (the first element of each vector, the second element of each, and so on). Note that paste() will recycle an argument a non-whole number of times without a warning. Try paste0(c("a", "b", "c"), 1:2, "z") and notice how 1:2 is recycled to c(1, 2, 1) to have a length of three.

## 7.2 Working Within Character Values

A character value is an ordered collection of characters drawn from some alphabet. R is capable of working within a “local” alphabet, converting locales, or working in Unicode (a universal alphabet). The details of switching alphabets gets complicated quickly, so we will skip that here.

The most basic manipulations of character data values are selecting specific characters in a value (matching), removing selected characters, or adding characters.

Matching can be done either by position within a value, or by character.
In the character value “12:08pm” we could operate on the fourth and fifth characters (to find “08”), or we can specify that we want to operate on the character pair “08” (finding the fourth position). We are either looking for an arbitrary character that occupies a specific position, or we are looking for an arbitrary position occupied by a specific character.

## 7.3 Position Indexing

The substr() and substring() functions use positions and return characters values. The regexpr() function matches characters and returns starting positions (it is an index function).

x <- c("12:08pm", "12:10pm")
substr(x, start=4, stop=5)
[1] "08" "10"
regexpr("08", x)
[1]  4 -1
attr(,"match.length")
[1]  2 -1
attr(,"index.type")
[1] "chars"
attr(,"useBytes")
[1] TRUE

Here the character string sought is found at the fourth position in the first value, and not at all (-1) in the second value. (Everything after the first line of output is metadata, which makes this output hard to read.)

Although this example of regexpr() is very simple, a word of warning. Character matching functions is R rely on regular expressions, a system of specifying character patterns that includes literal characters (as in the example above), wildcards, and position anchors. We’ll come back to this, below.

To work with positions, if is often useful to know the length of a character value, for which we have the nchar() function.

If we wanted the last two characters of each character value, we could specify

x <- c("12:08pm", "12:10pm", "1:08pm")
x_len <- nchar(x)
substr(x, start=x_len-1, stop=x_len)
[1] "pm" "pm" "pm"

By using substr() on the left-hand side of the assignment operator, we can substitute in new substrings by position.

substr(x, start=x_len-1, stop=x_len) <- "am"
x
[1] "12:08am" "12:10am" "1:08am"

## 7.4 Character Matching

While some character value manipulations are easily handled by position indexing, many others are handled more gracefully through character matching.

### 7.4.1 Global Wildcards

You may already be familiar with the concept of wildcards to specify patterns - computer operating systems all allow wildcards in searching for file names. These are sometimes referred to as global or glob wildcards.

For example, on a Windows computer you could open a Command window and type the following command

dir /b 0*.rmd

to get a list of all the Rmarkdown (“rmd”) files in your current folder beginning with the character “0” (if there are any).

02_data_types.rmd
03_data_structures.rmd
04_data_class.Rmd
05_numeric.Rmd
06_logical.rmd
07_character.RMD
08_dates.RMD
09_categorical.rmd

The asterisk (*) wildcard matches any characters (zero or more) after the literal “0” at the beginning of file names, while the “.rmd” literal matches only files which end with that file extension.

Similarly, a question mark matches a single arbitrary character.

(In RStudio, you could open the Terminal tab, next to the Console, for a Unix-like shell. Drop the “/b”.)

With global wildcards, the pattern to match always specifies a string from beginning to end. So

• "02\*" matches any character string beginning with “02”
• "\*02" matches any string ending with “02”
• "\*02\*" matches any string containing “02”.

R does not use global wildcards directly, but the glob2rx function can translate this type of wildcard into a regular expression for you.

glob2rx("02*")
[1] "^02"
glob2rx("*02")
[1] "^.*02$" glob2rx("*02*") [1] "^.*02" ### 7.4.2 Regular Expression Wildcards Regular expressions expand on the concept of wildcards, and allow us to match elements of arbitrary character vectors with more precise patterns than global wildcards. We expand the concept of a wildcard by separating • what characters to match • how many characters to match • where to match (what position) A single arbitrary character is specified as a period, “.”, much like the global question mark, “?”. For example, one way to get a vector of column names that are at least four letters long using a regular expression would be cars <- mtcars grep("....", names(cars), value=TRUE) [1] "disp" "drat" "qsec" "gear" "carb" The grep() function searches a character vector for elements that match a pattern. It returns position indexes by default, or values that contain a match with the value=TRUE argument. The grepl() (grep logical) function returns a logical vector indicating which elements matched. These two functions give us all three methods of specifying indexes along a vector. In addition to wildcard characters, we can also match literal characters, and literal substrings. grep("a", names(cars), value=TRUE) [1] "drat" "am" "gear" "carb" grep("ar", names(cars), value=TRUE) [1] "gear" "carb" #### 7.4.2.1 Position In contrast to global wildcards, these patterns match anywhere within a character value - they are position-less. To specify positions We have two regular expression anchors that we can specify - tying a pattern to the beginning (“^”) or the end (“$”) of a string.

grep("m", names(cars), value=TRUE)  # any m
[1] "mpg" "am"
grep("^m", names(cars), value=TRUE) # begins with m
[1] "mpg"
grep("m$", names(cars), value=TRUE) # ends with m [1] "am" Although we have added position qualifiers to our patterns, notice that we are still specifying partial strings, not whole strings. To specify a complete string, we use both anchors! One way to find column names that are exactly two characters long would be grep("^..$", names(cars), value=TRUE)
[1] "hp" "wt" "vs" "am"

Without both anchors, this example would find all column names at least two characters long, including those with three and four characters.

#### 7.4.2.2 Repetition

So far we have been specifying one character at a time, but the regular expression syntax also includes the concept of repetition. There are six ways to specify how many matches are required:

• a question mark, “?”, matches zero or one time, making a character specification optional
• an asterisk, "*", matches zero or more times, a character is optional but also may be repeated
• a plus, “+”, matches one or more times, a character is required and may be repeated
• braces with a number, “{n}” matches exactly n times
• braces with two numbers, “{n,m}”, matches at least n times and no more than m times
• no repetition qualifier means match exactly once

So another way to get two-letter column names would be to specify

grep("^.{2}$", names(cars)) [1] 4 6 8 9 While the global wildcard “?” is replaced by the dot in regular expressions, the global wildcard “" is replaced with the regular expression ".”. #### 7.4.2.3 Character Class So far we have introduced arbitrary matches and literal matches, but regular expressions are able to work between these two extremes as well. We can specify classes (sets) of characters to match, and we can do this by itemizing the whole class, or using a shortcut name for some classes. Square brackets, “[ ]”, are used to itemize classes, and to specify shortcut names. As an arbitrary example, column names that begin with “a” or “b” or “c” could be specified grep("^[abc]", names(cars), value=TRUE) [1] "cyl" "am" "carb" Notice that this interacts with with the repetition qualifiers. To require the first two characters to belong to the same character class we would specify grep("^[abc]{2}", names(cars), value=TRUE) [1] "carb" The twelve shortcut names that are predefined in R are documented on the help("regex") page, and they include • [:alpha:], alphabetic characters • [:digit:], numerals • [:punct:], punctuation marks These shortcuts are specific to the use of regular expressions in R, and must themselves be used within class brackets. Contrast the first (correct) example with the second (incorrect) example. Why does the second line pick out “def”? grep("[[:digit:]]", c("abc", "123", "def"), value=TRUE) [1] "123" grep("[:digit:]", c("abc", "123", "def"), value=TRUE) [1] "def" #### 7.4.2.4 Regular Expression Metacharacters We can match arbitrary characters, specified classes of characters, and most literal characters. However, we are using some characters as metacharacters with special meaning in regular expression patterns. The dot (period, .), asterisk (*), question mark (?), plus sign (+), caret (^), dollar sign ($), square brackets ([, ]), braces ({, }), dash (-), and a few more we haven’t discussed, all have a non-literal meaning. What if we want to use these as literal characters?

There are generally two ways to take a metacharacter and use it as a literal. We can specify it within a square bracket class, or we can escape it. Either method comes with caveats.

To “escape” a character - ignore it’s special meaning and use it as a literal - we typically think of preceding it with a backslash, " \ ". However, it turns out that a backslash is also a regular expression metacharacter (that we have not discussed so far), so to use it as an escape character in regular expressions we double it. That is, to use an escape character, we need to first escape the escape character!

For example, to find a literal dollar sign contrast the correct specification with two mistakes.

grep("\\$", c("$12.95", "40502"), value=TRUE) # correct
grep("$", c("$12.95", "40502"), value=TRUE)  # wrong: no slash = end of string
grep("\$", c("$12.95", "40502"), value=TRUE) # wrong: error
Error: '\$' is an unrecognized escape in character string starting ""\$"

An alternative is to write (most) metacharacters within a class.

grep("[$]", c("$12.95", "40502"), value=TRUE) # correct
[1] "$12.95" The caveat here is that the caret, dash, and backslash all have special meaning within character classes. (A third approach is to turn off regular expression matching and use only literal matching. Use the fixed=TRUE argument.) #### 7.4.2.5 More There is more to regular expression specification: negation, alternation, substrings, and other less used elements. Whole books have been written about regular expressions. To learn more, a useful reference is Regular Expressions Info ## 7.5 Substitution Simply identifying matching values is useful for creating indicator variables (grepl) or for creating sets of variable names (grep). But for data values, we often want to manipulate the values we identify. One of the main tools we will use for this is substitution (sub) and repeated substitution (gsub). Substitution where a regular expression pattern appears at most once in each value is straightforward. Returning to the time example, which we previously solved by positional substitution, we can use the very simple regular expression “pm” to identify matching characters to replace with “am”. x <- c("12:08pm", "12:10pm", "1:08pm") sub("pm", "am", x) [1] "12:08am" "12:10am" "1:08am"  Substitution also works as a method of deleting matched characters, when the replacement is the null string (quotes with no space). ## 7.6 Exercises 1. Percent to fraction Given a character vector with values in percent form, convert these to numerical proportions, values between 0 and 1. x <- sprintf("%4.2f%%", runif(5)*100) x [1] "60.00%" "97.21%" "87.26%" "74.93%" "26.74%" 2. Currency Currency is sometimes denoted with both a currency symbol and commas. Convert these to numeric values. x <- c("$10", "$11.99", "$1,011.01")
3. Grouping values

There are really just three colors in this vector.

colors <- c("red", "rot", "rouge",
"blue", "blau", "bleue", "bleu",
"green")

Recode these data into three color categories.

4. Capitalization

Inconsistent capitalization is a problem with some alphabets.

colors2 <- c("Red", "blue", "red", "blue", "RED")

A quick solution is often to make everything either lower case or upper case. R has dedicated functions for this, see help("tolower") and fix this vector.

5. Translating wildcard patterns

The glob2rx() function translates character strings with wildcards (* for any string, ? for a single character) into regular expressions. We can translate "Merc *" (a string starting with “Merc” and a space, followed by anything) into “^Merc”. Combining this with grep() allows us to select rows from our mtcars dataframe. (Note that value = TRUE returns values, while the default value = FALSE returns positions.)

glob2rx("Merc *")
[1] "^Merc "
grep(glob2rx("Merc *"), row.names(mtcars), value=TRUE)
[1] "Merc 240D"   "Merc 230"    "Merc 280"    "Merc 280C"   "Merc 450SE"
[6] "Merc 450SL"  "Merc 450SLC"
grep(glob2rx("Merc *"), row.names(mtcars))
[1]  8  9 10 11 12 13 14
mtcars[grep(glob2rx("Merc *"), row.names(mtcars)), ]
             mpg cyl  disp  hp drat   wt qsec vs am gear carb
Merc 240D   24.4   4 146.7  62 3.69 3.19 20.0  1  0    4    2
Merc 230    22.8   4 140.8  95 3.92 3.15 22.9  1  0    4    2
Merc 280    19.2   6 167.6 123 3.92 3.44 18.3  1  0    4    4
Merc 280C   17.8   6 167.6 123 3.92 3.44 18.9  1  0    4    4
Merc 450SE  16.4   8 275.8 180 3.07 4.07 17.4  0  0    3    3
Merc 450SL  17.3   8 275.8 180 3.07 3.73 17.6  0  0    3    3
Merc 450SLC 15.2   8 275.8 180 3.07 3.78 18.0  0  0    3    3

Now, try selecting rows from mtcars where the row name… * starts with “Toyota” and a space * starts with any four characters and then a space * ends with 0 * ends with a space and then any three characters


text:
# Dirichlet Processes: the basics

Published:

The Dirichlet process (DP) is one of the most common – and one of the most simple – prior distributions used in Bayesian nonparametric models. In this post, we’ll review a couple different interpretations of DPs.

## Dirichlet process

The Dirichlet process is a stochastic process wherein any arbitrary partition of the probability space has a Dirichlet distribution. One can think of it as a probability distribution that forms a prior for other probability distributions. In other words, a sample from a DP, $F \sim DP$, will be a probability distribution. For this reason, DPs can be used as a prior when estimating a CDF.

A DP, denoted $DP(\alpha, F_0)$ has two parameters: $F_0$ is a prior guess at the distribution in question, and $\alpha$ is a dispersion parameter that dictates how closely samples follow $F_0$.

To sample from a DP, it’s common to use a “stick-breaking” process. We present the generative process here, and explain the intuition below.

Stick-breaking process for sampling $F \sim DP(\alpha, F_0)$:

1. Draw $s_1, \dots, s_\infty$ from $F_0$.
2. Draw $V_1, \dots, V_\infty$ from $\text{Beta}(1, \alpha)$.
3. Let $w_1 = V_1$ and $w_j = V_j \prod\limits_{i = 1}^{j - 1} (1 - V_i)$ for $j = 2, \dots \infty$.
4. Let $F$ be the discrete distribution that puts mass $w_j$ at $s_j$, that is, $F = \sum\limits_{j = 1}^\infty w_j \delta_{s_j}$ where $\delta_{s_j}$ is a point mass at $s_j$.

The intuition is as follows: Imagine a stick with length $1$. Then, iteratively, we break the stick. On iteration $1$, we sample a value $V_1 \in [0, 1]$ from a Beta distribution, and break the stick at $V_1$. Now our stick has length $1 - V_1$. On iteration $2$, we sample another value $V_2 \in [0, 1]$, and break off $V_2$ fraction of the stick, leaving us with a stick of size $V_2(1 - V_1)$. We continue this process, letting the stick get smaller and smaller. The size of the stick at each iteration determines the weight placed on each point in the distribution.

Notice that the DP places more weight on samples that are early in the list of samples, and also places more weight on samples that are favored by the particular Beta distribution.

In practice, because we can’t draw infinite samples from $F_0$ in step 1, to demonstrate sampling from a DP, we truncate this to drawing $k$ samples.

We show how to perform this process in Python below, setting $F_0 = \mathcal{N}(0, 1)$ and $\alpha = 10$:

import numpy as np
from scipy.stats import beta
from scipy.stats import norm

k = 100
alpha = 10

# F_0
s = np.random.normal(0, 1, size=k)
# Beta for drawing V's
v = np.random.beta(1, alpha, size=k)

ws = [v[0]]
for ii in range(1, k):
ws.append(v[ii] * np.prod(1 - v[:ii]))
ws = np.array(ws)

Looking at the resulting draws, we can see that it starts to resemble $F_0$, but as a discrete distribution. Here we plot the draws and their CDF.

Drawing multiple samples we can start to get a sense for how much variability this prior (and our current hyperparameter settings) are allowing:

Finally, we can test how different values for $\alpha$ affect the variability of our DP draws:

Tags: