text:
## The curse of dimensionality and finite vs. asymptotic convergence results

Related to our (Aki, Andrew, Jonah) Pareto smoothed importance sampling paper I (Aki) received a few times a comment that why bother with Pareto smoothing as you can always choose the proposal distribution so that importance ratios are bounded and then central limit theorem holds. The curse of dimensionality here is that the papers they refer used small dimensional experiments and the results do not work so well in high dimensions. Readers of this blog should not be surprised that things look not the same in high dimensions. In high dimensions the probability mass is far from the mode. It’s spread thinly in surface of high-dimensional sphere. See, e.g. Mike’s paper Bob’s case study, and blog post

In importance sampling one working solution in low dimensions is to use mixture of two proposals. One component tries to match the mode, and the other takes care that tails go down slower than the tails of the target ensuring bounded ratios. In the following I only look at the behavior with one component which has thicker tail and thus importance ratios are bounded (but I have made the similar experiment with mixture, too).

The target distribution is multidimensional normal distribution with zero mean and unit covariance matrix. In the first case the proposal distribution is also normal, but with scale 1.1 in each dimension. The scale is just slightly larger than for the proposal and we are often lucky if we can guess the scale of proposal with 10% accuracy. I draw 100000 draws from the proposal distribution.

The following figure shows when the number of dimensions go from 1 to 1024.

The upper subplot shows the estimated effective sample size. By D=512 importance weighted 100000 draws have only a few practically non-zero weights. The middle subplot shows the convergence rate compared to independent sampling, ie, how fast the variance goes down. By D=1024 the convergence rate has dramatically dropped and getting any improvement in the accuracy requires more and more draws. The bottom subplot shows Pareto khat diagnostic (see the paper for details). Dashed line is k=0.5, which the limit for variance being finite and dotted line is our suggestion for practically useful performance when using PSIS. But how can khat be larger than 0.5 when we have bounded weights! Central limit theorem has not failed here, but we have just not reach yet the asymptotic regime to get CLT to kick in!

Since humans are lousy in looking at 1024 dimensional plots the top subplot shows the 1 dimensional marginal density of the target (blue) and the proposal (red) densities of the distance from the origo r=sqrt(sum_{d=1}^D x_d^2). The proposal density has only 1.1 larger scale than the target, but most of the draws from the proposal are away from the typical set of the target! The vertical dashed line shows 1e-6 quantile of the proposal, ie when we draw 100000 draws, 90% of time we don’t get any draws from there. The middle subplot shows the importance ratio function, and we can see that the highest value is at 0, but that value is larger than 2*10^42! That’s a big number. The bottom subplot scales the y axis show that we see importance ratios near that 1e-6 quantile. Check the y-axis: it’s still from 0 to 1e6. So if we are lucky we may get a draw below the dashed line, but then it’s likely to get all the weight. The importance function is practically as steep everywhere where we can get draws in a time of the age of the universe. 1e-80 quantile is at 21.5 (1e80 is the estimated number of atoms in the visible universe). and it’s still far away from the region where the boundedness of the importance ratio function starts to affect.

I have more similar plots with thick tailed Student’s t, mixture of proposals etc, but I’ll save you from more plots. As long as there is some difference in target and proposal taking the number of dimensions to high enough, IS and PSIS break (PSIS giving slight improvement in the performance, and more importantly can diagnose the problem and improves the Monte Carlo estimate).

In addition that we need to take into account that many methods which work in small dimensions can break in high dimensions, we need to focus more on finite case performance. As seen here it doesn’t help us that CLT holds if we never can reach that asymptotic regime (same as why Metropolis algorithm in high dimensions may require close to infinite time to produce useful results). Pareto diagnostics has been empirically shown to provide very good finite case convergence rate estimates which also match some theoretical bounds.

PS. There has been lot of discussion in comments about typical set vs. high probability set. In the end Carlos Ungil wrote

Your blog post seems to work equally well if you change
“most of the draws from the proposal are away from the typical set of the target!”
to
“most of the draws from the proposal are away from the high probability region of the target!”

I disagree with this, and to show evidence for this I add here couple more plots I didn’t include before (I hope that this blog post will not have eventually as many plots as the PSIS paper).

If the target is multivariate normal, we get bounded weights by using Student-t distribution with finite degrees of freedom nu, even if the scale of the Student-t is smaller than the scale of the target distribution. In the next example the target is same as above, and the proposal distribution is multivariate Student-t with degrees of freedom nu=7 and variance=1.

The following figure shows when the number of dimensions go from 1 to 1024.

The upper subplot shows the estimated effective sample size. By D=64 importance weighted 100000 draws have only a few practically non-zero weights. The middle subplot shows the convergence rate compared to independent sampling, ie, how fast the variance goes down. By D=256 the convergence rate has dramatically dropped and getting any improvement in the accuracy requires more and more draws. The bottom subplot shows Pareto khat diagnostic which predicts well the finite case convergence rate (even if asymptotically with bounded ratios k<1/2).

The top subplot shows the 1 dimensional marginal density of the target (blue) and the proposal (red) densities of the distance from the origo r=sqrt(sum_{d=1}^D x_d^2). The proposal density has the same variance but thicker tail. Most of the draws from the proposal are away from the typical set of the target and in this towards higher densities than the density in the typical set. The middle subplot shows the importance ratio function, and we can see that the highest value is close to 47.5 and then the ratio function starts to decrease again. The highest value of the ratio function is larger than 10^158. The region with highest value is far away from the typical set of the proposal distribution. The bottom subplot scales the y axis show that we see importance ratios near the proposal and the target distributions. The importance ratio goes from very small values up to 10^30 in very narrow range, and thus it’s likely that the largest draw from the proposal will get all the weight. It’s unlikely that in practical time we would get enough draws to get the asymptotic benefits of bounded ratios to kick in.

PPS. The purpose of this post was to illustrate that bounded ratios and asymptotic convergence results are not enough for the practically useful performance for IS, but there are also special cases where due to special structure of the posterior we can get practically good performance with IS, and especially with PSIS, also in high dimensions cases (PSIS paper has 400 dimensional example with khat<0.7).

1. ojm says:

> most of the draws from the proposal are away from the typical set of the target!

How are you defining typical set here?

• Intuitively, it’s where random samples from the posterior will land. The discrete case is covered in the Wikipedia page:

https://en.wikipedia.org/wiki/Typical_set

Cover and Thomas’s book cover the continuous case with differential entropy. I cover the basic idea in my case study with examples (borrowed from McKay’s great info theory book, which is also available free online in PDF form with permission of the publisher):

http://mc-stan.org/users/documentation/case-studies/curse-dims.html

• ojm says:

Hi Bob,

I asked because I was one of the people who felt your informal usage of typical set was incorrect in the technical sense. See also our discussion on your blog post

• ojm says:

For example, why not just say ‘high probability set’? Or do you actually mean ‘typical set’? If so, why the typical set and not the high probability set?

• Aki Vehtari says:

I’m not an expert on typical sets and high probability sets, but inspired by this question I wanted to learn what is the difference between these. I didn’t yet find an answer in continuous case, but with a quick google search I found this http://mat.hjg.com.ar/tic/img/lecture4.pdf, which says

To contrast, consider X_i be Bernoulli(.9). Then a typical set contains sequences in which the proportion of 1’s is close to .9. But it does not contain the single most probable sequence, the sequence of all 1s. The set B _\delta^{(n)} (the high probability set) will include all the most probable sequences, including the sequence of all 1s.

If the sequence is very long (ie X_i is high dimensional) it would be very unlikely to observe sequence of all 1s even if it is the single most probable sequence. In continuous case in high dimensions it’s unlikely we observe anything near the mode. So it seems typical set is more appropriate term.

If you have a better reference for typical set vs high probability set, please post in a comment.

• ojm says:

Bob’s references, in particular Cover and Thomas and Mackay are fine, and discuss the difference. The high probability set (of level alpha) is simple to define: the smallest set containing at least alpha of probability. The typical set is more convoluted and relies on the AEP to make sense of imo.

I think the comparison you give is somewhat misleading: you are comparing a single realisation to a set of realisations. The same argument can be made against any other single observation, however this would be even more unlikely than the mode point (in the continuous case all points have zero probability anyway, but let’s consider the discrete case for concreteness).

By definition any single observation has a higher probability of landing in the high probability set than in the typical set, though the difference is asymptotically negligible, as shown by the AEP.

• ojm says:

The most convincing argument for the typical set that I can think of is essentially the same as that given by Mackay: it is a convenient approximation to the high probability set.

Which is fine, but this undermines the selling point of the ‘counterintuitiveness’ of the typical set vs the high probability set. The mode is in the high probability set and so isn’t inherently ‘bad’ as seems to be the implication of some of Bob’s posts etc. It’s just that one can throw out the mode point and still get a set of approximately alpha of probability.

You can also throw out any other point and still get a set of approximately alpha of probability – in fact you will get one of (very slightly) higher probability.

• Dan Simpson says:

My view on this is that typical set is good because:
– it’s common English meaning is useful to people who don’t know information theory
– it is the set we want (or very close to it)
– it’s definition actually matches what we want (actually for LOO we want that the log of a nonlinear function of the sample converges to almost the right thing, but this is “close enough)
– it emphasises the link between the typical set and the small set you need for geometric ergodic. Intuitively they’re complements.

You’re argument appears to be that you think we’re not using the technical definition. I don’t agree with you. If your problem is communication, I really think typical set is unambiguous to both technical and non-technical communities.

• ojm says:

> it’s common English meaning

What’s wrong with ‘high probability set’?

> it is the set we want (or very close to it)

What’s wrong with ‘high probability set’?

> it’s definition actually matches what we want

What’s wrong with ‘high probability set’?

etc.

• ojm says:

Or do you actually want “a set of sequences whose probability is close to two raised to the negative power of the entropy of their source distribution”? I.e. the typical set?

• Dan Simpson says:

I’m not sure why you think repeating your self is making your point better. I prefer typical set for the reasons set out above.

• Carlos Ungil says:

> What’s wrong with ‘high probability set’?

