text:
# Model Comparison: How to constrain a model to positive mu-differences?

Hi! I’m new to PyMC3 and bayesian modelling in general. I’m not even sure my question is appropiate to my problem. I’m not sure about my modelling choices. So please be a bit patient with me, thanks!
There’s a code example following the explanations.

Background:
I’ve got data from an experiment with a repeated measures/within-subject design (the real data is from a mixed design, though the between factor doesn’t concern me right now). Each participant does 3 consecutive blocks of 30 trials. Due to outlier removal and missing data the resulting blocks do not have the same length. The data are lengths of projections of 2D points onto a vector (called parallel in the data) and onto its orthogonal vector. I’m interested in the difference in variance between these projection directions parallel (\parallel) and orthogonal (\perp).

Goal:
I’ve done the classic ANOVA approach (on synergy indices as a measure of the ratio between projections in an uncontrolled manifold analysis, for those who are interested), but I’d like to learn from the data how plausible certain theoretical models about the structure of variance are, given the observed data. My goal is to receive this kind of comparison between models. I used @junpenglao’s bridge sampling implementation, but didn’t parameterize the models very well.

Approach:
The mean of the projections is zero in each block. Therefore, the mean of the squared projections is the variance. So I could compare mean squared projections. Squared projections are all positive values, but there are a few extremes, the variance seems to increase with the mean, the squared data are skewed to the right, and the range of values is large. These were indicators to me that I should log-transform the squared data.

Here I’ll start with 2 simple models: The null model which proposes all data is generated by the same distribution. And a model for a main effect of projection direction (parallel > orthogonal). To come back to my question: How do I make sure that for the main effect model the estimated \mu_\parallel > \mu_\perp ?

The difference in the means of transformed projections between participants is so large that there’s overlap in my priors for \mu_\parallel & \mu_\perp . So it could happen that the model generates 2 almost identical \mu_\parallel & \mu_\perp, when in theory I want the likelihood of the data, given that \mu_\parallel > \mu_\perp, and subsequently P(\mu_\parallel > \mu_\perp|D) vs. P(\mu_\parallel = \mu_\perp|D).
In my mind I would need to standardize the data somehow, so I can make the probabilities for \mu_\parallel and \mu_\perp non-overlapping? Or somehow model the difference between the 2 to be at least > 0? You see, I’m a bit confused about it.

# %%
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import seaborn as sns

# %%
# Create some data conforming to null model: all deviations on average equal across blocks.
n_trials = 90 - np.random.randint(36)
df0 = pd.DataFrame({'user': 0, 'block': 1 + np.random.randint(3, size=n_trials),
'parallel': np.random.normal(0.0, 4.0, size=n_trials),
'orthogonal': np.random.normal(0.0, 4.0, size=n_trials)})
# Make sure mean of each block is zero as with real data.
df0[['parallel', 'orthogonal']] = df0.groupby('block')[['parallel', 'orthogonal']].transform(lambda x: x - x.mean())

# %%
# Create some data conforming to main effect of projection model: variance in parallel > orthogonal.
n_trials = 90 - np.random.randint(36)
df1 = pd.DataFrame({'user': 1, 'block': 1 + np.random.randint(3, size=n_trials),
'parallel': np.random.normal(0.0, 5.0, size=n_trials),
'orthogonal': np.random.normal(0.0, 1.0, size=n_trials)})
# Make sure mean of each block is zero as with real data.
df1[['parallel', 'orthogonal']] = df1.groupby('block')[['parallel', 'orthogonal']].transform(lambda x: x - x.mean())

# %%
df = pd.concat((df0, df1), axis='index')

# %%
df[['user', 'block']] = df[['user', 'block']].astype('category')
df.sort_values(by=['user', 'block'], inplace=True)

# log-transform squared deviations.
df[['parallel', 'orthogonal']] = df[['parallel', 'orthogonal']].transform('square').transform('log')

# %%
# Vizualize the tranformed data.
samples = df.melt(id_vars='user', value_vars=['parallel', 'orthogonal'], var_name='projection', value_name='ln(proj_sq)')
g = sns.FacetGrid(samples, col='user', hue="projection")
g.despine(offset=10)
g = (g.map(sns.distplot, "ln(proj_sq)").add_legend(title="Projection", loc='upper right'))
for ax in g.axes.flat:
ax.set_ylabel("Probability Density")
ax.set_xlabel("Ln(Squared Projection length)")
plt.tight_layout()
plt.show()

# %%
# Massage data into single distributions per user. Each column holds all data for 1 user.
# Since not all users have the same amount of data, we will have missing data at the ends of the columns.
null_data = df.melt(id_vars='user', value_vars=['parallel', 'orthogonal'], var_name='projection', value_name='ln(proj_sq)').drop('projection', axis='columns')
null_data['idx'] = null_data.groupby('user').cumcount()
null_data = null_data.pivot(index='idx', columns='user', values='ln(proj_sq)')
n_users = df['user'].nunique()

# %%
with pm.Model() as null_model:
null_model.name = "Null Model"
mu = pm.Normal("mu", mu=0.8, sigma=1.25, shape=n_users)
sigma = pm.Gamma("sigma", alpha=25.5, beta=9, shape=n_users)

obs = pm.Normal(f"obs", mu=mu, sigma=sigma, shape=df['user'].nunique(), observed=null_data)

# %%
with null_model:
idata0 = pm.sample(1000, tune=800, return_inferencedata=True)

# %%
az.summary(idata0, var_names=["Null Model_mu", "Null Model_sigma"])

# %%
# Massage data into single distributions per user and projection direction.
proj_effect_data = df.loc[:, ['user', 'parallel', 'orthogonal']]
proj_effect_data['idx'] = proj_effect_data.groupby('user').cumcount()
proj_effect_data = proj_effect_data.pivot(index='idx', columns='user', values=['parallel', 'orthogonal'])

# %%
with pm.Model() as projection_model:
projection_model.name = "Main Effect Projection"
mu_ortho = pm.Normal("mu_ortho", mu=0.09, sigma=2.55, shape=n_users)
mu_parallel = pm.Normal("mu_parallel", mu=1.5, sigma=2.5, shape=n_users)
sigma = pm.Gamma("sigma", mu=2.65, sigma=0.25, shape=n_users)

obs_parallel = pm.Normal(f"obs_parallel", mu=mu_parallel, sigma=sigma, shape=n_users,
observed=proj_effect_data['parallel'])
obs_orthogonal = pm.Normal(f"obs_orthogonal", mu=mu_ortho, sigma=sigma, shape=n_users,
observed=proj_effect_data['orthogonal'])

diff = pm.Deterministic('mu_diff', mu_parallel - mu_ortho)

# %%
with projection_model:
idata1 = pm.sample(1000, tune=800, return_inferencedata=True)

# %%
az.summary(idata1, var_names=["Main Effect Projection_mu", "Main Effect Projection_sigma"], filter_vars='like')

There are more models I want to evaluate for main effects of block and interactions (e.g. only parallel projection variance is smaller in block 2), so I need to understand the basics.

Previous approach:
Before I decided to log-transform the data, \mu_\perp and a \Delta\mu were modelled using gamma distributions and obs_parallel had \mu_\perp + \Delta\mu as its mu-parameter to construct \mu_\parallel.

I also created the models for each participant to feed the result to the bridge sampling. Then I found out I could vectorize the model. So now I would need to figure out, if and how the bridge sampling works with the vectorized form. The whole project can be found here.

1 Like

You can model the differences as positive value, and add it to mu_ortho:

with pm.Model() ...
...
mu_ortho = pm.Normal("mu_ortho", mu=0.09, sigma=2.55, shape=n_users)
mu_diff = pm.HalfNormal("mu_diff", 1., shape=n_users)
mu_parallel = pm.Deterministic("mu_parallel", mu_ortho+mu_diff)

However, instead of computing the variance of the 2 conditions and log transforming them, I would prefer modeling the single trial - in that case, you can give informative prior to the variance, and works much better for single subject modeling. In addition, I wouldnt even restrict parallel > orthogonal, and just look at the posterior samples count how many parallel > orthogonal

Also, for model comparison, generally LOO-CV based method is better.

2 Likes

Thanks for your help! I believe your answer is the solution to the question I had. Thanks for your recommendations on the modelling as well. As I mentioned I wasn’t sure about the model choices and methods to use. If you don’t mind, I have some follow-up questions.

I’m not sure what you mean by modeling single trials. Having a probability distribution for each single trial? The observed data I produce in the example below (and above) contains up to 90 trials per participant (trial count differs per block and participant). My previous code was transforming the single trials so their mean would give the variance.
Having a random variable for the standard deviation would give me what I need, right? Or did you mean to rather do the modelling for the untransformed projections? In that case I changed the example to the one below with untransformed projections. Though I’m not sure a normal distribution is the right fit for the real data. I think a Shapiro-Wilk test indicated non-normality.

short version:

with pm.Model() as projection_model:
projection_model.name = "Main Effect Projection"
sigma_ortho = pm.Gamma("sigma_ortho", mu=3.34, sigma=1.44, shape=n_users)
sigma_parallel = pm.Gamma("sigma_parallel", mu=7.0, sigma=3.17, shape=n_users)
# Alternatively: The sum of two gamma-distributed variables is gamma-distributed.
#sigma_diff = pm.Gamma('sigma_diff', mu=4.0, sigma=5.0, shape=n_users)
#sigma_parallel = pm.Deterministic("sigma_parallel", sigma_ortho + sigma_diff)

obs_parallel = pm.Normal(f"obs_parallel", mu=0.0, sigma=sigma_parallel, shape=n_users, observed=proj_effect_data['parallel'])
obs_orthogonal = pm.Normal(f"obs_orthogonal", mu=0.0, sigma=sigma_ortho, shape=n_users, observed=proj_effect_data['orthogonal'])
sigma_diff = pm.Deterministic('sigma_diff', sigma_parallel - sigma_ortho)

full example:

# %%
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import seaborn as sns

# %%
# Create some data conforming to null model: all deviations on average equal across blocks.
n_trials = 90 - np.random.randint(36)
df0 = pd.DataFrame({'user': 0, 'block': 1 + np.random.randint(3, size=n_trials),
'parallel': np.random.normal(0.0, 5.5, size=n_trials),
'orthogonal': np.random.normal(0.0, 5.5, size=n_trials)})
# Make sure mean of each block is zero as with real data.
df0[['parallel', 'orthogonal']] = df0.groupby('block')[['parallel', 'orthogonal']].transform(lambda x: x - x.mean())

# %%
# Create some data conforming to main effect of projection model: variance in parallel > orthogonal.
n_trials = 90 - np.random.randint(36)
df1 = pd.DataFrame({'user': 1, 'block': 1 + np.random.randint(3, size=n_trials),
'parallel': np.random.normal(0.0, 7.0, size=n_trials),
'orthogonal': np.random.normal(0.0, 3.3, size=n_trials)})
# Make sure mean of each block is zero as with real data.
df1[['parallel', 'orthogonal']] = df1.groupby('block')[['parallel', 'orthogonal']].transform(lambda x: x - x.mean())

# %%
df = pd.concat((df0, df1), axis='index')

# %%
df[['user', 'block']] = df[['user', 'block']].astype('category')
df.sort_values(by=['user', 'block'], inplace=True)