Recently I read a footnote (https://www.aqr.com/cliffs-perspective/wild-but-not-crazy) giving the following reason for the use of continuously compounded returns: “They are a bit better behaved than simple compound returns. Plus, I always feel smarter using logs.”

In the “typical set” vs “high probability set” choice I still don’t get what is the objective reason for preferring the former, given that in practice the refer to the same set.

• ojm says:

Carlos: I’m glad you comment around here!

• ojm says:

Dan – when I feel people don’t address my question clearly I’m not sure what else to do but ask it again. I think I’ll give up on arguing with the Stan folk about this, as the typical responses seem unhelpful with high probability

• Dan Simpson says:

If you don’t want to use “typical set”, don’t. I don’t like using the word “cognisant” because I prefer a synonym. I don’t feel the urge to interrogate people as to why they use it.

• Dan Simpson says:

Ojm: I literally gave you four reasons (the last two very specifically using the technical definition). I’m done with this now.

• ojm, Carlos: i’m glad you both stick around..

ojm: I think it’s an interesting fact that essentially all of the probability mass in a high dimensional posterior is in the set where log(p(x)) is within epsilon of some special value lp*

If we just *define* the word “typical set” to mean “the set of x values such that log(p(x)) = lp* +- epsilon”

and we define the high probability mass (HPM) set to mean “the smallest set such that integrate(p(x)) = 1-alpha”

then it’s an interesting fact that the two sets both have basically all of the probability mass in them. It seems the main difference between them is “on the boundary” where the HPM set includes the mode, and the typical set doesn’t include the mode and as a consequence includes a negligibly wider bit of the “tail” than the HPM set.

From a practical perspective it seems to me that arguing over these two is like arguing whether we should include a particular hair on your head, and a particular sliver of your fingernail as “part of your body”. Some people prefer including the head and leaving out the sliver of fingernail, some people prefer to include the fingernail, and leave out some of the hair… either way almost every portion of your body is included, so for people who want to calculate the mass of your liver or the volume of your blood or other “functionals of the set” they get the same answer regardless.

• ojm: my biggest problem is that the *usual* definition of typical set is about *sequences* of independent draws from a single univariate generating process (thereby producing a high dimensional product space). However in Stan/Bayes we’re interested in *points from a high dimensional probability space* which are not sequences of independent draws from a univariate generator. The fact about the log(p(x)) still being nearly constant is what now defines the subset of the space… so when you look up in wikipedia or elsewhere, it’s hard to map technical definitions as used in info theory to technical definitions as used by say Bob Carpenter in his case studies.

One “solution” to this communication problem would be to explicitly define the technical definition of the word “typical set” being used in the case studies.

• ojm says:

DL – true, it is splitting hairs somewhat. But aomw arguments seem to be based on this hair split. For example the Stan folk consistently maintain that they are in fact interested in the typical set rather than the high probability set. So it seems worth examining.

• ojm: If we want to get into the hair-split, the first thing we need is a proper definition of the typical set that applies to sampling from a high dimensional Bayesian posterior rather than from *sequences of independent draws from univariate distributions*. The sequences definition is very useful in information theory because there we are talking about communication across a channel… a sequence of symbols is how we communicate information.

In the Stan case, it’s very easy to see, practically speaking, that once the HMC process has settled into stationarity, all the samples are within epsilon of having the same value for log(p(x)). Basically, all the probability mass is in this constant lp region. It’s also all in the HPM region as well, it’s just that these sets overlap on a set of nearly measure 1, and the regions where they don’t overlap has very very small total measure, either because the volume is small (near the mode), or the density falls off dramatically quickly (towards the “tail”).

It seems to me that the reason the Stan crew focuses on the typical set (as defined by near constant log(p)) is that log(p(x)) is something they’re calculating all the time, and the samples all have nearly the same value, and this just happens because of the geometry of high dimensions, not because they have to “avoid” the mode intentionally. When you’re calculating a thing all the time, and using it to determine whether your algorithm is working properly, it’s more obvious to focus on that aspect.

Basically, if you see a bunch of samples that come from a region with dramatically different log(p(x)) it’s an indicator that either you haven’t converged to stationarity of the sampler, or there’s a bug in the code. You shouldn’t ever see a sample near the mode, simply because the volume of the region near the mode is so small you effectively can’t hit it, because you’re starting each iteration with a random velocity vector. Similarly, you shouldn’t see samples with lp much less (more negative) than the typical set value because even though the volume of space with values like that is effectively infinite, the “energy” required to get there is prohibitively large.

So, I think the fact that nearly all the probability mass is in the typical set, combined with the fact that Stan explicitly calculates log(p(x)) at every iteration and spits it out in the output… and the fact that naturally through correct sampling, in high dimensions, essentially every sample has a constant value for log(p(x))… leads to a focus on that set. Rather than any other fundamental reason.

• ojm says:

DL – yeah maybe man. If so I dunno why they don’t just say so though. I’m probably missing something here but I doubt I’m going get anywhere with discussing the nuances here.

• Chris Wilson says:

+1 (to Daniel Lakeland’s comments). Also, I think what Bob said is also easy to lose sight of: “intuitively, it’s where random samples from the posterior will land” (that is, think about this *generatively*).

• Dan Simpson says:

Daniel: yep. That’s basically exactly it (at least from my viewpoint. The Stan team is far from homogeneous and literally no one on earth has ever wanted met to speak for them). Nice explanation.

• ojm: I really think this is it, because when you’re using Stan and/or debugging Stan models, the near-constancy of the lp variable and traceplots of lp are extremely useful for figuring out how well Stan is doing at sampling. You know it’s not working if lp is still trending or oscillating or stuck somewhere and not changing at all.

But from a theoretical perspective, I agree with you that the mode is not “to be avoided” it’s just not very helpful for calculating expectations, because there’s a LOT of mass elsewhere you need to take into account to get good expectation calculations.

• Carlos Ungil says:

> If the sequence is very long (ie X_i is high dimensional) it would be very unlikely to observe sequence of all 1s even if it is the single most probable sequence. In continuous case in high dimensions it’s unlikely we observe anything near the mode. So it seems typical set is more appropriate term.

It is even more unlikely that we observe anything near any other given point.

When we look a a sequence of 0/1 bits from a one-dimensional binomial distribution, the concept of a typical sequence may be interesting because the permutations are interchangeable. For example, we are often interested in the proportion of 1s and 0s.

When we look at a single realization of a highly-dimensional distribution (corresponding to coordinates in the sample space) what sense does the idea of a typical sequence make? Unless all the coordinates are independently sampled from the same distribution, we cannot say that they are a sequence of realizations of a variable.

And even if we accept that the coordinates form some kind of sequence, the permutation of coordinates results in “similar” points only if the problem (probability distribution and function of interest) has spherical symmetry. If we calculate the expected value of an arbitrary function f(x1,..,x1000) for 1000-dimensional normal distribution two points in the typical set are in no meaningful way closer to each other than they are close to the origin.

2. Aki Vehtari says:

Based on the discussion here I think I have understood distinction between typical set and high probability set (at least approximately correctly) correctly, but let’s check that.

ojm wrote

The mode is in the high probability set and so isn’t inherently ‘bad’ as seems to be the implication of some of Bob’s posts etc.

Important part of my blog post was that in order to get in the asymptotic regime of bounded ratio benefits requires that we get draws from the proposal distribution which are close to the the mode. Close to the mode being maybe inside sphere with radius 5. I mentioned that probability mass close to the mode inside a sphere with radius approximately 21.5 is less than one per the number of atoms in the observable universe. For me it seems unlikely that with current computers we would get draws inside this sphere (r approx 21.5) in my lifetime so I would call that untypical set. If high probability set includes (in probability, epsilon etc.) inside of that sphere (r approx 21.5) and typical set does not include inside of that sphere (r approx 21.5) then I think typical set is more useful for me to describe what happens in finite case. For the theoretical (but completely impractical) asymptotic case I might be fine with high probability set, too.

• ojm says:

Pick a set of fixed volume. Then the high probability set is the one of that volume that you are most likely to see in your lifetime. When you compare sets of different volume then the comparison is not ‘fair’ so to speak. If there was no radial symmetry then would you lump all points of a fixed radius away together?

• If we pick a set of fixed volume large enough that its total probability is not negligible, then the high probability set with that volume will overlap the typical set, and if we sample until we hit this volume, the sample point will probably be in the typical set “portion” of this set just because the total fraction of the volume that is in the typical st portion is so much larger than the fraction that is near the mode.

In high dimensions, the HP set is basically equal to the typical set, plus a tiny volume near the mode, minus a miniscule wafer thin shell away from the mode. The euclidean distance to the mode from any point in the typical set can be quite large, and yet it doesn’t matter because the volume near the mode is so small, we can go “deep” into this “Modal tail” without adding much probability, hence we only need to shave a razor thin shell off the “outside”

• Chris Wilson says:

If you think a bit generatively – about drawing samples from the posterior distribution – it is clear that in high dimensions you will never draw samples near the mode. As you say, this is not due to any deliberate “modal exclusion”. Thus, algorithms to use MCMC to draw posterior samples (or compute expectations wrt posterior, as some say), should be guided by geometric information that aids in discovering the ‘typical set’, rather than this ‘HP set’, although their overlap is quite substantial.

• ojm says:

DL – yes, I agree. But note also that by talking about points in a shell a given distance away you are talking about a set.

As soon as you go to the second part of what you are saying you are dropping the constraint of comparing sets of fixed volume. Whenever you fix either total probability or volume, the high probability set is natural. Every argument for the typical set after this drops the fixed constraint and hence makes comparisons somewhat arbitrary. I would like to see a single example where the volume or probability is held fixed that makes the typical set more natural.

Of course the volume of a shell of given radius is larger in high dimensions than the small ball around the mode. But that’s not the point. (This is where eg Bob jumps in and says ‘that’s the point’ – but no, it is not my point.
Then Dan jumps in and says ‘well what is your point?’. Then I say something snarky and we unfollow each other on Twitter, and all is well with the world).

• Dan Simpson says:

I will buy you a beer and we can chat some day if we’re ever in the same place, but the internet is clearly not a medium that we can successfully communicate on.

• ojm says:

Or, how about an example without spherical symmetry? The symmetry makes it natural to think about the shell as ‘one thing’. But if you compare a small fixed segment of the shell to a volume around the mode of the same volume then the points are more likely to land around the mode than in the given segment, ie in the portion of the high probability set not the portion of the typical set.

• Chris Wilson says:

ojm, well what IS your point? :) In the bigger sense, isn’t the point to evaluate probability measures (requiring integration of densities over volumes)?

• ojm says:

Where it is stated

> Why ensemble methods fail: Executive summary

1. We want to draw a sample from the typical set
2. The typical set is a thin shell a fixed radius from the mode in a multivariate normal
3. Interpolating or extrapolating two points in this shell is unlikely to fall in this shell
4. The only steps that get accepted will be near one of the starting points
5. The samplers devolve to a random walk with poorly biased choice of direction

Without disagreeing with 2-5, I question the premise 1.

I think one reason I push back is a sense (there and elsewhere) of an emphasis on ‘wow counterintuitive result – the mode is bad’. But this is not so counterintuitive in the sense that the mode is no worse than any other single point. If the mode is included the appropriate number of times, along with other points being included an appropriate number of times, it’s fine. So it seems disingenuous to use that as premise one in criticisms of other methods – even if there are completely legitimate other criticisms to make!

• ojm says:

I think the pushback plays a role too – if people just said ‘yeah we want the high probability set but the typical set is convenient for us’ then cool. But when people actually try to make the argument that we really truly want the typical set and not the high probability set, I dunno, I guess it just provokes me to pointlessly argue on the internet about it.

• “If the mode is included the appropriate number of times, along with other points being included an appropriate number of times, it’s fine.”

But, due to the probability mass in the vicinity of the mode being nearly zero, the appropriate number of times to sample near the mode in any reasonable sample size… like thousands to millions of samples… is *zero*