# %%
# Vizualize the tranformed data.
samples = df.melt(id_vars='user', value_vars=['parallel', 'orthogonal'], var_name='direction', value_name='projection')
g = sns.FacetGrid(samples, col='user', hue="direction")
g.despine(offset=10)
g = (g.map(sns.distplot, "projection").add_legend(title="Projection", loc='upper right'))
for ax in g.axes.flat:
ax.set_ylabel("Probability Density")
ax.set_xlabel("Projection Length")
plt.tight_layout()
plt.show()

# %%
# Massage data into single distributions per user. Each column holds all data for 1 user.
# Since not all users have the same amount of data, we will have missing data at the ends of the columns.
null_data = df.melt(id_vars='user', value_vars=['parallel', 'orthogonal'], var_name='direction', value_name='projection').drop('direction', axis='columns')
null_data['idx'] = null_data.groupby('user').cumcount()
null_data = null_data.pivot(index='idx', columns='user', values='projection')
n_users = df['user'].nunique()

# %%
with pm.Model() as null_model:
null_model.name = "Null Model"
sigma = pm.Gamma("sigma", alpha=2.2, beta=0.4, shape=n_users)

obs = pm.Normal(f"obs", mu=0.0, sigma=sigma, shape=n_users, observed=null_data)

# %%
with null_model:
idata0 = pm.sample(1000, tune=800, return_inferencedata=True)

# %%
az.summary(idata0, var_names=["Null Model_sigma"])

# %%
# Massage data into single distributions per user and projection direction.
proj_effect_data = df.loc[:, ['user', 'parallel', 'orthogonal']]
proj_effect_data['idx'] = proj_effect_data.groupby('user').cumcount()
proj_effect_data = proj_effect_data.pivot(index='idx', columns='user', values=['parallel', 'orthogonal'])

# %%
with pm.Model() as projection_model:
projection_model.name = "Main Effect Projection"
sigma_ortho = pm.Gamma("sigma_ortho", mu=3.34, sigma=1.44, shape=n_users)
sigma_parallel = pm.Gamma("sigma_parallel", mu=7.0, sigma=3.17, shape=n_users)
# Alternatively: The sum of two gamma-distributed variables is gamma-distributed.
#sigma_diff = pm.Gamma('sigma_diff', mu=4.0, sigma=5.0, shape=n_users)
#sigma_parallel = pm.Deterministic("sigma_parallel", sigma_ortho + sigma_diff)

obs_parallel = pm.Normal(f"obs_parallel", mu=0.0, sigma=sigma_parallel, shape=n_users, observed=proj_effect_data['parallel'])
obs_orthogonal = pm.Normal(f"obs_orthogonal", mu=0.0, sigma=sigma_ortho, shape=n_users, observed=proj_effect_data['orthogonal'])
sigma_diff = pm.Deterministic('sigma_diff', sigma_parallel - sigma_ortho)

# %%
with projection_model:
idata1 = pm.sample(1000, tune=800, return_inferencedata=True)

# %%
az.summary(idata1, var_names=["Main Effect Projection_sigma"], filter_vars='like')

My final goal is to compare the marginal likelihoods of several different theoretically driven models of directed hypotheses (null, main effect parallel > orthogonal, main effect of block 2 on both directions, different interaction effects, …)
I don’t yet understand how the parallel > orthogonal counts would lead me to the marginal likelihoods of the models.
I was previously using your bridge sampling implementation this kind of way (not exactly):

def get_marginal_likelihood(model, trace):
with model:
marginal_log_likelihood = bridgesampler.marginal_llk(trace, model=model, maxiter=10000)["logml"]
return marginal_log_likelihood

# models and traces for each particpant, I hadn't vectorized it yet at the time.
users = [1,2,3,4,5]  # and so on.
model_posterior = dict()
for user in users:
model_posterior[user] = np.array([get_marginal_likelihood(*model) for model in zip(models[user], traces[user])])  # models[user] would contain all (unvectorized) models for this user.
# Exponentiate and normalize to turn into a probability distribution.

So you’re saying using a LOO-CV based method would yield more accurate results?
I tried az.compare({null_model.name: idata0, projection_model.name: idata1}), but it fails with:

TypeError: Found several log likelihood arrays ['Projection Effect_obs_parallel', 'Projection Effect_obs_orthogonal'], var_name cannot be None

I read about the error in this answer, but I don’t understand it enough to translate it to my application.

1 Like

I think I have to clarify some more why I turned to pymc3 and I’m not just doing a Bayesian ANOVA in JASP or jamovi.
I’m especially interested in predictions from 2 hypotheses I haven’t mentioned before.
1: low variability for orthogonal in all blocks, high variability for parallel projection in block 1 & 3, but not in block 2.

sigma
^
| .       . parallel
|  \     /
|   \   /
|    \./
| .___.___. orthogonal
|------------> block
1   2   3

2: same as before except that orthogonal variability rises in block 2.