whereas the appropriate number of times to have samples whose log(p(x)) = lp* is basically N, the number of samples you’re drawing. Sure, any *given* point in the typical set shouldn’t appear, but the fraction of points which are in the typical set should be essentially all of them.

• Carlos Ungil says:

> the fraction of points which are in the typical set should be essentially all of them.

The fraction of points which are in the high density region is also esentially all of them, so if there is a good reason to prefer the typical set it may be somewhere else :-)

• “The fraction of points which are in the high density region is also esentially all of them”

No, that’s not true. In the spherical symmetry normal distribution case, the highest density region is the mode near x=0, let’s say with radius 1 or so, but in a 1000 dimensional normal, the radius of such randomly selected points has a chi distribution with 1000 degrees of freedom and in such a distribution 99% of the points are outside radius 30 or so.

so points are not near the “high density” region, but they are in the “high probability” region because the “high probability” region includes everything from radius 0 to radius say 36, however almost all the points are radius 30-36 not anywhere “near” radius 0

• Carlos Ungil says:

I use High Density Region and High Probability Region as synonyms: the region where the probability density function is above some threshold. For a fixed volume or probability mass, what difference do you see between the HDR and the HPR?

• Carlos, take the 1000 dimensional unit normal. I do

vec = rnorm(1000)
sum(sapply(vec,function(x),dnorm(x,log=TRUE)))
and get
-1411.434

and sqrt(sum(vec*vec)) = 31.38

now I investigate the density near the mode:

I do sum(sapply(rep(0,1000),function(x) dnorm(x,log=TRUE)))
and get
-918.9385

so the density near my randomly chosen point, which is in the typical set, is exp(-492.5) times less dense than the density near the mode. So I don’t see how we can call that a high *density*

On the other hand, the probability, which involves integrating density * volume is very large for shells of radius 30-35 because the volume of those shells is enormous, not because the density is anything much more than 0 to 200 decimal places.

• Carlos Ungil says:

Daniel, what is the 95% HDR in your example?

• The usual sense of 95% region (trimming on both ends) we can get by

sqrt(qchisq(.125,1000)), sqrt(qchisq(.975,1000))

which is the shells for radius between 30.8 and 33.008

But ojm’s set, is the ball of radius less than

sqrt(qchisq(.95,1000)) = 32.7823

so includes the whole region from radius 0 to 32.7823 and this volume is slightly smaller than the volume between 30.8 and 33.008

but in both cases, almost all the points are in the range r = 31 to 32 and essentially none of the points are inside radius 28

sqrt(qchisq(1e-6,1000))
[1] 28.31297

• Carlos Ungil says:

The HD in “95% HDR” means that the points within the region have higher probability density than the points without the region.

The 95% HDR *is* the ball of radius 32.8 (assuming your calculation is correct).

The 95% HDR *is not* the shell between 30.8 and 33.0 (obviously, because the probability density is higher for the points in the inner “hole” than for the points in the shell).

• Carlos Ungil says:

without the region => outside the region

• Ok, fine ojm’s set and your “95 HDR” are both the ball of radius less than 32.8 but note that this is predicated on *both* high density *and* high total probability not just “the region where the probability density function is above some threshold”

If we took the region where probability density in the 1000 dimensional space is above say 0.001 times the peak density… which is a pure density based definition of a set, then this set is the ball of radius less than around 3.7

and this ball is inside your 95% High Probability ball around 0 but *does not intersect the typical set*, and this ball has essentially zero points in it precisely because it doesn’t intersect the typical set.

my calculation:

dnorm(3.7)/dnorm(0)
[1] 0.001064766

pchisq(3.7^2,1000,log.p=TRUE)
[1] -1656.402

so the total probability inside the set of density greater than 0.001 fraction of the peak density is exp(-1656)

• Carlos Ungil says:

> this is predicated on *both* high density

That’s the HD.

> *and* high total probability

That’s the 95%.

> not just “the region where the probability density function is above some threshold”.

More precisely the region where the probability density function is higher than the threshold that makes the probability mass contained in the region 95%.

I’m not sure what do you find so dificult about the concept of the HDR. It is the natural extension of the one-dimensional HDI. I’m sure you know and understand HPD intervals (the shortest credible intervals). How is this different?

• In 1000 dimensions the *geometry* of the HDR is very weird. In the unit normal case we can get a little intuition, but in the general case, it’s like trying to describe the surface of a human brain or something.

I do fully understand when you say the 95% HDR region you are talking about what I and OJM were calling the High Probability region. I guess the point is that when you use the “typical set” or you use the 95% HDR region, essentially all the points are *in both* and so the distinction is mainly occurring in regions that actually don’t matter, like in my example the radius less than 3.7 ball…

In the more general high dimensional posterior geometry, defining the 95%HDR region as “the smallest volume set containing 95% of the probability” is fine, but constructing it is basically impossible, in other words, you can’t compute with it…. But you *can* compute with the typical set because it’s the “set where density is essentially constant and equal to exp(lp*)+-epsilon for some particular lp* and epsilon sufficient to get you 95% of the mass.

HMC naturally is attracted to this set because it does correct unbiased sampling, and essentially all the samples are in this set (which is typically near the boundary of the 95%HPD set)

• Carlos Ungil says:

> HMC naturally is attracted to this set because it does correct unbiased sampling, and essentially all the samples are in this set

It will also be naturally attracted to the HPD set, for the same reason. Actually it will be slightly more attracted to the HPD set, given than it is slightly more compact than the typical set (for the same probability mass).

> (which is typically near the boundary of the 95%HPD set)

In high dimensions, the whole HPD set is also near its boundary.

• Carlos Ungil says:

“set where density is essentially constant”

*constant* as in “ranging over ~35 orders of magnitude” (for the 95% set in the 1000-dimensional standard gaussian).

• “As soon as you go to the second part of what you are saying you are dropping the constraint of comparing sets of fixed volume”

I wasn’t intending that. Suppose we start with the typical set and it’s volume is V.

Now, we want to compare it to the high probability set of *equal volume*. Let’s also just for ease of discussion stick to the independent symmetric normal…

So, the typical set is basically a shell at radius r with half-width w.

The high probability set is the entire ball within radius r+w-epsilon where we choose epsilon so that now the volume we added inside the ball is removed from the exterior “skin” retaining the constant volume V. epsilon is a *very small* number because there is almost no volume inside the ball.

Now the two sets have the same volume by construction, the total probability is basically 1 for both of them to zeroth order, but to first order, the HP set has an epsilon of higher probability. Still, if we randomly sprinkle points inside the HP set according to the probability mass, *essentially all of those points are in the portion of the HP set that intersects the original Typical Set*

If we take the sprinkled points and select from them uniformly at random, with a reasonable finite sample size, say N=100000 then *every one of them will be in the Typical Set portion* virtually guaranteed.

So if we have an algorithm that tends to bias us towards moving in the direction of the mode, to get a correct sample we need to actually *undo this bias* or we will over-sample the mode.

• Carlos Ungil says:

> If we take the sprinkled points and select from them uniformly at random, with a reasonable finite sample size, say N=100000 then *every one of them will be in the Typical Set portion* virtually guaranteed.

Every one of them will be also in the High Density Region portion, virtually guaranteed. And if you calculate how many points can you expect to miss (which we agree will be a very small quantity), the typical set will get worse coverage than the HDR by construction.

> So if we have an algorithm that tends to bias us towards moving in the direction of the mode, to get a correct sample we need to actually *undo this bias* or we will over-sample the mode.

If we had an algorithm that tended to bias us towards moving far from the mode we would also need to undo that bias to get a correct sample. Using unbiased algorithms is probably better :-)

• “If we had an algorithm that tended to bias us towards moving far from the mode we would also need to undo that bias to get a correct sample.”

Yes, but the point of the original post that brought ojms concern to the forefront was that people were proposing algorithms such as sampling along a line between two points which did bias things towards moving outside the typical set (and in the unit normal high dim case, towards the mode), and the metropolis hastings correction required to get unbiased samples caused the algorithms to devolve to tiny-step-size diffusions along the “surfaces” defined by the typical set, and hence they were unhelpful.

• Carlos Ungil says:

> algorithms such as sampling along a line between two points which did bias things towards moving outside the typical set (and in the unit normal high dim case, towards the mode)

Are you sure this bias exists? As the number of dimensions increases, the probability that the segment connecting two random points from a hypersphere passes near the origin vanishes.

• Carlos: I wasn’t the one who posted that, it was Bob who did the case study, so no I’m not sure the bias exists, but I can say that it doesn’t have to pass near the origin to go into useless territory, it just needs to deviate off the razor thin curved shell. Basically, if you take any two points on an eggshell and you draw a line between them what’s the probability that a randomly chosen point on the line will be inside the material of the calcified eggshell? If the points are very close together then it works ok. But if they’re far apart essentially none of the points are inside the eggshell material. If you can take a sequence of points that are close together and continue to travel in a consistent direction, you get a trajectory along the eggshell, if you are constantly changing directions, you get a random walk on the eggshell, and you never go very far.

In practice, fitting models in Stan, my complicated models often need very small timesteps and this is because they have to follow a complicated curved “surface” in N dimensions, a directed trajectory “along the eggshell”

• Carlos Ungil says:

The “algorithm” that was shown to fail (taking the median from two random points from the distribution) was a straw man. It’s quite obvious that it’s going to result in a more concentrated distribution and there is no need to go beyond one or two dimensions to realize that. Now we can discuss how in high dimensions the problem gets worse, but that’s completely irrelevant because the algorithm was never supossed to work.

The funny thing is that if you have a multi-dimensional standard normal distribution the median of two random points is distributed as a multi-dimensional normal distribution of variance 1/2. Therefore, it’s enough to modify the algorithm slightly from (a+b)/2 to (a+b)/sqrt(2) to fix the problem (whatever the number of dimensions).

Maybe ensemble algorithms have issues, but I don’t think that blog post provided any indication of it.

• Carlos Ungil says:

Your blog post seems to work equally well if you change

“most of the draws from the proposal are away from the typical set of the target!”

to

“most of the draws from the proposal are away from the high probability region of the target!”

• Aki Vehtari says:

No. We can also have a such proposal distribution (e.g. Student-t with nu=7, variance=1) that we have bounded ratios and “most of the draws from the proposal are *towards* to the high probability region of the target” and still have “most of the draws from the proposal are away from the typical set of the target”. This is illustrated in the one of the plots I left out, but I can add it later today.

• Carlos Ungil says:

Of course my comment was about what you wrote, not about other things you could have written.

For the case considered in your PS the problem is also easy to see in terms of high density regions. The HDR for the proposal covers only a small part of the HDR for the target so it is not surprising that importance sampling behaves badly.

• Aki Vehtari says:

> For the case considered in your PS the problem is also easy to see in terms of high density regions. The HDR for the proposal covers only a small part of the HDR for the target so it is not surprising that importance sampling behaves badly.

But this holds only for the second example, but not in the first example.