sigma
^
|.       . parallel
| \     /
|  \   /
|   \./
|   /*\
| ./   \.  orthogonal
|------------> block
1  2  3

The ANOVA on the Synergy Index (a measure of the variability ratio) can’t discriminate these 2 predictions, only that there is some difference in block 2 vs 1 & 3, which doesn’t help much. Therefore, I hope to model these 2 predictions in pymc3, amongst other ones.
That’s why I think counting the samples where parallel > orthogonal would not be able to give me the answer I seek, because it would give me the same result for these 2 examples?

1 Like

Ah I see, the observed is a square array (n_subject * n_trials). In that case, my suggestion would be to flatten the observed and use indexing to construct the linear predictors, it is better for missing trials etc - something like:

with pm.Model() as projection_model:
projection_model.name = "Main Effect Projection"
# better to also model mu
mu_ortho = pm.Normal(...)
mu_parallel = pm.Normal(...)
sigma_ortho = pm.Gamma("sigma_ortho", mu=3.34, sigma=1.44, shape=n_users)
sigma_parallel = pm.Gamma("sigma_parallel", mu=7.0, sigma=3.17, shape=n_users)

obs_parallel = pm.Normal(f"obs_parallel",
mu=mu_parallel[parallel_idx],
sigma=sigma_parallel[parallel_idx],
observed=parallel_data_flatten)
obs_orthogonal = pm.Normal(f"obs_orthogonal",
mu=mu_ortho[orthogonal_idx],
sigma=sigma_ortho[orthogonal_idx],
observed=orthogonal_data_flatten)
sigma_diff = pm.Deterministic('sigma_diff', sigma_parallel - sigma_ortho)

This is reaction time data? maybe try log-normal? I am not sure what is the recommended best practice nowadays as well.

Well, the language of hypothesis testing is always problematic to me as it usually dont answer what we really want to know - is this conclusion going to be useful (for making decision, for advancing our understanding, etc). As it is only useful when you have a non-tiny effect, looking at things like effect size is much more natural. In a linear model (like ANOVA), all the hypothesis you want to test are ultimately comparisons of parameters/coefficients, so instead of fitting multiple models and compare marginal likelihood or Bayes Factor, I usually recommend to fit the most comprehensive model with all the information (in this case would mean fit both conditions and coded for all subject, trial and block), and look at the distribution of the posterior contrasts (an example: https://docs.pymc.io/notebooks/BEST.html)

You need to coded the block in your model as well - for example, you can have a indexing of [0, 1, 2, …, 5] to indicate parallel_block_1, parallel_block_2, …, orthogonal_block_3, and index to a shape=6 vector of parameter.

2 Likes

Thank you very much for taking the time again to answer my questions. I try to understand it and why it may be the better approach.

I agree, that’s why I’ve come here. Too bad I’m still a noob at Bayesian inference.

No. It’s a motoric task measuring 2 end effector positions, thus forming a trial-to-trial point cloud in a cartesian 2D-coordinate system. The demeaned 2D-coordinate vectors get projected onto a vector and its orthogonal, basically just a change of basis. “parallel” and “orthogonal” data are merely the coefficients in \vec{v} = a_\parallel\hat{i} + a_\perp \hat{j}.

Yes, in my first post this was the case. That’s why I modelled mu which then was the parameter I was interested in, since in that case mu resembled the variance of the untransformed projection data. It was also log-transformed to pull in extremes (due to the squares) and stabilize the variance of mu itself (variance of variance) etc. It may have been more complex than it needs to be.
For the non-squared data in my second post, however, the mean is guaranteed to be 0, that’s why I didn’t model mu, only sigma.

By flattening you mean throwing all participants together? I don’t quite get your code example. mu and sigma are of shape n_users, but you only estimate these for 1 user and using the data of all users at the same time? I’m confused.
Since participants may deploy different strategies in the motoric task (assigned to different conditions as between factor and each participant may come up with their own solution), I want to estimate model parameters per participant. I started with having each single model per each participant before I discovered I could vectorize the parameters. Having all models per particpant eliminated missing data. Only after vectorizing parameters instead of having each single model for each participant I was facing the issue of missing data, since participants differ in their trial count per block (1 column per participant and trials being the rows in the observed data).

Yes, I only left out the blocks from my previous examples for brevity, since they weren’t relevant to those simple models, and for trying to understand the basics first. If I understand you correctly, this is the approach where I have 1 uber-model per participant? And for variables I’m interested in (9 differences) I would create deterministic variables,e.g. sigma[2] - sigma[1] for difference of parallel deviation block2-block1, and sigma[1]-sigma[4] for difference between parallel and orthogonal deviation in block 2?
But if I had all blocks in 1 vectorized form (trials as rows, block_<projection> as columns) I’d again have missing data, since the trial count differs between blocks? This can be solved though by having block1, block2, block3 with shape 2 each (parallel & orthogonal).

Having those 9 estimates and their uncertainty would be a good thing. Thanks!
Only, having all 9 estimates isn’t fully satisfying my question, since there are competing theoretical predictions about the sizes of those differences, see above ASCII diagrams. I’d like to know whose theory is most compatible with the observed data. Wouldn’t that involve a secondary step of a model selection with the estimated differences as observed data, including their uncertainty somehow?

The approach I had with several models was to state from the beginning what I’d expect the differences to look like for each theory. Null model would have no differences whatsoever. Main effect projection would state that there is a constant unknown difference (>0) between the projections across blocks. And then more sofisticated models with statements about the differences between individual blocks/projections. I feel like this was the other way around of what you suggested, first getting the estimates ofthose differences and then making statements about them.

1 Like

Since they are orthogonal you can model them with univariate distribution that has similar shape and the same support (e.g., positive only) as your observed.

Usually, you can coded your data as a table like:

subject_id | condition1 | condition2 | observation
0 |   parallel |    block_1 |       .3452
...

Instead of having 1 column per participant and trials being the rows in the observed data, you have 1 column for all your data, and participant and trials are indicator (subject_id and trial_id).

Then in your model, you can either use GLM like approach to construct a design matrix, or use indexing:

subject_mu = pm.Normal(..., shape=n_subject)  <= shape is n_subject
subject_mu_all_trial = subject_mu[subject_id] <= shape is now number of row in your data frame

which estimate the parameter independently for each subject.

No matter what approach you decide to go with in terms of model comparison / hypothesis testing, I think you should use this formulation as it is a more general way to do so, and connect nicely with classical GLM formulation. You can also consult the recent book Regression and other stories

2 Likes

Thanks for the clarification on the data format. So just have it in long format.

I’m sorry, I still don’t understand what type of variable subject_id in subject_mu_all_trial = subject_mu[subject_id] shall be in your example. Is that the single ID of a participant (e.g. single integer)? But then why have subject_mu with shape n_subjects, if I’d have to loop over all subject_id anyway and call with pm.Model() as model: ... within the loop? Or is subject_id an array of all unique subject IDs [0,1,2,3,…]?
Or is it the first column of the data with duplicates [0,0,0,0,0,…,1,1,1,1,1,…,2,2,2,2,2,2,…]? If the shape of subject_mu_all_trials was the number of rows in the long format dataframe (total number of all trials*2) it would be n_users*n_user_trials*n_directions (e.g. 40x90x2=7200)? Wouldn’t that just give me the values of every single trial? I obviously don’t understand what you are trying to tell me.

Or did you mean to do it in 2 steps:

df = ... # long format dataframe with column 'user' as ID for participants.
n_users = df['user'].nunique()
mus = list()
with pm.Model() as model:
user_mu = pm.Normal(..., shape=n_users)

for user_id in df['user'].unique():
with model:
user_mu_all_trials = user_mu[user_id]
mus.append(user_mu)  # Otherwise it would be overwritten each iteration.

That doesn’t seem to be what you intended. Thanks for your effort and patience with me!

1 Like

Yes

It wouldn’t as you dont have parameter related to each trial - this is essentially what IID data or repeat measurement means. Maybe it is clearer in code:

data1 = np.random.normal(loc=0., scale=1., size=10)
data2 = np.random.normal(loc=10., scale=5., size=10)
data_all = np.concatenate([data1, data2], axis=0)
sbj_idx = np.concatenate([np.zeros(10), np.ones(10)], axis=0).astype(int)

with pm.Model() as model1:
mu = pm.Normal('mu', 0., 10., shape=2)
sigma = pm.HalfNormal('sigma', 5., shape=2)
obs = pm.Normal('obs', mu, sigma, observed=np.stack([data1, data2]).T)

with pm.Model() as model2:
mu = pm.Normal('mu', 0., 10., shape=2)
sigma = pm.HalfNormal('sigma', 5., shape=2)
obs = pm.Normal('obs', mu[sbj_idx], sigma[sbj_idx], observed=data_all)

model1 and model2 are identical but model2 is much more flexible as it allows different length of data among subjects.

4 Likes

Thank you so much for the code example. This is great info. I don’t know if that is even anywhere in the documentation. Using the graphviz plot I see what’s going on. I imagine it now as mapping the observed data to the mus. I scoured the internet in an effort to find any explanation for how to use it, to not occupy too much more of your time, but I didn’t find a single example using the indexing method. So I am once again kindly asking for your help.

I found a well illustrated primer on different shape types. And a blog post which explains that shapes in PyMC3 are, unfortunately, one hot mess. It also states:

distribution_shape == batch_shape + event_shape
[observed] data should have a shape that looks like sample_shape + distribution_shape .

If I understood correctly, for a model where I estimate 6 parameters (e.g. \sigma for n_blocks(3)*n_projections(2)) for each user by drawing from independent univariate distributions (same distribution family), my sample shape is n_trials (per block), batch shape is n_users*n_blocks*n_projections, and event shape = ().

I have the data in long format like you suggested:

user block  direction    projection
0    1      parallel     7.129121
0    1      parallel     3.414878
0    1      parallel     -6.330880
0    1      parallel     0.428940
0    1      parallel     -3.650233
...  ...    ...          ...
1    3      orthogonal   -2.223337
1    3      orthogonal   -1.275683
1    3      orthogonal   -0.796934
1    3      orthogonal   2.221719
1    3      orthogonal   -5.626795

I could unstack the ‘direction’ if that was of any use, as these are guaranteed to be the same length for each individual block.
I tried different shape settings and I don’t know which one is correct. For example,

with pm.Model() as model3:
sigma = pm.HalfNormal('sigma', 5., shape=12)  # 2 users * 6 parameters. batch shape?
obs = pm.Normal('obs', 0.0, sigma[samples['user']], shape=12, observed=samples['projection'])

But it’s hard to identify which parameter (sigma 0-11) belongs to which user/block/direction coordinate.
Whatever else I tried with the indexing and shape argument, I always got this error:

RecursionError: maximum recursion depth exceeded in comparison

If I manage to use it as intended I may write something about how to vectorize a model properly without missing data using your indexing method to ease the struggle for other beginners.

1 Like

So, you can still use indexing to do this, but it will be quite inefficient as you need to create a index of direction + block, for example, index 0 would be parallel block 1.

A better approach is to formulate a design matrix:

from patsy import dmatrices
outcome, predictors = dmatrices("projection ~ block + direction + block * direction", data)

which is what usually call mixed effect model or hierarchical model. The pymc3 doc have quite an extensive session on this (e.g., https://docs.pymc.io/notebooks/GLM-hierarchical.html)

The differences is that, here you your case the linear prediction will goes into sigma instead of mu

2 Likes

Hi, we know cross validation and model comparison can be tricky because in many cases there are several correct ways of comparing the models. We are actually working on some resources to help with that, they are still a work in progress, but it should be useful already. Here is the notebook covering multiple likelihoods. The link is to my fork which has the open PR, feedback will be most appreciated. You may also be interested in the other notebooks of the section: https://nbviewer.jupyter.org/github/arviz-devs/Exploratory-Analysis-of-Bayesian-Models/tree/master/content/Section_04/

Having said that, there is also this example in R which I think will be much closer to your model and could also be helpful: https://avehtari.github.io/modelselection/rats_kcv.html

I will try to read this topic with more time to see if I can help further. I think this multiple measurement and multiple subject models are ideal to teach the nooks and crannies of CV.

By any chance, do you know about some public dataset other than the rats one with similar characteristics?

2 Likes

I needed some time to go through the ressources you both provided, dust off my knowledge about HLM (Mixed Model) and think about how this would apply to my actual problem. I’ll discuss now the 2 different approaches, one being mixed model, and the other I call the intuitive approach.

First I want to stress and illustrate again what my research question is. There are several hypothesis by different authors about how the variances of the data should differ between the blocks and direction of projection.
How the data is nested:

3 between-conditions
N users (N varies between conditions)
3 within-conditions: blocks
X trials (varies per block and user)
2 outcome directions: parallel, orthogonal
outcome: projection (scalar)

The question of interest is: Which prediction about the variance structure in the data is more probable?
To illustrate this I present again some ASCII-diagrams for different hypotheses. This first illustration makes it easy to think about the data in terms of its collection process (one block after the other, 2 outcomes). These only show block and direction as independent categorical variables and the standard devation as the depedent variableand neglect the between-factor and different users.

Main Effect Projection  Main Effect Block2    Main Effects Projection + Block2
sigma                   sigma                 sigma
^                       ^                     ^     o        o - parallel
|                       |                     |    / \       x - orthogonal
| o__o__o               |     o               |   /   \
|                       |    /x\              |  /  x  \
|                       |   // \\             | o  / \  o
|                       |  //   \\            |   /   \
| x__x__x               | //     \\           |  /     \
|                       |ox       ox          | x       x
|----------> block      |------------> block  |------------> block
1  2  3                 1   2   3            1   2   3

Block x Projection v1   Block x Proj. v2      Block x Projection v3
sigma                   sigma                 sigma
^                       ^                     ^
|o       o              | o       o           | o
| \     /               |  \     /            |  \
|  \   /                |   \   /             |   \
|   \o/                 |    \ /              |    \
|   /x\                 |     V               |     o___o
|  /   \                |     o               |
| x     x               | x___x___x           | x___x___x
|----------> block      |------------> block  |------------> block
1  2  3                 1   2   3             1   2   3

Block x Projection v4   Null Model
sigma                   sigma
^                       ^
|                       |
|                       |
| o__o__o               |
|    x                  |
|   / \                 | xo__xo__xo
|  /   \                |
| x     x               |
|----------> block      |------------> block
1  2  3                 1   2   3

Multiple users (e.g. 3) may look like this:

Block2 x Projection v1  Block2 x Proj. v2     Main Effects Block2 + Projection
sigma                   sigma                 sigma o
^                       ^                     ^
|                       |                     |     oo         o - parallel
|                       |         o           |     o          x - orthogonal
| oo      o             | oo      o           | oo      o
| o       oo            | o       o           |         oo
|     oo                |                     | o   xx
|     xox               |                     |     x
| x   x                 | x   oo  xx          | x
| xx      xxx           | xx  xxx x           | xx      xxx
|------------> block    |------------> block  |------------> block
1   2   3               1   2   3             1   2   3

These would be intuitive to formulate as separate models with restrictions on average sigma per block x direction.

So let’s look at how I would go about the intuitive approach. Here’s a code example with 3 models and 3 users behaving differently.

# %%
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc3 as pm
import seaborn as sns

#################################
#          Dummy Data           #
#################################
# %%
# Create some data conforming to null model: all deviations on average equal across blocks.
n_trials = 90 - np.random.randint(36)
df0 = pd.DataFrame({'user': 0, 'block': 1 + np.random.randint(3, size=n_trials),
'parallel': np.random.normal(0.0, 5.5, size=n_trials),
'orthogonal': np.random.normal(0.0, 5.5, size=n_trials)})
# Make sure mean of each block is zero as with real data.
df0[['parallel', 'orthogonal']] = df0.groupby('block')[['parallel', 'orthogonal']].transform(lambda x: x - x.mean())

# %%
# Create some data conforming to main effect of projection model: variance in parallel > orthogonal.
n_trials = 90 - np.random.randint(36)
df1 = pd.DataFrame({'user': 1, 'block': 1 + np.random.randint(3, size=n_trials),
'parallel': np.random.normal(0.0, 7.0, size=n_trials),
'orthogonal': np.random.normal(0.0, 3.3, size=n_trials)})
# Make sure mean of each block is zero as with real data.
df1[['parallel', 'orthogonal']] = df1.groupby('block')[['parallel', 'orthogonal']].transform(lambda x: x - x.mean())

# %%
# Create some data conforming to interaction effect version 1:
# blocks 1&3 have smallest orthogonal deviations, but largest parallel deviations. Strong synergy.
n_trials = 30 - np.random.randint(12, size=3)
block1 = pd.DataFrame({'user': 2, 'block': 1,
'parallel': np.random.normal(0.0, 7.0, size=n_trials[0]),
'orthogonal': np.random.normal(0.0, 3.0, size=n_trials[0])})
block3 = pd.DataFrame({'user': 2, 'block': 3,
'parallel': np.random.normal(0.0, 7.0, size=n_trials[1]),
'orthogonal': np.random.normal(0.0, 3.0, size=n_trials[1])})
# Block 2 has smaller parallel deviations than blocks 1&3, but higher orthogonal deviations than blocks 1&3.
# In this case sigma_block2_parallel - sigma_block2_orthogonal = 0
# No synergy.
block2 = pd.DataFrame({'user': 2, 'block': 2,
'parallel': np.random.normal(0.0, 5.0, size=n_trials[2]),
'orthogonal': np.random.normal(0.0, 5.0, size=n_trials[2])})
df2 = pd.concat((block1, block2, block3), axis='index')
df2[['parallel', 'orthogonal']] = df2.groupby('block')[['parallel', 'orthogonal']].transform(lambda x: x - x.mean())

# %%
df = pd.concat((df0, df1, df2), axis='index')

# %%
df[['user', 'block']] = df[['user', 'block']].astype('category')
df.sort_values(by=['user', 'block'], inplace=True)

# %%
# Reshape data into long format.
samples = df.melt(id_vars=['user', 'block'], value_vars=['parallel', 'orthogonal'],
var_name='direction', value_name='projection')
# Dimensions.
n_users = samples['user'].nunique()
n_blocks = samples['block'].nunique()
n_directions = samples['direction'].nunique()

# %%
# Vizualize the data.
g = sns.FacetGrid(samples, col='block', row='user', hue="direction")
g.despine(offset=10)
g = (g.map(sns.distplot, "projection").add_legend(title="Projection", loc='upper right'))
for ax in g.axes.flat:
ax.set_ylabel("Probability Density")
ax.set_xlabel("Projection Length")
plt.tight_layout()
plt.show()

# %%
#################################
#          Null Model           #
#################################
with pm.Model() as null_model:
null_model.name = "Null Model"
sigma = pm.Gamma("sigma", alpha=2.2, beta=0.4, shape=n_users)

obs = pm.Normal(f"obs", mu=0.0, sigma=sigma[samples['user']], observed=samples['projection'])

# %%
pm.model_to_graphviz(null_model)

# %%
with null_model:
idata0 = pm.sample(1000, tune=800, return_inferencedata=True)
az.summary(idata0, var_names=["Null Model_sigma"])
az.loo(idata0)

# %%
#################################
# Main Effect Projection Model  #
#################################
with pm.Model() as projection_model:
projection_model.name = "Main Effect Projection"
sigma_ortho = pm.Gamma("sigma_ortho", mu=3.34, sigma=1.44, shape=n_users)
# The sum of two gamma-distributed variables is gamma-distributed.
sigma_diff = pm.Gamma('sigma_diff', mu=4.0, sigma=5.0, shape=n_users)
sigma_parallel = pm.Deterministic("sigma_parallel", sigma_ortho + sigma_diff)

obs_parallel = pm.Normal("obs_parallel", mu=0.0,
sigma=sigma_parallel[samples.loc[samples['direction'] == 'parallel', 'user']],
observed=samples.loc[samples['direction'] == 'parallel', 'projection'])
obs_orthogonal = pm.Normal("obs_orthogonal", mu=0.0,
sigma=sigma_ortho[samples.loc[samples['direction'] == 'orthogonal', 'user']],
observed=samples.loc[samples['direction'] == 'orthogonal', 'projection'])

# %%
pm.model_to_graphviz(projection_model)

# %%
with projection_model:
idata1 = pm.sample(1000, tune=1000, return_inferencedata=True)
az.summary(idata1, var_names=["Main Effect Projection_sigma"], filter_vars='like')

# %%
#idata1.log_likelihood['projection_diff'] = (idata1.log_likelihood['Main Effect Projection_obs_parallel']
#                                            - idata1.log_likelihood['Main Effect Projection_obs_orthogonal'])
#az.loo(idata1, var_name='projection_diff')

# %%
##########################################
# Main Effects Block2 + Projection Model #
##########################################
with pm.Model() as both_main_effects_model:
both_main_effects_model.name = "2 Main Effects"
sigma_ortho13 = pm.Gamma("sigma_ortho_blocks_1_3", mu=3.34, sigma=1.44, shape=n_users)
# The sum of two gamma-distributed variables is gamma-distributed.
sigma_diff_proj = pm.Gamma('sigma_diff_projection', alpha=4., beta=1.0, shape=n_users)
sigma_diff_block2 = pm.Gamma('sigma_diff_block2', alpha=1., beta=1.0, shape=n_users)

sigma_ortho2 = pm.Deterministic('sigma_ortho_block_2', sigma_ortho13 + sigma_diff_block2)
sigma_parallel13 = pm.Deterministic('sigma_parallel_blocks_1_3', sigma_ortho13 + sigma_diff_proj)
sigma_parallel2 = pm.Deterministic('sigma_parallel_block_2', sigma_parallel13 + sigma_diff_block2)

data_parallel13 = samples.query('direction=="parallel" & block in [1, 3]')
data_ortho13 = samples.query('direction=="orthogonal" & block in [1, 3]')
data_parallel2 = samples.query('direction=="parallel" & block==2')
data_ortho2 = samples.query('direction=="orthogonal" & block==2')

obs_parallel13 = pm.Normal("obs_blocks_1_3_parallel", mu=0.0, sigma=sigma_parallel13[data_parallel13['user']],
observed=data_parallel13['projection'])
obs_parallel2 = pm.Normal("obs_block2_parallel", mu=0.0, sigma=sigma_parallel2[data_parallel2['user']],
observed=data_parallel2['projection'])
obs_orthogonal13 = pm.Normal("obs_blocks_1_3_orthogonal", mu=0.0, sigma=sigma_ortho13[data_ortho13['user']],
observed=data_ortho13['projection'])
obs_orthogonal2 = pm.Normal("obs_block2_orthogonal", mu=0.0, sigma=sigma_ortho2[data_ortho2['user']],
observed=data_ortho2['projection'])

# %%
pm.model_to_graphviz(both_main_effects_model)

# %%
with both_main_effects_model:
idata2 = pm.sample(1000, tune=1000, return_inferencedata=True)
az.summary(idata2, var_names=["2 Main Effects_sigma_diff"], filter_vars='like')

What’s missing in these models is the partial pooling, but I guess this could be easily done by just adding more priors on top. To get from these models to answering the research question is where the bridge sampling should come into play. I have to figure out how to get the traces and model in the vectorized form suited for the marginal_llk-function.

The Mixed Model approach does have advantages like partial-pooling and being able to include within- and between-conditions. The disadvantage is that it is not as intuitive for me for the research question at hand and fitting one model without restrictions also takes me only halfway to what I really want to know.

Let me illustrate the same predictions from above as regression models. The i stands for intercept and the b are the slopes. i12 is abbreviation for intercept block 1 & intercept block 2. The orthogonal direction is coded as 0 and parallel as 1.

    Main Effect Projection  Main Effect Block2    Main Effects Projection + Block2
sigma                   sigma                 sigma  / b2
^                       ^                     ^     /
|      / b123           |                     |    /
|     /               i2|_________b2          |   /  / b13
|    /                  |                     |  /  /
|   /                   |                     | /  /
|  /                    |                   i2|/  /
| /                     |                     |  /
i123|/                   i13|_________b13         | /
|                       |                  i13|/
|---------->            |---------->          |---------->
0      1                0      1              0      1
ortho  parallel         ortho  parallel       ortho  parallel

Block x Projection v1   Block x Proj. v2      Block x Projection v3
sigma                   sigma                 sigma
^                       ^                     ^
|      / b13            |      / b13          |      / b1
|     /                 |     /               |     /
|    /                  |    /                |    /
i2|___/___ b2             |   /                 |   /
|  /                    |  /                  |  /
| /                     | /                   | /
i13|/                  i123|/______ b2       i123|/______ b23
|                       |                     |
|---------->            |---------->          |---------->
0      1                0      1              0      1
ortho  parallel         ortho  parallel       ortho  parallel

Block x Projection v4   Null Model
sigma                   sigma
^                       ^
i2|______                 |
| b2  /                 |
|    /                  |
|   / b13               |
|  /                i123|________
| /                     |  b123
i13|/                      |
|                       |
|---------->            |---------->
0      1                0      1
ortho  parallel         ortho  parallel

If I fitted a model without restrictions on intercepts and slopes I would have to compare the estimates I get for all of them to these hypothesis which describe the illustrations above in order to answer my research question:

Main Effect Projection:
intercepts and slopes are fixed
b1 = b2 = b3 > 0
i1 = i2 = i3

Main Effect Block2:
b1 = b2 = b3 (fixed slopes)
i2 > i1, i3
i3 - i1 = 0

Main Effects Projection + Block2:
i2 > i1, i3
b1 = b2 = b3

Block x Projection v1:
i3 - i1 = 0
i2 > i1, i3
b3 - b1 = 0
b3, b1 > 0
b3 > i2 - i3
b1 > i2 - i1
b2 = 0

Block x Proj. v2:
i1 = i2 = i3
b1, b3 > 0
b3 - b1 = 0
b2 = 0

Block x Projection v3:
i1 = i2 = i3
b1 > 0
b2 = b3 = 0

Block x Projection v4:
i1 = i3
i2 > i1, i3
b2 = 0
b1 = i2 - i1
b3 = i2 - i3

Null Model:
i1 = i2 = i3
b1 = b2 = b3 = 0

One way I could imagine to do it quantitavely is to take the estimates as input to models with these restrictions and get the likelihoods. But the “intuitve” solution with bridge sampling appears more straight forward to me.
What do you think?

Thanks for the explanation, personally, I am not sure the data could answer the question like:“hey after comparing 8 different model here, we are sure that model X is the generative process”. And building a 8*8 comparison table of all the models seems a bit difficult to reason. The mixed effect model is more like: “we assume every conditions are different, and build a full model the encodes all the information, then look at the posterior to see if they are really all different (i.e., no overlapping)” - basically you plot the posterior and see which ASCII-diagrams for different hypotheses it fits the best.
But since you have strong theoretical motivation, it certainly sound like a reasonable approach.

I agree with your opinion about a less restrictive hierarchical model having the potential to reveal more about the underlying data. I’ve tried to work my way up to the full model after reading this informative Primer on Bayesian Methods for Multilevel Modeling, but I’m stuck early on after getting to the user-level and trying to add the block-level.
Using the dims keyword argument and coords seems more elegant than the indexing scheme from this post.

Example data with multilevel predictors

    user	experience condition  block rating	direction	projection
0	9	    0.0        A      	  1	    2.0     parallel	-1.058687
1	9	    0.0        A      	  1	    2.0     parallel	0.201144
2	9	    0.0        A      	  1	    2.0     parallel	-0.806420
...

experience is a user-level predictor and rating a block-level “predictor” (actually, rating is measured afterwards, so it’s correlated with block).

The indexing/design matrix still poses a problem to me. I don’t understand how to use the predictors when doing outcome, predictors = dmatrices("projection ~ condition * user * block * direction", df). The example didn’t utilize the design matrix, stating the glm() function does not play nice with hierarchical models yet.
This is an excerpt from what I get with aforementioned formula.

Terms:
‘Intercept’ (column 0)
‘condition’ (columns 1:3)
‘user’ (columns 3:32)
‘block’ (columns 32:34)
‘direction’ (column 34)
‘condition:user’ (columns 35:93)
‘condition:block’ (columns 93:97)
‘user:block’ (columns 97:155)
‘condition:direction’ (columns 155:157)
‘user:direction’ (columns 157:186)
‘block:direction’ (columns 186:188)
‘condition:user:block’ (columns 188:304)
‘condition:user:direction’ (columns 304:362)
‘condition:block:direction’ (columns 362:366)
‘user:block:direction’ (columns 366:424)
‘condition:user:block:direction’ (columns 424:540)

From what I gathered from a tutorial in R about rat growth that @OriolAbril provided, in R this formula would be kind of this \sigma_{proj} \sim (condition * block * direction)/(1|trial) + (1|user), meaning trials are viewed as nested random effects, and users as random effects. It seems patsy doesn’t support nesting and random effects.

Following the primer I’ve come so far as to include users, as they are what counties are in the radon example. The difference is that I model the standard deviation sigma, not the mu. When trying to add rating as a pendent to uranium in the radon example, I ran into broadcasting issues.

Model with only direction x user:

conditions = df['condition'].cat.categories.tolist()
n_conditions = df['condition'].nunique()
condition_indices = df['condition'].cat.codes.values # Condition index for every single value.
condition_lookup = pd.Series(*pd.factorize(df['condition'].cat.categories), name='condition_idx_lookup')

# We need to be able to index by user. But the user IDs of the sample don't start at 0. We create an index and a mapping between them.
users = df['user'].cat.categories.tolist()
n_users = df['user'].nunique()
user_indices = df['user'].cat.codes.values  # User index for every single value.
user_lookup = pd.Series(*pd.factorize(df['user'].cat.categories), name='user_idx_lookup')
# User level predictors. More options: age group, gender.
# experience is not on an interval scale (no equidistance between items), but ordinal.
experience = df.drop_duplicates(subset="user")['experience']  # Shape = user_indices

blocks = df['block'].cat.categories.tolist()
n_blocks = df['block'].nunique()  # Same as using cat.categories.size
block_indices = df['block'].cat.codes.values # Block index for every single value.
block_lookup = pd.Series(*pd.factorize(df['block'].cat.categories), name='block_idx_lookup')
# Block level predictors.
# Note: Ratings are correlated with block index.
ratings = df.drop_duplicates(subset=["user", 'block'])['rating']  # Shape = block_indices

# Coded direction of projection: 0 = orthogonal, 1 = parallel.
directions = df['direction'].cat.categories.tolist()
n_directions = df['direction'].nunique()
direction_indices = df['direction'].cat.codes.values
direction_lookup = pd.Series(*pd.factorize(df['direction'].cat.categories), name='direction_idx_lookup')

coords = {'Direction': directions, 'User': users, 'Block': blocks, 'Condition': conditions, 'obs_id': np.arange(direction_indices
.size)}

with pm.Model(coords=coords) as varying_intercept_slope_by_user:
direction_idx = pm.Data('direction_idx', direction_indices, dims='obs_id')
user_idx = pm.Data('user_idx', user_indices, dims='obs_id')

# Hyperpriors:
intercept = pm.Gamma('intercept', mu=4.0, sigma=2.0)
intercept_sd = pm.Exponential('sigma_intercept', 1.0)
slope = pm.Normal('slope', mu=0.0, sigma=1.0)
slope_sd = pm.Exponential('sigma_slope', 0.5)

# Varying intercepts:
intercept_user = pm.Gamma('intercept_user', mu=intercept, sigma=intercept_sd, dims='User')
# Varying slopes:
slope_user = pm.Normal('slope_user', mu=slope, sigma=slope_sd, dims='User')

# Expected deviation per user:
theta = intercept_user[user_idx] + slope_user[user_idx] * direction_idx

projection = pm.Normal('projection', 0.0, sigma=theta, observed=df['projection'], dims='obs_id')

How do I add the other levels? Condition would be like state in the radon example. I try, but I keep hitting walls. In the meantime I try again and share my progress. I guess I have to formulate linear expressions for each level and use those as inputs to the levels below.

I just realized that using the ratings (for blocks) as I did there isn’t going to work with how I defined blocks. Block indices are only 0,1,2, but each user gave an individual rating for each block leading to n_users*n_blocks ratings. If it was a block-level predictor, there would only be 3 values, right? I could make an index for each block (0-89), but then I wouldn’t be able to compare each first, second and third block across users.

I tried adding a new level, but I’m stuck with the broadcasting issue. I really don’t know what I’m supposed to do with the dims or shape keyword arguments and the slicing of the parameters. Unfortunately, I didn’t find a comprehensible documentation on how this works anywhere. All the examples just have models with 1 or 2 levels, nothing nested deeper than that.

What I’d like to get is an estimate for intercept and slope of projection variance for Blocks 1,2, and 3, while each user can have their own intercepts and slopes in each block as different baselines (random effect).

In addition to the above I tried with blocks (0,1,2) and an ID for each single block.

coords['Block_ID'] = block_ids  # 0..n_users*n_blocks

with pm.Model(coords=coords) as by_block_model:
user_idx = pm.Data('user_idx', user_indices, dims='obs_id')
block_id_idx = pm.Data('block_id_idx', block_id_indices, dims='obs_id')
block_idx = pm.Data('block_idx', block_indices, dims='obs_id')
direction_idx = pm.Data('direction_idx', direction_indices, dims='obs_id')
experience = pm.Data('experience', gaming_exp, dims='User')
rating = pm.Data('rating', ratings, dims='Block_ID')  # Block has indices 0,1,2. But ratings are size of n_users*n_blocks. Try block_id?

# Hyperpriors:
global_intercept = pm.Gamma('global_intercept', mu=4.0, sigma=2.0)
global_intercept_sd = pm.Exponential('global_intercept_sd', 1.0)
global_slope = pm.Normal('global_slope', mu=1.5, sigma=1.0)  # We expect the slope to be positive most of the time.
global_slope_sd = pm.Exponential('global_slope_sd', 0.5)

# Users.
# Varying intercepts per user:
intercept_user = pm.Gamma('intercept_user', mu=global_intercept, sigma=global_intercept_sd, dims='User')
intercept_user_sd = pm.Exponential('intercept_user_sd', 1.0, dims='User')
# Varying slopes per user:
slope_user = pm.Normal('slope_user', mu=global_slope, sigma=global_slope_sd, dims='User')
slope_user_sd = pm.Exponential('slope_user_sd', 1.0, dims='User')

# Blocks.
# Varying intercepts per block:
intercept_block = pm.Gamma('intercept_block', mu=intercept_user, sigma=intercept_user_sd, dims='Block')  # FixMe: ValueError: operands could not be broadcast together with shapes (3,) (30,)
# Varying slopes per block:
slope_block = pm.Normal('slope_block', mu=slope_user, sigma=slope_user_sd, dims='Block')

# Expected deviation per block:
theta = intercept_block[block_idx] + slope_block[block_idx] * direction_idx

projection = pm.Normal('projection', 0.0, sigma=theta, observed=df['projection'], dims='obs_id')

Once I know how to add more levels, I’d like to cluster users by condition, which may have an effect on standard deviation.

I tried resolving the broadcasting issues by repeating the intercept_user and slope_user parameters.
block_id_idx 0, 1, 2 belong to user 0, block_id_idx 3, 4, 5 to user 1 and so on. intercept_user should then be like 0, 0, 0, 1, 1, 1, 2, 2, 2 etc.

import theano.tensor as tt

with pm.Model(coords=coords) as by_block_model:
user_idx = pm.Data('user_idx', user_indices, dims='obs_id')
block_id_idx = pm.Data('block_id_idx', block_id_indices, dims='obs_id')
block_idx = pm.Data('block_idx', block_indices, dims='obs_id')
direction_idx = pm.Data('direction_idx', direction_indices, dims='obs_id')
experience = pm.Data('experience', gaming_exp, dims='User')
rating = pm.Data('rating', ratings, dims='Block_ID')  # Block has indices 0,1,2. But ratings are size of n_users*n_blocks. Try block_id?

# Hyperpriors:
global_intercept = pm.Gamma('global_intercept', mu=4.0, sigma=2.0)
global_intercept_sd = pm.Exponential('global_intercept_sd', 1.0)
global_slope = pm.Normal('global_slope', mu=1.5, sigma=1.0)  # We expect the slope to be positive most of the time.
global_slope_sd = pm.Exponential('global_slope_sd', 0.5)

# Users.
# Varying intercepts per user:
intercept_user = tt.repeat(pm.Gamma('intercept_user', mu=global_intercept, sigma=global_intercept_sd, dims='User'), n_blocks)
intercept_user_sd = tt.repeat(pm.Exponential('intercept_user_sd', 1.0, dims='User'), n_blocks)
# Varying slopes per user:
slope_user = tt.repeat(pm.Normal('slope_user', mu=global_slope, sigma=global_slope_sd, dims='User'), n_blocks)
slope_user_sd = tt.repeat(pm.Exponential('slope_user_sd', 1.0, dims='User'), n_blocks)

# Blocks.
# Varying intercepts per block:
intercept_block = pm.Gamma('intercept_block', mu=intercept_user, sigma=intercept_user_sd, dims='Block_ID')
# Varying slopes per block:
slope_block = pm.Normal('slope_block', mu=slope_user, sigma=slope_user_sd, dims='Block_ID')

# Expected deviation per block:
theta = intercept_block[block_idx] + slope_block[block_idx] * direction_idx

projection = pm.Normal('projection', 0.0, sigma=theta, observed=df['projection'], dims='obs_id')

As you can see the model gets created now. There are a few Data containers floating around. I’m using the Block_ID (0…89) to create a slope and intercept for each single block in the data (30 users * block1, 30 * block2, 30 * block3). But since I want to estimate intercept and slope across users for each block (0,1,2), I used block_idx (0,1, or 2 for each data point) in the linear regression for the expected standard deviation in the data theta = intercept_block[block_idx] + slope_block[block_idx] * direction_idx. When using block_idx (either 0, 1, or 2) instead of block_id_idx (0…89), sampling time increases like 10x. Regardless of whether I choose block_idx or block_id_idx, I get every sampling error you can imagine: thousands of divergences, rhat > 1.05, acceptance probability, number of effective samples is smaller than 200.
If I used block_id_idx in the linear expression, is there a way to sample from the posterior using the block_idx? Like group every first, second and third block and get their group estimates?

How can I reparameterize the model? The between factor of condition isn’t even in there yet.
I just want estimates on the differences between standard deviations for blocks 1, 2 and 3, which each user does (within), but users are in different conditions (between). As a bonus it would be nice to know if experience (user-level) has an influence, or the individual rating of blocks. Doesn’t sound too hard.

I often read a common answer to some topics that everything should be straight forward to do in PyMC3. Unfortunately I don’t have a PhD in Bayesian Statistics and 15 years of computer science under the belt, so I appreciate any help! To me it’s everything but straight forward. I’m on a tight schedule, desperate and on the verge of going back to frequentist statistics. The devil you know…

IIUC, you have condition, block, direction which potentially be different, and user could also display different strategy - in that case, I recommend you to try parsing the linear formula σproj∼(condition∗block∗direction*user) in patsy, and use indexing for trial (the random effect.)

In general, frequentist testing like ANOVA are GLMs, and GLMs are pretty straightforward to express in Bayesian - I would love to give you a bit more concrete suggestion, maybe it is easier you share a jupyter notebook with some simulation data?

That sounds about right. I don’t really care about trials, e.g. if all trials with number 1 differ or not. I did try to use patsy and got an outcome and a DesignMatrix. Unfortunately, I don’t have the slightest clue about what to do with the DesignMatrix object. I tried multiplying it with a distribution and only got an AssertionError. I also don’t know what to index where?

Yes, and you only have to drag some variables into the right slots in a GUI software most of the time.

That is very kind of you and much appreciated! I prepared a gist with some mocked up data. There’s not much of a model, because I don’t know how to use the DesignMatrix and indexing.


text:
2014-08-08 - ImageJ 2.0.0-rc-11 released

Today, the ImageJ team is pleased to announce a new public release candidate for ImageJ2: version 2.0.0-rc-11.

New features

• Usage tracking! We can now track and upload (anonymously) use counts at the plugin level, along with information about the environment of use (country, java version, language, operating system, time zone). Although we are taking care not to expose any user data with these statistics, this feature can be disabled via a new Edit › Options › Privacy... setting. Expect great things coming to the usage page! (note: “great things” = statistics graphs)
• All images, whether opened with the SCIFIO option enabled or not, should now be automatically populated in @Dataset annotations in scripts and plugins. This will make parameterized scripting much more consistent and useful. See the scripting guide for more information and tutorials.
• More plugin types! Continuing our dedication to making ImageJ2 a flexible and extensible platform, we converted some lower level utilities to pluggable services: type conversion and preference storage. This won’t mean much to users right now, but it gives developers new ways to tap in to the functionality of ImageJ, and plays a key role in extensions like our upcoming MATLAB scripting language!

Bug-fixes

• We resolved a significant memory leak due to images not being closed properly. This was also causing naming errors (e.g., the second time you opened “Data” it would display a title “Data-1” even after closing the first image) and other terrible gotchas. As a consequence, the linkage between ImageJ2 and ImageJ 1.x data structures maintained by the legacy layer is now more robust.
• We fixed several bugs in SCIFIO; e.g., ImageJ slice labels should now be preserved after saving and reopening in TIFF format.
• Quite a few other small bug-fixes; for example, the languages shown in the Languages menu of the Script Editor now appear in a consistent and alphabetical order.

How to test the release candidate?

Use the HelpUpdate command to update your ImageJ installation.

Please send any comments to the ImageJ mailing list. And thank you very much for trying ImageJ2!

format: plain


text:
What does the big cap notation mean?

I'm trying to understand How can an ordered pair be expressed as a set? and I don't know what the big cap/cup notations mean when placed next to an ordered pair: $\bigcap(a,b)$ and $\bigcup(a,b)$.

-
It mean union $\cup$ or intersection $\cap$ of sets in $(a,b)$. – Sasha Oct 5 '11 at 18:29
@Sasha So it only makes sense if $a$ and $b$ are sets. :) – Paul Manta Oct 5 '11 at 18:30
The cap is intersection, the cup is union. Remember that the ordered pair $(a,b)$ is just shorthand notation for the set $\{ \{ a \} , \{a,b\} \}$. – Ragib Zaman Oct 5 '11 at 18:32
In the question you linked to $(a,b)$ represented a set of sets. Then $\cup (a,b)$ denoted the union of those sets. – Sasha Oct 5 '11 at 18:33

This is actually answered in the linked question, but for clarification, if by definition $(a,b)=\{\{a\},\{a,b\}\}$ then $$\bigcap(a,b) = \bigcap\{\{a\},\{a,b\}\} = \{a\} \cap \{a,b\} = \{a\}$$ and $$\bigcup(a,b) = \bigcup\{\{a\},\{a,b\}\} = \{a\} \cup \{a,b\} = \{a,b\}.$$


text:
# What's wrong with this proposed resolution to the St Petersburg Paradox?

We have a game where your payout is $2^k$ where $k$ is the number of times you flipped a coin to land on heads (if your first flip is a head, then $k=1$). Then the expected payout is: $$E = \frac{1}{2}(2) + \frac{1}{4}(4) + \frac{1}{8}(8)+...$$ $$E=1+1+1+...$$ $$E=\infty$$

How much should I pay to play this game?

Well we know from the geometric distribution the expected number of coins I'll flip until getting a head is:

$$\frac{1}{P(HEAD)} = \frac{1}{.5}=2$$

So I will pay anything less than $2^k$ with $k=2$:

i.e. < 4 dollars

• "well" needs references. Where does "expected" come from - give a single example before describing a sum. Bold the actual question. Find a way to augment the "expected times to get a head" ... I recommend simulation. Does the flip that made the head count? If you flip it once and get a head, is that zero before, or one before? Can you prove that you accounted for the "max" and not just the "mean" - where is the right tail? Kelly Criterion? (en.wikipedia.org/wiki/Kelly_criterion) – EngrStudent Dec 5 '16 at 1:38
• If you value the gamble at \$4, how about I pay you \$5 to play a shortened version of the St. Petersburg gamble where I can go at most 10 rounds? Deal? – Matthew Gunn Dec 5 '16 at 3:02
• @EngrStudent The "expected" number of flips is equal to the mean of a Geometric Random Variable. I think you're right in that I should be more explicit about which range I'm using {0,1,2,...} or {1,2,3...}. – timwiz Dec 5 '16 at 12:38
• @MatthewGunn touché – timwiz Dec 5 '16 at 12:39

• Let $K$ be some random variable.

• In your problem, $K$ is number of times you flip before getting heads.
• Let $f(k)$ be some payoff function.

• In your problem $f(k) = 2^k$.
• Let $f(K)$ be the payoff

You're saying that a reasonable valuation of the gamble $f(K)$ is given by $f(\mathrm{E[}K])$. This is an entirely ad-hoc, rather unprincipled heuristic. Perhaps fine in some situations (eg. where $K$ is small and $f$ near linear), but it's easy to construct an example where it suggests something non-sensical.

## Example where your system makes absolutely no sense

Let $K$ be a draw from the normal distribution $\mathcal{N}(0,10000000000000)$ and let the payoff function be $f(K) = K^2$. Your system says I shouldn't pay more than $0$ for this gamble because $f(\mathrm{E}[K]) = 0^2 = 0$. But shouldn't you assign some positive value to this gamble?! There is a 100% probability the payoff is greater than zero!

## A more classic resolution of the St. Petersburg Paradox

One approach is to add risk-aversion. If you're sufficiently risk averse, what you're willing to pay to play this infinite expectation gamble will be finite. If you accept the Von Neumann-Morgernstern axioms, then the certainty equivalent of playing the game is given by $z$ where:

$$u(w + z) = \mathrm{E}[ u(w + f(K)) ]$$

and where $w$ is your wealth and $u$ is a concave function (in jargon, a Bernoulli utility function) which captures your level of risk aversion. If $u$ is sufficiently concave, the valuation of $2^K$ will be finite.

A Bernoulli utility function with some nice properties turns out to be $u(x) = \log(x)$. Maximizing expected utility where the Bernoulli utility function is the log of your wealth is equivalent to maximizing the expected growth rate of your wealth. For simple binary bets, this gives you Kelly Criterion betting.

An important other point is that the risk aversion approach leads to different certainty equivalents depending on what side of the gamble you are on.

• (+1) Furthermore, specifying the distribution of a payoff does not uniquely determine $K$ and $f$, so the heuristic described in the first paragraph would need some way to deduce the canonical $f(K)$ representation. – Juho Kokkala Dec 6 '16 at 12:26

There's nothing wrong with that proposed resolution.

In the original paradox we look at the expected value (mean) of the profit which is infinite and therefore you should stake an infinite amount. However, after the first flip of the coin there's a 50% chance that you've lost money and that is why people don't like it. Your resolution just formalizes this, instead of looking at the mean profit you are looking at the median profit. Unlike the mean profit, the median profit is finite and the paradox goes away.

• Why median profit would be relevant? Imagine following gamble: 1% of getting nothing, 1% of getting a million $and 98% of geting$1. How much would you pay to play that? – Daerdemandt Dec 5 '16 at 7:58
• @Daerdemandt overreliance on mean is what causes this "paradox", If you bet a huge amount on a game like this then there's a 50% chance of being bankrupt and a chance of winning millions. It's not irrational to try to avoid bankruptcy. – Hugh Dec 5 '16 at 10:24
• Unlimited money > huge money. Hence, money-wise, it would be pretty rational to play if you're fine with the game taking arbitrarily long time. These are assumptions we make to frame statistical problem as more real-sounding. Realistically, you'd have trouble with waiting several lifetimes to get your unlimited money, your marginal utility for money can depend on amount of money you have (heck, you can even lack utility function) - and you'd be pretty skeptical with strangers offering you to gamble to win unlimited money. But all of this is outside of statistics' scope. – Daerdemandt Dec 5 '16 at 11:50
• In general no, we don't have to. In statistical problem of similar kind we are expected to. (here's the game, how much would you pay to play?). Money is just a fluff to make the problem less boring. If you want to challenge the fluff, you're free to do so but that's out of statistics scope. – Daerdemandt Dec 5 '16 at 13:23
• I'm not sure exactly if people in general attach specific enough meaning to the word "paradox". This is a hypothetical situation that stretches some of used abstractions thin. It is useful to demonstrate that these abstractions are not perfect. With that reading, I personally would lump this situation with other so-called "paradoxes", your usage may be different. – Daerdemandt Dec 5 '16 at 14:28

If I understand correctly, your analysis is:

1. Calculate the expected number of coin-flips required to get a head.
2. Calculate the payout for the outcome where you get exactly the expected number.
3. Value the game equal to that payout.

...OK, let's modify that game a little bit. Just like the original version, I will flip a coin and keep flipping until I throw heads. Only the payouts have changed:

• If I flip heads on the second throw, you get four dollars.
• On any other outcome, you lose everything you own and have to come work for me forever, for free.

How many coins do we expect to flip before we get a head? 2, exactly the same as before.

What is the payout for the outcome where we flip two coins to get a head? $4.00, exactly the same as before. How much would you be willing to pay for the 'privilege' of paying this game that has a 75% chance of bankrupting you and a 25% chance of returning$4.00?

I suspect the answer is not "up to four dollars, exactly the same as before". Which means there's a hole in your logic.

Taking a broader perspective, expected winnings are not necessarily enough information to answer this sort of question; usually it depends on some additional context. Is this a one-off opportunity or are you expecting to be offered this gamble many times? How much money do you have on hand? And how much money do you need to be happy?

For example, if my total wealth is $100 but I urgently need a million dollars for a life-saving operation, I would be willing to pay all my money for a single shot at the St. Petersburg gamble. It only gives me a 1/2^19 chance of winning the money I need, but if I don't play I have no chance at all. On the other hand, if my total wealth is$1000,000 and I need exactly a million dollars for that operation, the most I'd be willing to pay for a single game is two dollars (which I'm guaranteed to win back). Anything more, and I have a 1/2 chance of ending up short of the million bucks I need to save my life.

If I'm expecting to have many chances to play such games, then I probably want to choose a strategy that gives me a high probability of having lots of money at the end of all those games. For example:

Game A is guaranteed to increase my wealth by 10% every time I play it. (Expected winning: +10% of my current wealth.) Game B has a 90% chance of doubling my wealth, and a 10% chance of bankrupting me. (Expected winning: +70% of my current wealth.) [edit: actually +80% because I fail at basic arithmetic, but the argument still holds.]

If I play 100 iterations of Game A, I'm certain to multiply my wealth by 13,780 times.

If I play 100 iterations of Game B, I have a 0.0027% chance of becoming unimaginably wealthy (about 10^30 x what I started with)... and a 99.73% chance of being bankrupted. Even though the average is better than for Game A, it's not a good option.

For this sort of heavily-iterated game, rather than trying to maximise my expected winnings in each game, I'm better off trying to maximise expected value of ln(total wealth after game/total wealth before game). This ensures long-term growth without getting wiped out.

If the stakes for every game are small relative to my total wealth, then this is approximately equivalent to maximising expected winnings in each game.

So, if you're playing lots of games and never risking a large portion of your current wealth, then the expected value of the gamble tells you all you need to know. In just about any other situation, you need to think about other things too.

• +1 Thank you for a thoughtful answer that really gets to some of the fundamental underlying issues. Welcome to our site! – whuber Dec 7 '16 at 16:29


text:
By Topic

# Electronics Letters

## Filter Results

Displaying Results 1 - 25 of 37
• ### Optical DPSK heterodyne detection experiments using DBR laser diodes with external optical feedback

Publication Year: 1985, Page(s):1121 - 1122
Cited by:  Papers (6)
| | PDF (279 KB)

High receiver sensitivity, with 7 dB improvement over a direct detection system, has been achieved in a 400 Mbit/s optical DPSK heterodyne detection experiment using DBR laser diodes with external optical feedback. The influence of the laser output phase noise was evaluated experimentally in good agreement with theory. View full abstract»

• ### Automatic timing alignment for regenerative repeaters

Publication Year: 1985, Page(s):1122 - 1123
Cited by:  Papers (3)
| | PDF (284 KB)

A novel circuit technique that provides automatic timing alighnment of digital regenerative repeaters is described. The scheme employs correlation feedback to maintain the correct clock and data phase relationship. Theoretical and practical results are presented. View full abstract»

• ### High-rate fabrication of wholly synthesised fibre preforms by the multiflame VAD method using SiHCI3 raw materials

Publication Year: 1985, Page(s):1123 - 1124
Cited by:  Papers (3)
| | PDF (315 KB)

High-rate fabrication techniques for wholly synthesised graded-index fibre preforms were developed. An overall deposition rate of 5.6 g/min was attained using a multiflame VAD burner and using SiHCl3. A low loss of 0.25 dB/km at 1.59 Â¿m and a small eccentricity of less than 0.4% were obtained. View full abstract»

• ### GaAlAs/GaAs heterojunction bipolar phototransistors grown by LPE with a current gain of 50 000

Publication Year: 1985, Page(s):1124 - 1126
Cited by:  Papers (2)
| | PDF (493 KB)

A simple mesa GaAlAs/GaAs heterojunction bipolar phototransistor fabricated by LPE growth technique is presented. A high current gain of the order of 50 000 has been obtained and is believed to be the highest ever reported for a bipolar device. The corresponding high sensitivity is about 15 000 A/W for an 800 nm illumination wavelength. View full abstract»

• ### Floating integrator technique for switched-capacitor circuits employing reset cycles

Publication Year: 1985, Page(s):1126 - 1128
| | PDF (412 KB)

A simple floatingÂ¿ integrator technique is presented that is applicable to both single-ended and differential switched-capacitor circuits which employ reset cycles. This method is useful in eliminating the effects of any input offset or interstage differential or common-mode offsets in DC-coupled circuits. A monolithic circuit is demonstrated which uses this technique to achieve accurate differe... View full abstract»

• ### Choice of PRF for error-tolerant range destaggering of radar targets

Publication Year: 1985, Page(s):1128 - 1129
| | PDF (308 KB)

High and medium PRF radars perform ranging by means of multiple PRFs chosen according to the Chinese remainder theorem. Unfortunately this theory does not cope with errors in range measurements. A new rule for the choice of the PRFs is proposed which is tolerant of a prescribed error amplitude in range measurements. View full abstract»

• ### Linewidth of 1.3μm cleaved-coupled-cavity lasers

Publication Year: 1985, Page(s):1129 - 1131
| | PDF (429 KB)

The linewidth and side-mode suppression ratio of three variable-gap C3 lasers were measured as a function of output power and drive currents. With feedback carefully eliminated, the narrowest linewidth obtained was 8 MHz at an output power of 13.6 mW. There was little correlation between line-width and suppression ratio. View full abstract»

• ### Fabrication of long single-mode and multimode fluoride glass fibres by the double-crucible technique

Publication Year: 1985, Page(s):1131 - 1132
Cited by:  Papers (1)
| | PDF (352 KB)

Long single-mode and multimode fluoride glass fibres with core-clad structure are fabricated for the first time by the double-crucible technique. 500 m-long single-mode fibre and 700 m-long multimode fibre are obtained using ZrF4-BaF2-NaF-LaF3-AlF3 core and ZrF4-BaF2-NaF-LaF3-AlF3-HfF4 clad glasses. View full abstract»

• ### Meteor burst polarisation trials

Publication Year: 1985, Page(s):1132 - 1133
Cited by:  Papers (1)
| | PDF (355 KB)

Measurements over an 850 km path show that radio waves reflected by ionised meteor trails are typically received with the same polarisation as they were transmitted. Furthermore, a vertically polarised wave seems to suffer a 3Â¿4 dB larger loss than a horizontally polarised one. View full abstract»

• ### Mutual interaction between solitons of unequal amplitudes in optical fibre

Publication Year: 1985, Page(s):1133 - 1134
Cited by:  Papers (10)
| | PDF (334 KB)

An exact solution to the mutual interaction of two solitons propagating in a nonlinear single-mode fibre is given. From this, it is shown that the soliton separation can be effectively maintained by allowing the solitons to have unequal amplitudes. This provides a simple way to maintain the high bandwidth of a nonlinear fibre. View full abstract»

• ### The SINFET: a new high conductance, high switching speed MOS-gated transistor

Publication Year: 1985, Page(s):1134 - 1136
Cited by:  Papers (6)
| | PDF (522 KB)

A new MOS power semiconductor device with a very low on-resistance and a switching speed comparable to conventional n-channel power MOSFETs is described. The fabrication process is similar to that of an n-channel lateral DMOS transistor but with the conventional high-low ohmicÂ¿ drain contact replaced by a Schottky contact. In operation, the Schottky contact injects minority carriers to conductiv... View full abstract»

• ### Dark noise of Ga0.47 In0.53As photoconductive detectors

Publication Year: 1985, Page(s):1136 - 1138
| | PDF (492 KB)

We have studied the dark noise of planar, interdigitated Ga0.47In0.53As photoconductive detectors by measuring the statistical distribution of the dark current under DC bias. The measurements reveal two interesting results: (i) the probability distribution of the dark current around its DC level is Gaussian and (ii) the standard deviation of the probability distribution grows exponentially with in... View full abstract»

• ### Molecular-beam-epitaxy-grown single-longitudinal-mode GaAs-AlGaAs interferometric lasers on stepped-channelled substrates

Publication Year: 1985, Page(s):1138 - 1140
| | PDF (356 KB)

The GaAs-AlGaAs double-heterostructure interferometric lasers were successfully grown by molecular beam epitaxy on [011]-oriented-stepped-channelled GaAs(100) substrates. These lasers have a side-mode suppression ratio of greater than 200:1. View full abstract»

• ### Detection of fixed targets embedded in log-normal clutter

Publication Year: 1985, Page(s):1140 - 1141
| | PDF (288 KB)

Measurements of radar clutter using high-resolution systems show that clutter echoes cannot be modelled as Gaussian-distributed. A log-normal distribution is assumed and single-pulse detection performances are evaluated according to the Marcum-Swerling approach as well as the signal losses with respect to the Gaussian case. View full abstract»

• ### Measurements of crosspolarisation and gain reduction due to dust on a microwave reflector antenna

Publication Year: 1985, Page(s):1141 - 1143
Cited by:  Papers (1)
| | PDF (448 KB)

9 GHz measurements of effects of dust accretion on the surface of a reflector antenna confirm an earlier prediction of a larger gain reduction when the dust layer covers half the area. Crosspolarisation levels were higher than expected theoretically. Deposits on the primary feed also cause high crosspolarisation. View full abstract»

• ### Analysis of reduced chirping of semiconductor lasers for improved single-mode-fibre transmission capacity

Publication Year: 1985, Page(s):1143 - 1145
Cited by:  Papers (3)
| | PDF (425 KB)

It is shown that direct modulation of semiconductor lasers yielding optical pulses with an exponential shape corresponds to a nearly chirp-free operation. With such pulses the maximum transmission capacity at Â¿ = 1.55 Â¿m is about 55 Gbit/s Â¿km for a nondispersion-shifted single-mode fibre. View full abstract»

• ### Low-OH-content fibre fabrication using particle-size control sol-gel method

Publication Year: 1985, Page(s):1145 - 1146
| | PDF (272 KB)

By the sol-gel method using a new `particle-size controlÂ¿ technique, a low-OH content at a level of 0.1 parts in 106 and a minimum optical loss of 1.8 dB/km at 1.6 Â¿m wavelength has been achieved in pure silica core single-mode optical fibres. Improvement of dry gel properties using particle-size control is essential to bubble-free fibre drawing form dehydrated gel-derived glasses. View full abstract»

• ### Linewidth reduction of cleaved-coupled-cavity lasers by optical feedback from a single-mode polarisation-preserving fibre external cavity

Publication Year: 1985, Page(s):1146 - 1148
Cited by:  Papers (1)
| | PDF (430 KB)

Linewidth characteristics of a cleaved-coupled-cavity laser with a single-mode polarisation-preserving fibre external cavity have been investigated as a function of feedback power level and side-mode suppression ratio. The optimum power feedback ratio is near 2Ã10Â¿4 for a minimum linewidth of about 30 kHz. A weak dependence of linewidth on side-mode suppression ratio is observed as long as the s... View full abstract»

• ### Automatic threshold-level control (ATC) for optimum decision operation in optical repeaters

Publication Year: 1985, Page(s):1148 - 1149
| | PDF (308 KB)

An automatic threshold-level control (ATC) method for optimum decision operation in an optical repeater is studied using in-service error monitoring in the line code with the bit-inversion rule. Results of an experiment with a 110 Mbit/s optical repeater and a 10B1C line coder clarify that the ATC method is effective for maintaining the optimum decision threshold level. View full abstract»

• ### Measurement of wavelength variation of mode field radius using far-field pattern method

Publication Year: 1985, Page(s):1149 - 1151
Cited by:  Papers (2)
| | PDF (426 KB)

The far-field pattern measuring technique is used to measure the wavelength dependence of the mode field radius in standard single-mode fibres. Chromatic dispersion calculated from these results is shown to be accurate using as few as three lasers. Both theoretical and experimental results are presented. View full abstract»

• ### Dependence of single-longitudinal-mode probability on DFB laser facet structure

Publication Year: 1985, Page(s):1151 - 1152
Cited by:  Papers (6)
| | PDF (277 KB)

The probability of single-longitudinal-mode (SLM) operation is calculated for three kinds of DFB facet structures. To obtain high SLM probability, a DFB laser with cleaved facets requires relatively large kL, but the introduction of asymmetry to facet reflectivity greatly relaxes the kL requirements. View full abstract»

• ### Varactor diode for tuning high-power solid-state microwave oscillators

Publication Year: 1985, Page(s):1152 - 1153
| | PDF (257 KB)

A novel varactor diode has been developed for tuning high-power microwave solid-state oscillators. The resulting output power against frequency remains much flatter than conventional high-power varactor-tuned oscillators. View full abstract»

• ### Double-resonance technique for C/V measurements of semiconductor devices

Publication Year: 1985, Page(s):1153 - 1154
| | PDF (269 KB)

A double-resonance technique has been developed for measuring the capacitance/voltage characteristics of junction semiconductor devices. It can be used in regions of forward bias where currently used methods are often unreliable. View full abstract»

• ### High reliability of 1.5 μm two-step VPE-LPE-grown DC-PBH lasers

Publication Year: 1985, Page(s):1154 - 1156
Cited by:  Papers (1)
| | PDF (389 KB)

The reliability of 1.5 μm InGaAsP/InP DC-PBH lasers prepared by two-step hydride VPE and LPE was estimated at the 70°C, 5 mW constant-power mode. The low degradation rate (8Ã10-5/h) indicated high reliability of the two-step VPE-LPE-grown LDs with the VPE-grown active layer. View full abstract»

• ### Switched-capacitor integrators with low finite-gain sensitivity

Publication Year: 1985, Page(s):1156 - 1157
Cited by:  Papers (42)
| | PDF (357 KB)

Compensated SC integrator stages are analysed and compared with respect to their sensitivities to finite amplifier gain effects. Simple single-stage amplifiers may be used in such circuits even in highly selective filtering applications. This may allow an extension of the present frequency range of SC circuits. View full abstract»

## Aims & Scope

Electronics Letters is internationally renowned for its rapid communication of new developments and emerging topics across the broad and interdisciplinary field of modern electronics and electrical engineering.

Full Aims & Scope

## Meet Our Editors

Editors-in-Chief

Professor Chris Toumazou
Director of the Biomedical Institute
Imperial College, London, UK

Professor Ian White


text:
# Should I define default register values in start up code?

It has come to my attention that some ARM programmers set some CPU registers with values that are the same default values of the register on microcontroller reset.

For example: the default value of IOCON_PIO2_9, parameter MODE on LCP1114 is 0x02 (pull-up resistor enabled). In Olimex's example code this value is set to 0x02 in software initialization but this is already the default value.

This behavior occurs with many other registers.

Why would someone do that? In my opinion, this extra code just pollutes everything and make code readability worse. Is this a good programming practice?

-
I'm not a person to be talking about "best practice", but especially on a platform like ARM, where all peripheral functionality is on the manufacturer and you might find yourself jumping between manufacturers, it seems like this practice would enhance portability. –  Scott Seidman Jan 20 at 14:11
You ever heard the quote "Assumption is the mother of all f***kups"? Never assume anything, even the default values can change if there is a brownout or somesuch so never rely on anything to be in any state unless you set it yourself. –  John U Jan 20 at 19:07

It may be worth noting that if FOO->CTL has a default value of 0x12345678, and one only cares about having the lower bits set to 8, there's a clarity/code-space/performance trade-off between writing the entire register to 0x12345678, versus using something like maskedwrite32(&(FOO->CTL), BLAH_MASK, BLAH_MODE8);. If the masked-write routine includes code to make it atomic in case of interrupts, that may make it slower than it needs to be, but improve the code's generality. –  supercat Jan 20 at 17:27


text:
# Identifying Links Between the S&P500 and VIX Derivatives

Posted on 4/5/16 in Research Articles

#### By Andrew Papanicolaou

The technique of volatility trading has been common practice since the 1970′s. Typically, a long (short) position in volatility included a long (short) position in options. In 2003 there came a more standardized way of trading volatility, as the VIX formula became the universally-accepted predictor of volatility. The VIX is a 30-day predictor of volatility given by a weighted portfolio of out-of-the-money European options on the S&P500 (see Demeterfi, et al. 1999),

$$VIX= \sqrt{\frac{2e^{r \,\! \tau}}{\tau} \left(\int_{0}^{F} \frac{P(K)}{K^2} dK + \int_F^{\infty} \frac{C(K)}{K^2} dK \right)}$$

where $$τ = 30$$ days, $$r ≥ 0$$ is the risk-free rate, F is the 30-day forward price on the S&P500, and with $$P(K)$$ and $$C(K)$$ being respective prices for puts and calls with strike and 30 days to maturity. An important thing to notice in the VIX formula is how it’s composed entirely of option premia (i.e. out-of-the-money options have no intrinsic value) and hence trading VIX is like trading the premium that is priced into the options market. A common practice is to hedge high volatility with long positions in derivatives on the VIX itself. For example, the Chicago Board Options Exchange uses long positions in VIX call options to produce their VIX Tail Hedge Index, an index which historically has performed better than several benchmarks (Chicago Board Options Exchange (CBOE) 2011).

The VIX was among the topics discussed at IPAM’s Spring 2015 program “Broad Perspectives and New Directions in Financial Mathematics”, with the goal of finding ways to bring greater stability into banking systems, algorithmic trading, and commodities markets. There are countless models for pricing and hedging of derivative products, but for IPAM participants it was the model parameters implied by market prices that were of primary interest. In particular, VIX is an implied volatility parameter that quantifies fear/optimism among traders, and in the past has been a strong indication of any stress or impending crisis.

A contemporary problem in this area during Spring of 2015 was in the establishment of links between markets for VIX derivatives and the primary securities (i.e. no VIX products; only S&P500, S&P500 options, low-risk bonds, etc….) (Carr and Wu 2006). Relative to the market sub-space spanned by primary securities, VIX derivatives provide a new dimension of hedgability because buying and selling of VIX previously was not possible. Figure 1 illustrates the added dimension relative to the sub-space of primary securities. A new result that was produced at IPAM was the linking of the two markets with the so-called time-spread portfolio –a portfolio that is long in the spreads between 1-month and 2-month S&P500 put options (Papanicolaou 2016). The time-spread portfolio is part of a term-structure on volatility composed entirely with primary market securities, and is a non-null component of the primary market sub-space that is also present in VIX derivatives (see the red line in Figure 1). Hence, the time-spread portfolio links the two markets by providing a robust/model-free component that is common to both market subspaces.

Figure 1: The horizontal place is the market sub-space spanned by primary securities, and the vertical direction is the market sub-space created by VIX options. The red line represents the projection of the VIX-options component onto the primary sub-space.

The 2015 program was a movement toward a deeper analysis of financial data and making new discoveries about information implied by financial markets. Thanks to the efforts of the program’s participants, the future of the financial system will be less uncertain with an increased understanding of market dynamics and financial risk.

#### Works Cited

Carr, P., and L. Wu. “A Tale of Two Indices.” The Journal of Derivatives 13, no. 3 (2006): 13-29.

Chicago Board Options Exchange (CBOE). White paper: CBOE VIX tail hedge index (ticker VXTH). Available at http://www.cboe.com/micro/vxth/ documents/vxthwhitepaper.pdf., 2011.

Demeterfi, K., E. Derman, M. Kamal, and J. Zou. “More Then You Ever Wanted to Know About Volatility Swaps.” The Journal of Deriviatives 6, no. 4 (1999): 9-32.

Papanicolaou, A. “Analysis of VIX Markets with a Time-Spread Portfolio.” (under review with journal) (available at http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2618933.), 2016: 1-35.


text:
# List of combinatorial problems

I have a list of problems that I am currently thinking about. It’s quite difficult to judge which one is easy or whether a conjecture is true. But all of them seems to be important enough, in the sense of the possibility to initiate a whole new branch of combinatorics. I also mark their importance (in my opinion) by stars. Certainly there are many other important open problems, here is just a very small collection of them, those I spent at least a couple of hours on.

A. Graph coloring:

1. Erdos-Faber-Lovasz conjecture  ✭✭✭✭
2. Hedetniemi conjecture ✭✭✭✭
3. Chromatic number of Cayley graph ✭✭✭✭
4. Doubly critical graph conjecture ✭✭

Common points: Chromatic number is NP-hard, but for graphs with nice algebraic structures, maybe we can characterize it. It’s also possible to find some other also-NP-hard equivalent description of chromatic number, which is easier to study for specific problems.

1. Seymour’s second neighborhood conjecture ✭✭✭✭

C. Correlation Inequality:

1. Bunkbed conjecture ✭✭✭✭

D. Extremal Combinatorics:

1. Generalization of Sperner Lemma ✭✭✭

E. Combinatorial Geometry:

1. Unit distance problem ✭✭✭✭
2. Kusner conjecture (l1 distance) ✭✭✭✭

# News

I guess now I know how to prove generalized bunkbed conjecture when there are two vertical edges. This is the first nontrivial case since when there is only one vertical edge, the statement follows directly from FKG inequality about increasing events.

I will not uncover more details until I found a way to solve the general case (when there are more than two vertical edges), or I decide to give up on this conjecture and submit a junk paper.

Anyway, I feel quite happy since I am finally able to prove something, rather than giving counterexamples…

# Ahlswede-Daykin Inequality

Is there any good generalization of Ahlswede-Daykin inequality?

There is a paper by R. Aharoni. But I believe A-D inequality (or FKG) definitely have many different versions of generalizations. One of them can possibly lead to the proof of Bunkbed Conjecture.

# How to give a good talk

Actually I am only a third-year graduate student now and have very little experience to give a public talk. So here I basically only collect what other people said and wish that I can learn something from them.

1. Presentation

If you are  happened to talk in SIAM or other place, this will be very useful:
How to give a good 20-minute math talk

2. Preparing for the slides

# My mathematical diseases

Richard Lipton has an article on mathematical disease. My own definition of mathematical disease is quite similar: it must be simple in form, easily explainable to a non-specialist in a few minutes. Moreover, there should be more than one person exclaiming “how can such a problem be open?”. Here I also list some problems I consider as my favourite “diseases”.

1. Erdos-Faber-Lovasz Conjecture

The main difficulty of this conjecture is, how to bound the chromatic number from above effectively? There are some lower bounds for chromatic number, like using eigenvalues of adjacency matrix, or topological bounds developed by Lovasz in proof of Kneser Conjecture. However for upper bound, we knew very little except giving an explicit coloring. Although an asymptotic version of this conjecture has been proven by J.Kahn, I believe the exact version, if indeed correct, requires more understanding about chromatic number.

2. Hedetniemi Conjecture

A reformulation of this conjecture can reduce it to proving the chromatic number of some graphs can be bounded from above. So maybe it has some connection with the previous problem.

3. Bunkbed Conjecture

There is no article on this conjecture in wikipedia on mathworld. But it is very easy in appearance, and once proved, it can be written as a theorem in any textbook of probability in graduate level. The most general form is as follows, given a graph $G$, consider its Cartesian product with complete graph of two vertices $K_2$. Now consider the model such that every edge in this graph appears independetly in some probability. The only restriction is the probability of any edge $uv$ in $G$ and that of its corresponding edge $u'v'$ are equal. The conjectures says the probability that there is a path connecting $u$ and $v$ is greater than or equal to the probability that there is a path connecting $u$ and $v'$.

4. Unit distance problem

More precisely, I am talking about the problem of counting the maximum number of pairs of unit distances for $n$ points in $\mathbb{R}^2$. The other, perhaps more famous, problem of computing chromatic number of unit distance graph of $\mathbb{R}^2$ seems much more difficult.

5. Kusner Conjecture

I also don’t find a survey paper on this conjecture. It’s so natural and beautiful that every middle school student can have some thought on it. The conjecture is: consider the $l^1$ metric (sum of absolute values) in the d-dimensional real space $\mathbb{R}^d$, there can be at most $2d$ points whose pairwise distances are exactly 1. This bound can be reached by looking at the octahedron with vertcies $(\frac{1}{2},0,0), (-\frac{1}{2},0,0), (0,\frac{1}{2},0), (0,-\frac{1}{2},0), (0,0,\frac{1}{2}), (0,0,-\frac{1}{2})$.

# A digraph whose longest directed path<=k is k-colorable

To prove this theorem, we have the following steps: (first I will assume the graph G is weakly connected)

(i) Find a maximal subgraph H of G which contains no directed cycles.
(ii) Prove that this subgraph H can be k-colored such that if v->u, then the color number c(v) is strictly less than c(u).
(iii) If we add any directed edges into H, it must result a cycle. Use (ii) to show the two end vertices of this edge has different colors.

The main difficulty lies in step (ii). Let’s also divide it into several small steps: (remember that now our graph contains no directed cycles.)

(ii.1) Find a longest directed path (its length is k), choose its starting vertex v. Then for every vertex u which can be reached from v, give it the color c(u)= maximal length of directed path from v to u.
(ii.2) Prove that for every u1, u2 s.t. c(u1)=c(u2), there is no edge between them. Therefore, so far our coloring is proper.
(ii.3) From the remaining vertices, find a vertex w such that it maximize the following thing:
the maximal length of paths from w which doesn’t touch any colored vertices except the last one – the color number of the last vertex.
Then for any uncolored x, let c(x)=maximal length of paths from w to x which doesn’t touch any colored vertices. Like in step (ii.1) and (ii.2), we can also prove that so far the coloring is proper.
(ii.4) Continue this process, since G is weakly connected. We will finally get a proper coloring of all the vertices.