3. The point I was trying to make isn’t that we need to exclude the mode in some way. It’s just that a lot of people develop intuitions for algorithms based on one or two dimensional examples where there’s a lot of volume around the mode. They think the point is to try to find a mode and sample near there. That won’t actually work, because if you take draws near the mode, estimates like posterior variances are going to be wrong. The point is to take draws that let you estimate the expectations you care about.

Of course, the sense in which the mode is excluded from the typical set is also relevant—its log density is atypically high. If you take lots of draws from the posterior, their log density will be bounded pretty far away from the mode’s log density. Again, nothing needs to exclude the mode—draws near it in log density just don’t come up.

• ojm says:

What if your target functional is the maximum log probability? Will sampling from the typical set be good for this?

• Chris Wilson says:

ojm, next time you run a Stan model, you might have a look at the pairs() plot. You can examine lp_, etc. no problem, if that is what you are after here.

• Sampling won’t get you max log prob, because with sampling you will never see the mode, it’s just a fact about probability. You need optimization for this.

• Chris Wilson says:

Daniel, you will see the mode of the posterior for the marginal parameter distributions though…

• If you want a max posterior mode you need to use an optimizer. It’s why Stan has both sampling and optimization built in.

Bayesians don’t do optimization when it can be avoided. What we care about for Bayesian inference is posterior expectations. For example, point estimates are typically calculated as posterior means, because that minimizes expected squared error (relative to the model, of course); poserior median minimizes expected absolute error. The posterior mode doesn’t have any such properties. So we don’t pay much attention to the mode.

Having said that, Andrew likes to think of a Laplace approximation on the unconstrained scale as a kind of approximate Bayesian posterior. That makes it very much like variational inference, which drops an approximate distribution on an approximate posterior mean rather than posterior mode.

• ojm says:

I would think my comments are impossible to understand if it wasn’t for Carlos…

• I think it comes down to this: if you want to approximate an independent identically distributed sample from a high dimensional posterior distribution, the best available algorithm is Stan’s version of HMC. That algorithm samples in the region where you’d find your samples if they were IID and that region is roughly the region where log(p(x)) is near-constant and equal to some value that is much lower than the value at the mode.

so, either IID sampling, or HMC sampling will both give you points that are in the typical set and *not* particularly near the mode in euclidean distance.

If you want something that can’t be calculated effectively by IID sampling, such as the location of the mode, then you also won’t be able to calculate it effectively with HMC type sampling.

The point of HMC sampling is that it gives you a sample which is “close to IID” in the sense of effective sample size is nearly the same as the actual sample size. Anything you can calculate with IID samples you can calculate pretty well with HMC samples…

• Chris Wilson says:

ojm, what Daniel and Bob said is definitely technically correct, but it feels like you are trying to ask something else. My comment was only tangentially related in that I am pointing out for any marginal of interest, sampling from the typical set definitely does enable one to visually find the mode…


text:
Which multiphoton states are related via linear optics?

We investigate which pure states of $n$ photons in $d$ modes can be transformed into each other via linear optics, without post-selection. In other words, we study the local unitary (LU) equivalence classes of symmetric many-qudit states. Writing our state as $f^\dagger|\Omega\rangle$, with $f^\dagger$ a homogeneous polynomial in the mode creation operators, we propose two sets of LU-invariants: (a) spectral invariants, which are the eigenvalues of the operator $ff^\dagger$, and (b) moments, each given by the norm of the symmetric component of a tensor power of the initial state, which can be computed as vacuum expectation values of $f^k(f^\dagger)^k$. We provide scheme for experimental measurement of the later, as related to the post-selection probability of creating state $f^{\dagger k}|\Omega\rangle$ from $k$ copies of $f^{\dagger}|\Omega\rangle$.
Submitted 12 Mar 2014 to Quantum Physics [quant-ph]
Published 13 Mar 2014
Updated 14 Apr 2014
Author comments: 14 pages, 3 figures; minor changes (3 new citations, a remark on fermionic systems)
http://arxiv.org/abs/1403.3069
http://arxiv.org/pdf/1403.3069.pdf

format: plain


text:
# beamer : embedded video (.avi) in presentation

It's the first time that I'm using beamer to prepare my defense presentation. This is the beamer template that I'm using : template (pdflatex version). In fact, I need to include several videos (.avi) in my presentation. I've already tried the different ways proposed in this question. None of them worked :

\begin{frame}{test.video}
\includemedia[activate=pageopen,
width=200pt,height=150pt]{}{resources/video.avi}
\end{frame}

or

\begin{frame}{test.video}
\movie[width=8cm,height=4.5cm]{play the video}{video.avi}
\end{frame}

or

\begin{frame}{test.video}
\includemedia[
width=0.4\linewidth,
totalheight=0.225\linewidth,
activate=pageopen,
flashvars={
%important: same path as in addresource'
source=video.avi
}
]{\fbox{Click!}}{VPlayer.swf}
\end{frame}

\begin{frame}{test.video}
\includemedia[
activate=pageopen,
width=200pt,height=150pt,
flashvars={%
source=video.avi% same path as in addresource!
%   &autoPlay=true%    % optional configuration
%   &loop=true%        % variables
}
]{}{VPlayer.swf}
\end{frame}

and the last one that I tried and was unsuccessful :

\begin{frame}{test.video}
\includemovie{.85\textheight}{.85\textheight}{video.mp4}%
\end{frame}

They all create a white empty box that does not do anything. I have converted the videos to .mp4, this did not help either. I have installed flash-player for adobe (version 2019.010.20091). My operating system is : macOS High Sierra V. 10.13.6. Your help would be greatly appreciated.

• Which PDF viewer do you use? Do the video examples in the media9 manual play? – AlexG Feb 20 '19 at 11:25
• I use Adobe Acrobat Reader DC (version 2019.010.20091). No, they do not work unfortunately. – Basilique Feb 20 '19 at 11:28
• If you click the link, the media9 manual most likely opens in the Web browser, not in AR DC. Download the PDF and open it in AR DC directly. Then, the video examples should play, provided that Flash Player was properly installed. – AlexG Feb 20 '19 at 11:31
• Yes, they work perfectly when I open the file using AR DC. – Basilique Feb 20 '19 at 11:41
• Ok. Now, convert your avi to MP4. Take the code from the first box in tex.stackexchange.com/a/345433/36296 and paste it into an empty *.tex file. In the code, replace the m4v file with your MP4 file, compile and open PDF in AR DC. – AlexG Feb 20 '19 at 11:54

You mentioned in a comment that you generated the video from a series of images. Instead of trying to include the video, I would simply include the images directly.

A short example:

\documentclass{beamer}
\usepackage{xmpmulti}
\begin{document}

%asuming you images are called "something-0.png" up to "something-16.png"
\begin{frame}
\transduration<0-16>{0}
\multiinclude[<+->][format=png, graphics={width=\textwidth}]{something}
\end{frame}

\end{document}
`

• I tried but the pdf do not act as a continuous image. Your image is a gif. How can I get a gif from the generated pdf? – juan Jul 12 at 10:05


text:
# Work done in non-uniform circular motion

In non-uniform circular motion, the work done by force should not be zero because tangential force acts on body. Is this correct? If I am wrong, please give a proper explanation.

Of course. . . Work done is stored in the body as its rotational kinetic energy i.e, $$\tau \theta = \Delta KE_{rot}$$ If torque is zero, then $\tau=0$ and thus there is no change in rotational kinetic energy. Then the motion is said to be uniform circular motion.


text:
# Can PA prove very fast growing functions to be total?

The Goodstein-sequence is a total function, but PA cannot prove this.

Is this true for any other function with growth rate at least $f_{\epsilon_0}$ or are there functions growing at least as fast as the Goodstein-sequence that PA can prove to be total ?

I heard that the "power" of the PA is below the $f_{\epsilon_0}$-level, but I do not know if this answers my question.

• This is the first question I've seen correctly tagged with [proof-theory] in a very long time. It makes me sad to realize how often I have to remove this tag, and even sadder that I'm happy not to remove it here. – Asaf Karagila Dec 25 '14 at 21:40
• Thanks for this compliment. – Peter Dec 25 '14 at 21:41
• You can find some discussion in the final half of this answer and its links. – Bill Dubuque Dec 25 '14 at 23:53

Every function which eventually outgrows $f_{\varepsilon_0}$ cannot be proven to be total in Peano arithmetic. This is implied by more general result:
PA can prove a computable function $F$ to be total if and only if $F$ is primitively recursive in $f_{\omega\uparrow\uparrow n}$ for some finite $n$, where $\omega\uparrow\uparrow n=\omega^{...^\omega}$ with $\omega$ $n$ times.
In particular, if $F$ outgrows all of $f_{\omega\uparrow\uparrow n}$ for finite $n$, then PA cannot prove $F$ total.
• The result I mentioned is related to the fact that $\varepsilon_0$ is so called proof theoretic ordinal of PA. For similar result for ZFC we would need ordinal which is at least its proof theoretic ordinal, which is far beyond the reach of all known notations. Even then, we are not guaranteed such classification exists. – Wojowu Dec 25 '14 at 21:58
• @Peter ZFC can actually handle fast growing hierarchy up to $\Gamma_0$ (a.k.a. Feferman-Schutte ordinal). – Wojowu Dec 25 '14 at 22:02
• @Peter In fact, the proof theoretic ordinal of ZFC is much, much larger than $\Gamma_0$. Even the Takeuti-Feferman-Buchholz ordinal $\psi(\varepsilon_{\Omega_\omega+1})$ is not coming close. Nor is the proof theoretic ordinal of KPI, which is even larger. The comparison of these ordinals with the proof theoretic ordinal of ZFC is like the comparison of a googolplex with Graham's number, to put it mildly. – wythagoras Jul 27 '15 at 8:14


text:
# Math Help - help with very small equation

1. ## help with very small equation

i have
.50x+80=.-80x=68

do i need to subract the negative .8x from the .5x or the other way around

or do i need to subract 68 from 80 and get 12
but after i geat 12 i will need to divide by -.30x and my answer will be negative for my problem and it wont make any sense
apreciate any help

2. can you rewrite the equation a bit neater please.

because what im getting from it is:

$

0.50x + 80 = -0.80x = 68
$

... :/

3. ## sorry

i meant to say
.5x + 80 +.-80x =68
what are steps i need to do?

4. Equation:

$
0.5x + 80 - 0.8x = 68
$

Combine Like Terms:

$
80 - 0.3x = 68
$

Move The Constant To The Other Side:

$
-0.3x = -12
$

Divide Through By -0.3:

$
x = 40
$

Understand?


text:
## Intermediate Algebra (12th Edition)

Published by Pearson

# Chapter R - Section R.1 - Basic Concepts - R.1 Exercises: 10

#### Answer

$k=\{ ...,-7,-5,-3,-1 \}$

#### Work Step by Step

Integers are the negative, zero, and positive numbers which include $\{...,-3,-2,-1,0,1,2,3...\}$. Odd integers are integers not divisible by $2$. Hence, the odd integers, $k$, less than $1$ are \begin{array}{l}\require{cancel} k=\{ ...,-7,-5,-3,-1 \} .\end{array}

After you claim an answer you’ll have 24 hours to send in a draft. An editor will review the submission and either publish your submission or provide feedback.


text:
# Letters of recommendation

Sometimes high school students ask me for reference letters. (Some of them even got in!) So here are some details. Inspired by Ravi Vakil.

## Programs I can write for#

I write letters most commonly for RSI and typical undergraduate colleges (MIT and Common App). For the latter, most people prefer to invite me as a supplemental letter rather than as a math teacher proper. I can also write letters for MIT PRIMES, although I am of the impression that the awards and entrance exam hold much more weight than rec letters for PRIMES (compared to, say, college apps).

I’ve also written letters for summer math camps such as Canada/USA MathCamp, PROMYS/ROSS, etc.

Letters from me will necessarily be a bit unusual, so you should think about whether I’m the right person to do this. Here are a few reasons why:

• I am not formally a teacher or professor, which may confuse any clueless admissions officers1.
• Also, if I only worked with you at MOP or OTIS, I may not have spent enough time with you in-person to “paint a picture” compared to a teacher at your school, so to speak.
• You are going to be compared to other OTIS/MOP students. This is not to say the comparison will necessarily be unfavorable. However, the letter is probably not going to say you are the strongest student I’ve ever met, either, whereas many of you are the best student your local high school has ever seen.

Usually people will ask me for letters if they want to add some context about their math contest experiences that’s more than just score, or to have me advocate for their mathematical ability. I try to also, when I have enough data, say a bit about what you’re like as a person when interacting with me or other OTIS students.

Since I’m a nice person, I promise to not intentionally write a negative letter.2 If I feel I cannot write you a letter with nonnegative value, I will decline your request instead.3

## Sample letter#

Here is a sample letter for an imaginary student named Eva Chan, so you know roughly what you’re signing up for. (And those of you writing your first reference letter are welcome to use this as an example.)

## Instructions#

If you do ask a letter for me, I need the following information.

1. Please remind me which years I worked with you (in OTIS/MOP/SPARC/etc.). I have been doing this for so long that the years are all mixed up in my head, so I’d like to double-check!
2. Provide any recent math olympiad results, like USAMO / TSTST / TST / IMO. (I may not use them, but I want to have them if I need them).
3. Is there anything in particular I should discuss? In particular, how math-focused do you want me to be?
4. Please send me any application essays or completed forms you are willing to share (drafts of essays are okay). This way, I can tailor my letter to complement what you’ve already written.
5. If I’ve written a letter for you in the past, please remind me about it.
6. Please also specify which email address you applied to the program with.
7. Specify what “type” of letter I am writing, if applicable. For example, colleges commonly distinguish between, say “required” letters from a math teacher, or optional “supplements” (see e.g. MIT, second-to-last paragraph). It is useful for me to know whether I am the main voice or if I am expecting to add on to what other (usually more traditional) teachers have to say.
8. State explicitly when the deadline is.

• Ideally, ask at least a month in advance. (If it’s less than one month, you can still ask, but I may decline.)
• Requests, and accompanying materials/info, must be sent by email. Not Discord direct message. I can’t keep track of those.
• When inviting me to whatever application the school or program uses, it’s best if you invite the address evan at evanchen.cc.
• Exception: for some reason, MathPrograms.org (which is used by REU’s and MIT PRIMES) does not seem to work well with my professional email. You should use evanchen at mit.edu there.
• Usually, I will always submit letters at least a week in advance if you asked ahead of time. So if it’s less than a week until the deadline and the letter I promised hasn’t been submitted yet, then please double-check on it.

1. Although I am gradually picking up increasingly fancy titles.

2. Regrettably, I speak from experience here. When applying to graduate school and fellowships, a sympathetic person tipped me off that one of my letters was “not doing me any favors” and warned me to delete it from my graduate school applications if possible. I am still grateful to him today since I would likely not have been admitted to MIT’s math PhD program otherwise. Alas, this tip was too late for all the other schools, as I had already finalized my application at that point.

3. This is because I believe a good mentor should be unconditionally supportive, much like a parent, and that it’s not my place to judge students.

Updated Tue 27 Dec 2022, 19:05:13 UTC by eeeb3f2b8b87


text:
# how did mr. shimerda death have an effect on the usual routine of the courtry?

book 1 chapters 15-19


text:
Question

# The refractive index of diamond is 2.42 and the speed of light in air is $3 \times {10^8}m{s^{ - 1}}$ calculate the speed of light in diamond.

Verified
129.3k+ views
Hint: One of the most important properties of the wave is refraction. In this phenomenon, the rays of light get bent when they pass through different mediums and when we provide the ratio of the speed of light (wave) in various media, it is called a refractive index.

Complete step by step solution:
The refractive index $\left( n \right)$ is given by,
$n = \dfrac{{Speed\;of\;light\,\;in\;vacuum}}{{Speed\,\,of\;light\;in\,\;medium\,}}$
Refractive index of diamond = 2.42
Speed of light = $3 \times {10^8}m{s^{ - 1}}$
Speed of the light in diamond = ?
As we all know that,
$n = \dfrac{{Speed\,of\,light\,in\,vacuum}}{{Speed\,\,of\,Light\,in\,medium\,}}$
$n = \dfrac{c}{V}$
$2.42 = \dfrac{{3 \times {{10}^8}}}{{speed\,\,of\,\,light\;in\;medium}}$
$Speed\,\,of\,light\;in\;medium\, = \,\dfrac{{3 \times {{10}^8}}}{{2.42}} \\ = 1.25\, \times {10^8}\,m/s \\$

Note: Speed of light in vacuum is the fastest speed in the universe that is $3 \times {10^8}\;m/s$
Refractive index is the ratio of the velocity of light in a vacuum to velocity of light in media.


text:
## FANDOM

10,835 Pages

There may be some informal or invalid step in my proof. Just point it out. There's still a lot for me to learn.

Feeding FGH into itself was considered before. I just want to find the limit of it and probably some usage of it.

It looks like this:

1. $f_{f_\alpha(n)}(n)$ , or

2. $f_{f_\alpha(\omega)}(n)$

1. and 2. is very different. 1. comes to an limit very early, while 2. can go very far.

## The first situation

the limit of $f_{f_\alpha(n)}(n)$ is $f_{f_{f_{...}(n)}(n)}(n)$, while n is an integer.

when m,n are integers, for integer n large enough , $2\uparrow^{m}n>f_m(n)>2\uparrow^{m-1}n$.

Therefore, $2\uparrow^{2\uparrow^{2\uparrow^{...}n}n}n>f_{f_{f_{...}(n)}(n)}(n)>2\uparrow^{2\uparrow^{2\uparrow^{...}n-1}n-1}n$

Keeping the integer n big enough ,then $f_{\omega+1}(n)>{\{n,n,1,2}\}=n\uparrow^{n\uparrow^{n\uparrow^{...}n}n}n>2\uparrow^{2\uparrow^{2\uparrow^{...}n}n}n$

Here we know that $f_{\omega+1}(n)>f_{f_{f_{...}(n)}(n)}(n)>2\uparrow^nn>f_\omega(n)$

So, $f_{\omega+1}(n)>f_{\alpha\mapsto{f_\alpha(n)}}(n)>f_\omega(n)$

## The second situation

The second situation is more fierce and powerful. Since it involves FGH with n replaced by ordinals. It will be harder to define.

Here, I will not define the FGH with transfinite ordinals, and will not prove the limit of it. I'll just write down my idea about it and the possible limit of it. I will keep studying about further definition and proof in the future.

Slow Growing Hierarchy, another hierarchy that grows more slowly, has very unique properties.

for specific function $f(n)$, it looks like that:

$g_{f(\omega)}(n)\approx f(n)$

$g_{f(\omega)}(\omega)\approx f(\omega)$

So,

$\alpha=f_\alpha(\omega)=f_{f_\alpha(\omega)}(\omega)\approx f_{g_\alpha(\omega)}(\omega)$

is actually where the $f_\alpha(\omega)\approx g_\alpha(\omega)$, the point that SGH(n) catches up FGH(n).

There are various results of when it occurs, but most commonly it's $\psi_0(\Omega_\omega)$.

Thus, $f_{\alpha\mapsto{f_\alpha(\omega)}}(n)\approx f_{\psi_0(\Omega_\omega)}(n)$. It's theoretically the limit of FGH itself. But with Ordinal Collapsing Function that is more powerful and inaccessible ordinals, it is possible for FGH to go beyond its limit.

--D57799 (talk) 06:24, October 5, 2014 (UTC)


text:
Publication 2011 Issue No. 1 - January-February Abstract - A Distributed Algorithm for Finding All Best Swap Edges of a Minimum-Diameter Spanning Tree
 This Article Share Bibliographic References Add to: Digg Furl Spurl Blink Simpy Google Del.icio.us Y!MyWeb Search Similar Articles Articles by Beat Gfeller Articles by Nicola Santoro Articles by Peter Widmayer
A Distributed Algorithm for Finding All Best Swap Edges of a Minimum-Diameter Spanning Tree
January-February 2011 (vol. 8 no. 1)
pp. 1-12
 ASCII Text x Beat Gfeller, Nicola Santoro, Peter Widmayer, "A Distributed Algorithm for Finding All Best Swap Edges of a Minimum-Diameter Spanning Tree," IEEE Transactions on Dependable and Secure Computing, vol. 8, no. 1, pp. 1-12, January-February, 2011.
 BibTex x @article{ 10.1109/TDSC.2009.17,author = {Beat Gfeller and Nicola Santoro and Peter Widmayer},title = {A Distributed Algorithm for Finding All Best Swap Edges of a Minimum-Diameter Spanning Tree},journal ={IEEE Transactions on Dependable and Secure Computing},volume = {8},number = {1},issn = {1545-5971},year = {2011},pages = {1-12},doi = {http://doi.ieeecomputersociety.org/10.1109/TDSC.2009.17},publisher = {IEEE Computer Society},address = {Los Alamitos, CA, USA},}
 RefWorks Procite/RefMan/Endnote x TY - JOURJO - IEEE Transactions on Dependable and Secure ComputingTI - A Distributed Algorithm for Finding All Best Swap Edges of a Minimum-Diameter Spanning TreeIS - 1SN - 1545-5971SP1EP12EPD - 1-12A1 - Beat Gfeller, A1 - Nicola Santoro, A1 - Peter Widmayer, PY - 2011KW - Fault-tolerant routingKW - swap edgesKW - minimum-diameter spanning treeKW - distributed algorithms.VL - 8JA - IEEE Transactions on Dependable and Secure ComputingER -
Beat Gfeller, ETH Zurich, Zurich
Nicola Santoro, Carleton University, Ottawa
Peter Widmayer, ETH Zurich, Zurich
Communication in networks suffers if a link fails. When the links are edges of a tree that has been chosen from an underlying graph of all possible links, a broken link even disconnects the network. Most often, the link is restored rapidly. A good policy to deal with this sort of transient link failures is swap rerouting, where the temporarily broken link is replaced by a single swap link from the underlying graph. A rapid replacement of a broken link by a swap link is only possible if all swap links have been precomputed. The selection of high-quality swap links is essential; it must follow the same objective as the originally chosen communication subnetwork. We are interested in a minimum-diameter tree in a graph with edge weights (so as to minimize the maximum travel time of messages). Hence, each swap link must minimize (among all possible swaps) the diameter of the tree that results from swapping. We propose a distributed algorithm that efficiently computes all of these swap links, and we explain how to route messages across swap edges with a compact routing scheme. Finally, we consider the computation of swap edges in an arbitrary spanning tree, where swap edges are chosen to minimize the time required to adapt routing in case of a failure, and give efficient distributed algorithms for two variants of this problem.

[1] B. Gfeller, N. Santoro, and P. Widmayer, "A Distributed Algorithm for Finding All Best Swap Edges of a Minimum Diameter Spanning Tree," Proc. 21st Int'l Symp. Distributed Computing (DISC), pp. 268-282, 2007.
[2] M. Bui, F. Butelle, and C. Lavault, "A Distributed Algorithm for Constructing a Minimum Diameter Spanning Tree," J. Parallel and Distributed Computing, vol. 64, pp. 571-577, 2004.
[3] H. Ito, K. Iwama, Y. Okabe, and T. Yoshihiro, "Single Backup Table Schemes for Shortest-Path Routing," Theoretical Computer Science, vol. 333, no. 3, pp. 347-353, 2005.
[4] E. Nardelli, G. Proietti, and P. Widmayer, "Finding All the Best Swaps of a Minimum Diameter Spanning Tree under Transient Edge Failures," J. Graph Algorithms and Applications, vol. 5, no. 5, pp. 39-57, 2001.
[5] E. Nardelli, G. Proietti, and P. Widmayer, "Swapping a Failing Edge of a Single Source Shortest Paths Tree Is Good and Fast," Algorithmica, vol. 35, no. 1, pp. 56-74, 2003.
[6] A.D. Salvo and G. Proietti, "Swapping a Failing Edge of a Shortest Paths Tree by Minimizing the Average Stretch Factor," Theoretical Computer Science, vol. 383, no. 1, pp. 23-33, 2007.
[7] P. Flocchini, A. Mesa Enriques, L. Pagli, G. Prencipe, and N. Santoro, "Point-of-Failure Shortest-Path Rerouting: Computing the Optimal Swap Edges Distributively," IEICE Trans. Information and Systems, vol. E89-D, no. 2, pp. 700-708, 2006.
[8] P. Flocchini, L. Pagli, G. Prencipe, N. Santoro, and P. Widmayer, "Computing All the Best Swap Edges Distributively," J. Parallel and Distributed Computing, vol. 68, no. 7, pp. 976-983, 2008.
[9] F.T. Leighton, Introduction to Parallel Algorithms and Architectures: Arrays $\cdot$ Trees $\cdot$ Hypercubes. Morgan Kaufmann Publishers, Inc., 1992.
[10] N. Santoro, Design and Analysis of Distributed Algorithms, Wiley, 2007.
[11] N. Santoro and R. Khatib, "Labelling and Implicit Routing in Networks," The Computer Journal, vol. 28, no. 1, pp. 5-8, 1985.
[12] M. Thorup and U. Zwick, "Compact Routing Schemes," Proc. 13th Ann. ACM Symp. Parallel Algorithms and Architectures (SPAA '01), pp. 1-10, 2001.

Index Terms:
Fault-tolerant routing, swap edges, minimum-diameter spanning tree, distributed algorithms.
Citation:
Beat Gfeller, Nicola Santoro, Peter Widmayer, "A Distributed Algorithm for Finding All Best Swap Edges of a Minimum-Diameter Spanning Tree," IEEE Transactions on Dependable and Secure Computing, vol. 8, no. 1, pp. 1-12, Jan.-Feb. 2011, doi:10.1109/TDSC.2009.17


text:
0
Chapter 23
Solidification in a Casting Mold

## Excerpt

Exact temperature distribution solutions to unsteady-state one-dimensional heat conduction equations in rectangular, cylindrical, and spherical coordinates, for initial condition and different boundary conditions, have been provided in literature (see References [1], [6], [10] and [16]). The reference by Carslaw, H. S. and J. C. Jaeger [1] extended these exact solutions to moving boundary conditions with phase change.

In this chapter, sensitivity analysis for the solidification front of a semi-infinite liquid in a semi-infinite mold is considered. The unsteady-state one-dimensional heat conduction equation in rectangular coordinates for the mold region, solidified cast material region and liquid cast material region, with constant property assumptions, is given by:

$∂2T∕∂x2=(1∕α)∂T∕∂θ$

\$25.00

### Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related Proceedings Articles
Related eBook Content
Topic Collections


text:
# Theory of proveit.physics.quantum.QPE¶

In this package, we prove outcome bounds for the Quantum Phase Estimation algorithm. Our proof is based on the proof in pp. 221--225 of Michael A. Nielsen and Isaac L. Chuang, "Quantum Computation and Quantum Information", (10th Anniversary Edition). Cambridge University Press, New York, NY, 2010.

In [1]:
import proveit
%theory # toggles between interactive and static modes


text:
ghc-check-0.1.0.2: detect mismatches between compile-time and run-time versions of the ghc api

Returns True if the run-time version of the ghc package matches the compile-time version
Returns the compile-time version of the ghc package
Returns the run-time version of the ghc package by looking up in the package database


text:
## Elementary and Intermediate Algebra: Concepts & Applications (6th Edition)

Published by Pearson

# Chapter 10 - Exponents and Radicals - Review Exercises: Chapter 10 - Page 694: 47

#### Answer

$5\sqrt{2} \text{ } cm$

#### Work Step by Step

The diagonal and 2 sides of the square form a right triangle. Let $c$ be the hypotenuse and $a,b$ be the equal legs. Using $a^2+b^2=c^2$ or the Pythagorean Theorem, where $c=10$ and $a=b,$ then \begin{array}{l}\require{cancel} a^2+b^2=c^2 \\\\ a^2+a^2=10^2 \\\\ 2a^2=100 \\\\ a^2=\dfrac{100}{2} \\\\ a^2=50 \\\\ a=\sqrt{50} \\\\ a=\sqrt{25\cdot2} \\\\ a=\sqrt{(5)^2\cdot2} \\\\ a=5\sqrt{2} .\end{array} Hence, the side of the square (which has the same value as the legs) is equal to $5\sqrt{2} \text{ } cm .$

After you claim an answer you’ll have 24 hours to send in a draft. An editor will review the submission and either publish your submission or provide feedback.


text:
# HP Forums

Full Version: SVD
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
I can't seem to get SVD to work for non-square matrices. This works perfectly on the 50g but fails on the Prime...unless I'm doing something wrong. Anyone else seen this?
On the latest firmware:

$SVD\left ( \begin{bmatrix} 1 & -1 \\ 1 & -1 \\ \sqrt{3} & 0 \end{bmatrix} \right )$

Gives "Error: Bad Argument Type" in both Home and CAS.

Looks like a definite bug to me.

Mark Hardman
If it is a bug, it is a minor bug.

Try entering the SVD command in lowercase:

$\textrm{simplify} \left ( \textrm{svd}\left ( \begin{bmatrix} 1 & -1 \\ 1 & -1 \\ \sqrt{3} & 0 \end{bmatrix} \right ) \right )$

Gives a result of:

$\begin{bmatrix} 0.5477 & 0.7071 & -0.4472 \\ 0.5477 & -0.7071 & -0.4472 \\ 0.6325 & 5.1840 \times 10^{-17} & 0.7746 \end{bmatrix}$

Mark Hardman
Mark,

Thanks. How did you think to do this? The toolbox menu inserts in uppercase and the manual defines SVD in uppercase. Did you just try at random or is there some logic to when this matters?

Also, using SVD() (uppercase) on a square matrix seems to work...very odd!
Search the old forum archives. The difference between UPPERCASE/lowercase for commands on the Prime has already generated several other bug reports. It is a matter of confusion that the Prime Team will hopefully address in a future release.

Mark Hardman
Thanks, I'll check it out. I admit I want to like this calculator but it's been a bit of a trial so far!
The current algorithm covered only square real matrices (as you can guess, implementing numeric algorithms is not my main business, and giac/xcas uses LAPACK for svd, therefore it was only required here), svd was also sometimes able to return an answer for rectangular matrices. I'm improving it, it will handle rectangular real or complex matrices that are well conditionned (0 or small singular values will not be handled, I'm not writing a dedicated numerical algorithm for that, just reusing existing routines).
Parisse,

Thanks for the clarification. Nice to know these things are being worked out. As long as there's a reasonable error message for ill-conditioned matrices I'll be happy

It would also be nice if it worked correctly in RPN: svd() should take a matrix from the stack and just push U, S, V back onto the stack.

--jeff
Actually, you don't really need a message for ill-conditionned matrix, just look at the ratio between the largest singular value and the smallest one in s.
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :


text:
# Tag Info

## Hot answers tagged deep-rl

Accepted

### Mathematically, what is happening differently in the neural net during exploration vs. exploitation?

Typically, the NN is trained the same way whether an action is chosen for exploration or exploitation. Look at the objective (AKA loss) function for any algorithm you're interested in and you'll ...
• 451
Accepted

### Why do Q-values diverge without a target network?

I don't see how the target Q-value gets updated when the current Q-value is changed. Without a separate target network, this happens because the approximator will generalise, and the generalisation ...
• 24.5k
Accepted

### How do I design the network for Deep Q-Network?

What is the strategy to get to a better network? There are a few different strategies that you can use to search for good hyperparameters in reinforcement learning RL, but you should be aware that ...
• 24.5k

### What is the difference between a loss function and reward/penalty in Deep Reinforcement Learning?

1. Question: The difference between loss and reward/penalty So I see both the loss function and the reward/penalty are the quantitative way of measuring the output/action and making the model to ...
• 671

### What is the difference between a loss function and reward/penalty in Deep Reinforcement Learning?

Ultimately, in RL, the policy is what you want to find. It's the solution to the Markov Decision Process (MDP). But you don't want to find any policy, but the optimal policy, i.e. the one that will ...
• 34.4k
Accepted

### What can we learn from AlphaZero in the development towards AGI?

Some learnings from AlphaZero: Self-play, and more generally sandbox training, is effective. This indicates that given the right enviroment and enough computational power we can build highly ...
• 361
Accepted

### Determining to terminate at a reward or not

I am trying to understand how you would determine whether it is better for the agent to terminate at the state with the number 3 or to continue to the state with a number 4 to collect the more reward? ...
• 24.5k

### How to deal with small reward values

The numbers that a value-based neural network will predict are usually based on expected returns (sum of rewards by end of an episode, or a discounted infinite sum), although in some cases they might ...
• 24.5k
Accepted

### What do equations 1 and 3 describe in the "Image Augmentation Is All You Need: Regularizing Deep Reinforcement Learning from Pixels" paper?

Equation 1 In normal Q-Learning your target is defined as $y_t = r_t + \gamma \mathrm{max_a}Q(s_{t+1}, a)$. Since you're training a regularized version, you construct the estimated value of the next ...
• 348

### Should I apply normalization to the observations in deep reinforcement learning?

The use of normalisation in neural networks and many other (but not all - decision trees are a notable exception) machine learning methods, is to improve the quality of the parameter space with ...
• 24.5k
Accepted

• 136
Accepted

### How should I initialize the weights of the neural network so that the initial policy is uniform?

You may be interested in section 3.2 of this paper What Matters In On-Policy Reinforcement Learning? A Large-Scale Empirical Study (2020) by Google Research. They claim that the initialization of the ...
• 931

### How Come My (D)DQN Fails To Learn?

I have two suggestions that you can look into. Based on my own work in RL, I believe the first one will require less work to implement. If the observability of the environment is not an issue, then ...
• 179

### Do I need to normalize all state-space variables? If so, how?

It's likely to train as long as they're reasonably within the orders of magnitude of other normalized variables. The network can adjust for that. But it might cause problems later, if the values move ...
• 451
Accepted

### Do I need to normalize all state-space variables? If so, how?

The way I've seen most codes treat the state normalization is that they simply take a running mean and standard deviation for each dimension of the state space. As you point out, this normalization ...
• 931

### Do I need to normalize all state-space variables? If so, how?

I'll start with the literal question in the title: Do I need to normalize all state-space variables? You don't strictly need to in theory. It's often really useful, or sometimes borderline necessary,...
• 9,419
Accepted

### Why do we use "true labels" that are based on the output of our network in Deep Q-Learning?

The labels in DQN, and in Q-learning in general, are not "true" in the sense that they represent optimal action value functions. Instead they represent approximate action values of a current ...
• 24.5k
Accepted

### How to perform the back propagation step in Semi-gradient SARSA using a deep neural network?

If so, then would it be correct to take the gradient of this expression of $\hat{q}(S,A,w)$ with respect to the weights of the network Yes, your expression for $\hat{q}(S,A,w)$ looks correct for your ...
• 24.5k
Accepted

### In RL, is it possible to design a multiplicative/exponential reward function? A reward func that depends on current accumulated reward?

The main thing you will need to do is add the accumulated reward (total_score_so_far) to the state. In order to predict future reward with any accuracy, the agent ...
• 24.5k

### Why is a large replay buffer inefficient?

I read the same thing recently, and my interpretation was this: If you only use the very-most recent data, you will overfit to that and things will break We'd like to train the network to predict ...
• 451
Accepted

### How do I create an AI controller for Pacman?

About the environments For the controller part of your question, I would advice looking at openAI gym. https://www.gymlibrary.ml/content/environment_creation/ #how to make your own gym enviroment ...
• 349
Accepted

### How to sample the tuples during the initial time steps of the DDPG algorithm?

You have a free choice to either: Wait until the replay buffer hits a minimum size for sampling. Take smaller samples from the buffer initially, until the buffer is large enough. On the first time ...
• 24.5k
Accepted

### Which paper describes the effect of learning_starts in Reinforcement Learning?

The replay buffer allows breaking the temporal dependence of the data and thus makes them more i.i.d. (which is what we want). The replay buffer needs to be at least filled with enough experience to ...
1 vote

### How to take gradient of log policy when actions are negative?

The policy gradient tells us that $$\nabla_\theta v_\pi(s) = \mathbb{E}_\pi\left[ G_t \nabla_\theta \log \pi_\theta(a | s) \right] \; ;$$ where as usual $G_t$ are the returns, $\pi$ is policy ...
• 4,180
1 vote
Accepted

### How to define actions on a list of values?

I had tried working on a problem similar to this using combinatorial scoring games. I ran into other issues with the players competing, but I think I can give some advice to how I handled this. In my ...
• 46
1 vote
Accepted

### Why does my actor-critic network always give either -1 or 1 at the output layer?

Most probably your network is underfitted. In that case, the network outputs values randomly. Hyperbolic tangent tanh converges very quickly towards $-1$ or $1$, so that is why you always find $-1$ ...
• 522
1 vote

### Does AlphaGo play random moves in a real competition?

The core mechanics of AlphaZero during selfplay and real tournament games are the same: something similar to Monte Carlo Tree Search is done but guided by the current neural network instead of random ...
• 146

Only top scored, non community-wiki answers of a minimum length are eligible


text:
# Wrapfigure side by side on two columns

I need to insert two figures such that the figures would be side by side, starting from the left column and the second figure partially spanned on the next column. I've tried both wrapfig and minipage with no luck. Example of the requested layout:

Any suggestions would be greatly appreciated!

If you make a figure at the top of the left column and put in contents (pictures) that is too wide they will expand into the right column. Then you can use e.g. wrapfig in the right column to make space for the extra width. It will robably be a bit tricky to start the wrapfigure at the top of the right column. Here I have not put any effort at ending the left column in a nice way, I just made a \newpage when I could not fit another paragraph.

\documentclass[twocolumn]{article}
\usepackage{lipsum}
\usepackage{graphicx}
\usepackage{wrapfig}
\begin{document}
\title{My title}
\author{Me}
\maketitle
\begin{abstract}
\lipsum[1]
\end{abstract}
\section{Introduction}
\lipsum[1-4]%% On left column second page

\begin{figure}[t]%% Put a picture that is too wide. It will continue in to the right column
\begin{minipage}[t]{0.7\linewidth}
\includegraphics[width=\linewidth]{example-image}
\caption{First picture}
\label{fig:first}
\end{minipage}\rule{1em}{0pt}%
\begin{minipage}[t]{0.7\linewidth}
\includegraphics[width=\linewidth]{example-image}
\caption{Second picture}
\label{fig:second}
\end{minipage}
\end{figure}

\lipsum[1]%% Now at the bottom of left column on second page

\newpage %% On top of right column
\begin{wrapfigure}{L}{0.45\linewidth}%% Make room for the too wide figure in the left column
\rule{0pt}{15em}% height of wrapfigure
\end{wrapfigure}

\lipsum[1-5]
\end{document}

As a note, I am not sure I would recommend this. If you have two fiugures that are too wide to put in one column I would either put them in seprate figures after each other or in on that spans both columns, e.g.

\begin{figure*}[t]
\begin{minipage}[t]{0.47\linewidth}
\includegraphics[width=\linewidth]{example-image}
\caption{First picture}
\label{fig:first}
\end{minipage}\hfill%
\begin{minipage}[t]{0.47\linewidth}
\includegraphics[width=\linewidth]{example-image}
\caption{Second picture}
\label{fig:second}
\end{minipage}
\end{figure*}

• what if my images doesn't show in the start of the page? eliminating the \newpage yields a wrong layout. – Guy Cohen Sep 8 '17 at 10:36
• @GuyCohen That is why it is at the top of page. There is no automatic alignment between the figure and wrapfigure so it is easier to do that way. You can of course put it anywhere but the positioning of the wrapfigure becomes harder. – StefanH Sep 8 '17 at 11:32


text:
# Introduction

The kinship matrix is foundational for random effects models with family data.
For $$n$$ subjects it is an $$n \times n$$ matrix whose $$ij$$ element contains the expected fraction of alleles that would be identical by descent if we sampled one from subject $$i$$ and another from subject $$j$$. Note that the diagonal elements of the matrix will be 0.5 not 1: when we randomly sample twice from the same subject (with replacement) we will get two copies of the gene inherited from the father 1/4 of the time, the maternal copy twice (1/4) or one of each 1/2 the time. The formal definition is $$K(i,i) = 1/4 + 1/4 + 1/2 K(m,f)$$ where $$m$$ and $$f$$ are the father and mother of subject $$i$$.

The algorithm used is found in K Lange, , Springer 1997, page 71–72.

The key idea of the recursive algorithm for $$K(i,j)$$ is to condition on the gene selection for the first index $$i$$. Let $$m(i)$$ and $$f(i)$$ be the indices of the mother and father of subject $$i$$ and $$g$$ be the allele randomly sampled from subject $$i$$, which may of either maternal or paternal origin.

\begin{align} K(i,j) &= P(\mbox{g maternal}) * K(m(i), j) + P(\mbox{g paternal}) * K(f(i), j) \label{recur0} \\ &= 1/2 K(m(i), j) + 1/2 K(f(i), j) \label{recur1} \\ K(i,i) &= 1/2(1 + K(m(i), f(i))) \label{self} \end{align}

The key step in equation is if $$g$$ has a maternal origin, then it is a random selection from the two maternal genes, and its IBD state with respect to subject $$j$$ is that of a random selection from m(i) to a random selection from $$j$$. This is precisely the definition of $$K(m(i), j)$$. The recursion does not work for $$K(i,i)$$ in equation since once we select a maternal gene the second choice from $$j$$ cannot use a different maternal gene.

For the recurrence algorithm to work properly we need to compute the values of $$K$$ for any parent before the calculations for their children. Pedigree founders (those with no parents) are assumed to be unassociated, so for these subjects we have \begin{align*} K(i,i) &= 1/2 K(i,j) &=0 \; i\ne j \end{align*}

The final formula slightly different for the $$X$$ chromosome. Equation still holds, but for males the probability that a selected $$X$$ chromosome is maternal is 1, so when $$i$$ a male the recurrence formula becomes $$K(i,j) = K(m(i),j)$$.
For females it is unchanged. All males will have $$K(i,i) = 1$$ for the $$X$$ chromosome.

In order to have already-defined terms on the right hand side of the recurrence formula for each element, subjects need to be processed in the following order

The kindepth routine assigns a plotting depth to each subject in such a way that parents are always above children.
For each depth we need to do the compuations of formula twice. The first time it will get the relationship between each subject and prior generations correct, the second will correctly compute the values between subjects on the same level. The computations within any stage of the above list can be vectorized, but not those between stages.

Let [[indx]] be the index of the rows for the generation currently being processed, say generation $$g$$.
We add correct computations to the matrix one row at a time; all of the calculations depend only on the prior rows with the exception of the [i,i] element. This approach leads to a for loop containing operations on single rows/columns.

At one point below we use a vectorized version. It looks like the snippet below

for (g in 1:max(depth)) {
indx <- which(depth==g)
kmat[indx,] <- (kmat[mother[indx],] + kmat[father[indx], ])/2
kmat[,indx] <- (kmat[,mother[indx]] + kmat[,father[indx],])/2
for (j in indx) kmat[j,j] <- (1 + kmat[mother[j], father[j]])/2
}

The first line computes all the values for a horizontal stripe of the matrix. It will be correct for columns in generations $$<g$$, unreliable for generation $$g$$ with itself because of incomplete parental relationships, and zero for higher generations. The second line does the vertical stripe, and because of the line before it does have the data it needs and so gets all the stripe correct. Except of course for the diagonal elements, for which formula does not hold. We fill those in last. We know that vectorized calculations are always faster in R and I was excited to figure this out. The unfortunate truth is that for this code it hardly makes a difference, and for the X chromosome calculation leads to impenetrable if-then-else logic.

The program can be called with a pedigree, a pedigree list, or raw data. The first argument is [[id]] instead of the more generic [[x]] for backwards compatability with an older version of the routine. We give founders a fake parent of subject $$n+1$$ who is not related to anybody (even themself); it avoids some if-then-else constructions.

For S3 method dispatch to work on pedigree and pedigreeList objects, we need to define a kinship(), kinship.default(), kinship.pedigree(), and kinship.pedigreeList(). The default version takes the vectors

## Pedigree Object

The method for a pedigree object is an almost trivial modification. Since the mother and father are already indexed into the id list it has two lines that are different, those that create mrow and drow. The other change is that now we potentially have information available on monozygotic twins. If there are any such, then when the second twin of a pair is added to the matrix, we need to ensure that the pairs kinship coefficient is set to the self-self value. This can be done after each level is complete, but before children for that level are computed. If there are monozygotic triples, quadruplets, etc. this computation gets more involved.

The total number of monozygotic twins is always small, so it is efficient to fix up all the monzygotic twins at each generation. A variable [[havemz]] is set to TRUE if there are any, and an index array [[mzindex]] is created for matrix subscripting.

## PedigreeList Object

For the Minnesota Family Cancer Study there are 461 families and 29114 subjects. The raw kinship matrix would be 29114 by 29114 which is over 5 terabytes of memory, something that clearly will not work within S.
The solution is to store the overall matrix as a sparse Matrix object. Each family forms a single block. For this study we have [[n <- table(minnbreast$famid); sum(n*(n+1)/2)]] or 1.07 million entries; assuming that only the lower half of each matrix is stored. The actual size is actually smaller than this, since each family matrix will have zeros in it — founders for instance are not related — and those zeros are also not stored. The result of each per-family call to kinship will be a symmetric matrix. We first turn each of these into a dsCMatrix object, a sparse symmetric form. The [[bdiag]] function is then used to paste all of these individual sparse matrices into a single large matrix. Why do we note use [[(i in famlist)]] below? A numeric subscript of [[[9]]] %’ selects the ninth family, not the family labeled as 9, so a numeric family id would not act as we wished. If all of the subject ids are unique, across all families, the final matrix is labeled with the subject id, otherwise it is labeled with family/subject. ## MakeKinship The older [makekinship] function, from before the creation of pedigreeList objects, accepts the raw identifier data, along with a special family code for unrelated subjects, as produced by the [[makefamid]] function. All the unrelated subjects are put at the front of the kinship matrix in this case rather than within the family. Because unrelateds get put into a fake family, we cannot create a rational family/subject identifier; the id must be unique across families. We include a copy of the routine for backwards compatability, but do not anticipate any new usage of it. Like most routines, this starts out with a collection of error checks. ## Monozygotic Twins, for Pedigree object Return now to the question of monzygotic sets, used specifically in the kinship for pedigree objects. Consider the following rather difficult example: Subjects 1, 2, 3, and 7 form a monozygotic quadruple, 5/6 and 9/10 are monzygotic pairs. First create a vector which contains for each subject the lowest index of a monozygotic twin for that subject. For non-twins it can have any value. For this example that vector is set to 1 for subjects 1, 2, 3, and 7, to 5 for 5 and 6, and to 9 for 9 and 10. Creating this requires a short while loop. Once this is in hand we can identify the sets. Next make a matrix that has a row for every possible pair. Finally, remove the rows that are identical. The result is a set of all pairs of observations in the matrix that correspond to monozygotic pairs. # Older routines For testing purposes we have a version of the kinship function prior to addition of the chrtype argument. oldkinship <- function(id, ...) { UseMethod('oldkinship') } oldkinship.default <- function(id, dadid, momid, ...) { n <- length(id) if (n==1) return(matrix(.5,1,1, dimnames=list(id, id))) if (any(duplicated(id))) stop("All id values must be unique") kmat <- diag(n+1) /2 kmat[n+1,n+1] <- 0 pdepth <- kindepth(id, dadid, momid) mrow <- match(momid, id, nomatch=n+1) #row number of the mother drow <- match(dadid, id, nomatch=n+1) #row number of the dad for (depth in 1:max(pdepth)) { indx <- (1:n)[pdepth==depth] for (i in indx) { mom <- mrow[i] dad <- drow[i] kmat[i,] <- kmat[,i] <- (kmat[mom,] + kmat[dad,])/2 kmat[i,i] <- (1+ kmat[mom,dad])/2 } } kmat <- kmat[1:n,1:n] dimnames(kmat) <- list(id, id) kmat } oldkinship.pedigree <- function(id, ...) { n <- length(id$id)
if (n==1)
return(matrix(.5,1,1, dimnames=list(id$id, id$id)))
if (any(duplicated(id$id))) stop("All id values must be unique") kmat <- diag(n+1) /2 kmat[n+1,n+1] <- 0 pdepth <- kindepth(id) mrow <- ifelse(id$mindex ==0, n+1, id$mindex) drow <- ifelse(id$findex ==0, n+1, id$findex) # Are there any MZ twins to worry about? if (!is.null(id$relation) && any(id$relation$code=="MZ twin")) {
havemz <- TRUE
temp <- which(id$relation$code=="MZ twin")
## drop=FALSE added in case only one MZ twin set
mzmat <- as.matrix(id$relation[,c("indx1", "indx2")])[temp,,drop=FALSE] # any triples, quads, etc? if (any(table(mzmat) > 1)) { #yes there are # each group id will be min(member id) mzgrp <- 1:max(mzmat) #each person a group indx <- sort(unique(as.vector(mzmat))) # The loop below will take k-1 iterations for a set labeled as # 1:2, 2:3, ...(k-1):k; this is the worst case. while(1) { z1 <- mzgrp[mzmat[,1]] z2 <- mzgrp[mzmat[,2]] if (all(z1 == z2)) break mzgrp[indx] <- tapply(c(z1, z1, z2, z2), c(mzmat,mzmat), min) } # Now mzgrp = min person id for each person in a set matlist <- tapply(mzmat, mzgrp[mzmat], function(x) { x <- sort(unique(x)) temp <- cbind(rep(x, each=length(x)), rep(x, length(x))) temp[temp[,1] != temp[,2],] }) } else { #no triples, easier case matlist <- tapply(mzmat, row(mzmat), function(x) matrix(x[c(1,2,2,1)],2), simplify=FALSE) } } else havemz <- FALSE for (depth in 1:max(pdepth)) { indx <- (1:n)[pdepth==depth] for (i in indx) { mom <- mrow[i] dad <- drow[i] kmat[i,] <- kmat[,i] <- (kmat[mom,] + kmat[dad,])/2 kmat[i,i] <- (1+ kmat[mom,dad])/2 } if (havemz) { for (i in 1:length(matlist)) { temp <- matlist[[i]] kmat[temp] <- kmat[temp[1], temp[1]] } } } kmat <- kmat[1:n,1:n] dimnames(kmat) <- list(id$id, id$id) kmat } oldkinship.pedigreeList <- function(id, ...) { famlist <- unique(id$famid)
nfam <- length(famlist)
matlist <- vector("list", nfam)
idlist  <- vector("list", nfam) #the possibly reorderd list of id values

for (i in 1:length(famlist)) {
tped <- id[i]  #pedigree for this family
temp <- try(oldkinship(tped, ...), silent=TRUE)
if (class(temp)=="try-error")
stop(paste("In family", famlist[i], ":", temp))
else matlist[[i]] <- as(forceSymmetric(temp), "dsCMatrix")
idlist[[i]] <- tped$id } result <- bdiag(matlist) if (any(duplicated(id$id)))
temp <-paste(rep(famlist, sapply(idlist, length)),
unlist(idlist), sep='/')
else temp <- unlist(idlist)

dimnames(result) <- list(temp, temp)
result
}


text:
# Acceleration of a car?

1. Sep 29, 2010

### ScienceGirl90

1. The problem statement, all variables and given/known data
A car is moving at 50.0 m/s and brakes to a halt in 6.00 seconds.
(1)From the time it starts braking, how far does the car travel before the car comes to a halt if its acceleration is constant?
(2)What is the car's speed 1.35 s after it starts to brake if the cars acceleration is constant?

2. Relevant equations
v=a*t+vo

3. The attempt at a solution

2. Sep 29, 2010

### fss

The equation you listed will get you the acceleration, but for the distance you'll need

$$v_f^2 = v_0^2 + 2ad$$

3. Sep 29, 2010

### ScienceGirl90

Ok thank you! I was able to find the distance for question one but I am still stuck as to what is the car's speed 1.35 s after it starts to brake if the cars acceleration is constant?

4. Sep 29, 2010

### fss

Just use the first formula you gave in your original post. Pay attention to your signs and you'll be golden.

5. Sep 29, 2010

### ScienceGirl90

Got it. Thanks!


text:
# Problem of the Week Problem A and Solution Tasty Tarts

## Problem

Genevieve and Dima decided to make peach tarts. According to their recipe, each tart requires three-quarters of a peach. If they bought a dozen peaches, how many tarts can they make?

## Solution

Here is one way to solve the problem. We can draw a diagram of a dozen peaches where each peach is divided into four quarters. Then we group three quarters together with a common label, corresponding to the tart that those quarters are used in. We continue to identify groups of three quarters over the dozen peaches. Finally, we count the number of groups we have identified.

From this diagram, we can see that Genevieve and Dima can make $$16$$ tarts with a dozen peaches.


text:
# nLab ultraroot

Contents

model theory

## Dimension, ranks, forking

• forking and dividing?

• Morley rank?

• Shelah 2-rank?

• Lascar U-rank?

• Vapnik–Chervonenkis dimension?

# Contents

## Idea

An elementary substructure $N \preceq M$ of a model of a first-order theory is an ultraroot of $M$ if $M$ is isomorphic to an ultrapower of $N$.

## Definition

Formally, recall that any structure $N$ admits a diagonal elementary embedding $N \overset{\Delta}{\to} N^{\mathcal{U}}$ to any ultrapower $N^{\mathcal{U}}$ of $N$.

If $i \colon N \to M$ is an elementary embedding, we say that $N$ is an ultraroot of $M$ if there is some ultrafilter $\mathcal{U}$ on some index set such that there is an isomorphism

$\left(N \overset{i}{\to} M\right) \simeq \left(N \overset{\Delta}{\to} N^{\mathcal{U}}\right)$

in the co-slice category under $N$.

## Examples

• $N$ is an ultraroot of $N^{\mathcal{U}}$.

• Assuming the continuum hypothesis, any countable model of a complete-first order theory will be an ultraroot of that theory’s $\aleph_1$-saturated continuum-sized model.

• Similarly, taking a countably-indexed ultrapower of a slightly-expanded copy (like a nonstandard model generated by throwing in a single infinitesimal) of the reals still yields the hyperreal numbers.

## Remarks

Closure under taking ultraroots, along with closure under ultraproducts and elementary embeddings, gives a criterion for detecting elementary classes. See also at Birkhoff's HSP theorem.

## References

• See e.g. page 7 of these notes by Ward Henson on continuous model theory,

Last revised on February 24, 2017 at 07:48:15. See the history of this page for a list of all contributions to it.