text:
# Refactor the code which performs "cross-product", "reduce", "product of list" and "sum of list"

I have come up with a sequence of steps to find the maximum product of y positive numbers which add up to x. But the program is highly inefficient.

Here's my code:

from itertools import product
from functools import reduce
from operator import mul

for i in range(int(input())):               # <-- number of pairs of input
lgr = 1
k = set()
x, y = tuple(map(int, input().split())) #  <-- the input pairs
p = product(range(1, x+1), repeat=y)    #  <-- gets the cross product
for i in p:
if sum(i) == x:                     #  <-- checks if sum is equal to x
for i in k:
tmp = reduce(mul, i, 1)             #  <-- multiplies the items in the tuple to get the largest
if tmp > lgr:
lgr = tmp
print(lgr)

But say for an input like 14 7 it just takes too long (around 20s). Is there a better way to do it?

Update 1: I've managed to improve a bit by using combinations_with_replacement and set & list comprehension... but I believe it can be done better... any suggestions...?

from itertools import combinations_with_replacement
from functools import reduce
from operator import mul

for i in range(int(input())):
lgr = 1
x, y = tuple(map(int, input().split()))
k = {i for i in combinations_with_replacement(range(1, x+1), y) if sum(i) == x}
lgr = max([max(reduce(mul, i, 1), lgr) for i in k])
print(lgr)

• I was able to solve the problem for (14, 7) = 128 in 5 seconds without the help of a computer. The product is biggest when all numbers are equal, which is easy to prove. Apr 6 '20 at 11:24
• @Peilonrayz - Start with an array of all ones, then repeatedly add one to the smallest number (because that gives the largest percentual increase) until you reach the maximum allowed sum. This can be fomulated as a proof by induction. Apr 6 '20 at 14:01
• @Peilonrayz, if the product contains two factors which differ by more than one, then it is not optimal: Say these are x and y with x < y, then replace them by x + 1 and y - 1. Now (x+1)(y-1) = xy + (y - x - 1) > xy. Apr 6 '20 at 14:11
• I do not think that this kind of algorithmic question is right for code review SE. Apr 6 '20 at 14:31
• I agree with Carsten that this is not a code review, which is also why I started this thread in the comments. OP brute forced a trivial math exercise and I wanted to highlight that. Apr 6 '20 at 15:01

The comments contain two proofs that a better solution exists.

1. Start with an array of all ones, then repeatedly add one to the smallest number (because that gives the largest percentual increase) until you reach the maximum allowed sum. This can be fomulated as a proof by induction. – Rainer P.

This can be proven rather simply by using Python, but it can also be proven with a basic formula. If we start with the array $$\y, y-1\$$ and we have to increase one of the values by one then we can deduce adding one to the $$\y-1\$$ value results in a larger product.

$$yy \gt (y + 1)(y-1)\\ y^2 \gt y^2 - 1\\$$

This change works regardless of what the rest of the array is. And so we know adding one to the lowest value always nets the biggest increase.
I used to not be very good with math, and so a more hands on proof using Python may be easier to follow and agree with.

import math

def inc_1(numbers, pos):
nums = numbers[:]
nums[pos] += 1
return nums

def proof(total, amount):
if total < amount:
return 0
numbers = [1] * amount
for i in range(total - amount):
d = i % amount
if d:
larger = math.prod(inc_1(numbers, d - 1))
smaller = math.prod(inc_1(numbers, d))
if smaller < larger:
raise ValueError('proof does not hold')
numbers[d] += 1
return numbers, math.prod(numbers)

print(proof(14, 7))
print(proof(302, 5))

2. if the product contains two factors which differ by more than one, then it is not optimal: Say these are $$\x\$$ and $$\y\$$ with $$\x < y\$$, then replace them by $$\x + 1\$$ and $$\y - 1\$$. Now $$\(x+1)(y-1) = xy + (y - x - 1) > xy\$$. – Carsten S

Changed to use mathjax

With this description the factors are the array of numbers. With this we can start with a random assortment of numbers that total the desired number, and we can continuously improve the product.

$$3, 8\\ 3 < 8\\ 4 \times 7 \gt 3 \times 8\\ 28 \gt 24\\ 5 \times 6 \gt 4 \times 7\\ 30 \gt 24$$

In all the fastest solution is to just use divmod, and produce the product once.

def maximum_product_with_sum(total, count):
q, r = divmod(total, count)
return (q+1)**r * q**(count-r)

print(maximum_product_with_sum(14, 7))
print(maximum_product_with_sum(302, 5))

One immediate improvement is realizing that for e.g the input 14 7, all combinations like x x x x x x 9 and above can never sum to 14, because even for the minimal case x = 1, 9 + 6*1 = 15. So you can vastly reduce the number of combinations you need to check by using

k = {i for i in combinations_with_replacement(range(1, x - y + 2), y) if sum(i) == x}

For the input 14 7 this constitutes a difference of 77520 possible combinations down to 3432.

This adds a constraint only on the last number, but you could take this to its logical extreme and write code that always aborts the combinations once the sum has been reached.

One other thing you should do is use better names and use functions. I have no idea that k is a list of combinations of numbers, or what this means. Also, using i (and j, k, n, too, for that matter) for anything other than an integer is not helpful either.

I would at least rewrite your code to this:

from itertools import combinations_with_replacement
try:  # Python 3.8+
from math import prod
except ImportError:
from functools import reduce
from operator import mul
def prod(x):
"""Returns the product of all elements of x, similar to sum."""
return reduce(mul, x, 1)

def maximum_product_with_sum(x, n):
"""Return the maximum product of all n positive numbers that sum to x."""
r = range(1, x - n + 2)
candidates = {candidate
for candidate in combinations_with_replacement(r, n)
if sum(candidate) == x}
return max(map(prod, candidates), default=1)

if __name__ == "__main__":
for _ in range(int(input())):
x, n = tuple(map(int, input().split()))
print(maximum_product_with_sum(x, n))

Note that I also added some short docstrings to describe what each function does and simplified getting the maximum, your inner call to max was not needed. The optional keyword argument argument default of max is the default value if the iterable is empty.

I also added a if __name__ == "__main__": guard to allow importing from this script from another script without the code being run.

Here is the performance measured against your second function (the first one is way too slow) once for the input x 7 and once for 14 n:

You can see that this code is not only faster in general, it also becomes much faster if n > x / 2.

• If you're on 3.8 you can use math.prod Apr 6 '20 at 12:43
• @MaartenFabré: Ah, yes, I had forgotten they finally added that :) Edited to include it, albeit under a guard since quite a few systems don't have 3.8 yet as a default (Ubuntu 18.04 LTS, looking at you). Apr 6 '20 at 15:12

Thanks @Graipher & @Peilonrayz for taking time to answer my question. Yesterday around 7PM IST I discussed this with my brilliant friend Abhay, who gave me a clue on how this can be solved with just $$\O(n)\$$ complexity.

To keep the question in mind:

Find the maximum product of y positive numbers which adds up to x.

This is how he explained it to me:

If we observe the pattern of tuples that qualify for the output, we can easily see that all the values lie closer to $$\\frac{x}{y}\$$. See for the pairs:

| Pairs    | Solution Tuple                                | Max Product |
|----------|-----------------------------------------------|-------------|
| (14, 7)  | (2, 2, 2, 2, 2, 2, 2)                         | 128         |
| (20, 15) | (1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2) | 32          |
| (13, 4)  | (3, 3, 3, 4)                                  | 108         |
| (5, 3)   | (1, 2, 2)                                     | 4           |

If we take a even closer look we find that those values in the tuple are either ceil($$\\frac{x}{y}\$$) or floor($$\\frac{x}{y}\$$). And the last individual element is what is left off in sum.

Hence the program according to him was:

from math import floor, ceil

for t in range(int(input()):
x, y = map(int, input().split())
psum, product = 0, 1
for i in range(y):
if i == y-1: product *= x - psum
else:
num = x/y
num1, num2 = floor(num), ceil(num)
psums, product = (psum+num1, product*num1) if abs(num - num1) < abs(num - num2) else (psum+num2, product*num2)
print(product)

@RBarryYoung's comment came close but I think x's and y's got messed up... (please correct me if it is not so).

I've never seen how well this answer fares @Peilonrayz's solution, but I would love to get some insight on it.

• Peilonrayz's solution is the same as yours. solve(17, 5) = product(4, 4, 3, 3, 3) = 4² * 3³ = (3+1)² * 3³ where q = 3, r = 2 Apr 7 '20 at 10:19


text:
About Tyler Folkman: I am currently pursuing my masters in computer science at UT Austin with a focus on machine learning and data mining. I am doing research at the intersection of healthcare and data mining as a member of the IDEAL lab and work part-time at Cognitive Scale.
This tutorial was originally featured here, and is available on Github

So I feel like there is not a lot of good information out their on how to use Bayes Theorem for modeling – especially with Python code.  Like try figuring out how to understand a Bayesian Linear Regression from just Google searches – not super easy. So I thought I would maybe do a series of posts working up to Bayesian Linear regression.

The first post in this series is an introduction to Bayes Theorem with Python. I hope this post helps some understand what Bayes Theorem is and why it is useful. As well as get a small insight into how it differs from frequentist methods.

If anything isn’t clear or you have any comments, please let me know! Also, if you have any great Bayes resources, please share, and I can add them to the site.

### What is Bayes Theorem?

Bayes theorem is what allows us to go from a sampling (or likelihood) distribution and a prior distribution to a posterior distribution.

### What is a Sampling Distribution?

A sampling distribution is the probability of seeing our data (X) given our parameters (θ). This is written as $p(X|θ).$

For example, we might have data on 1,000 coin flips. Where 1 indicates a head. This can be represented in python as

import numpy as np
data_coin_flips = np.random.randint(2, size=1000)
np.mean(data_coin_flips)

Out[2]: 0.46800000000000003

A sampling distribution allows us to specify how we think these data were generated. For our coin flips, we can think of our data as being generated from a Bernoulli Distribution. This distribution takes one parameter p which is the probability of getting a 1 (or a head for a coin flip). It then returns a value of 1 with probablility p and a value of 0 with probablility (1-p).

You can see how this is perfect for a coin flip. With a fair coin we know our p = .5 because we are equally likely to get a 1 (head) or 0 (tail). We can create samples from this distribution like this:

bernoulli_flips = np.random.binomial(n=1, p=.5, size=1000)
np.mean(bernoulli_flips)

0.46500000000000002

Now that we have defined how we believe our data were generated, we can calculate the probability of seeing our data given our parameters $p(X|θ)$. Since we have selected a Bernoulli distribution, we only have one parameter: p.

We can use the probability mass function (PMF) of the Bernoulli distribution to get our desired probability for a single coin flip. The PMF takes a single observed data point and then given the parameters (p in our case) returns the probablility of seeing that data point given those parameters. For a Bernoulli distribution it is simple: if the data point is a 1 the PMF returns p, if the data point is a 0 it returns (1-p). We could write a quick function to do this:

def bern_pmf(x, p):
if (x == 1):
return p
elif (x == 0):
return 1 - p
else:
return "Value Not in Support of Distribution"

We can now use this function to get the probability of a data point give our parameters. You probably see that with p = .5 this function always returns .5

print(bern_pmf(1, .5))
print(bern_pmf(0, .5))

0.5
0.5

This is a pretty simple PMF, but other distributions can get much more complicated. So it is good to know that Scipy has most of these built in. We can draw from the PMF as follows:

import scipy.stats as st
print(st.bernoulli.pmf(1, .5))
print(st.bernoulli.pmf(0, .5))

0.5
0.5

This is nice, but what we really want to know is the probability of see all 1,000 of our data points. How do we do that? The trick here is to assume that our data are independent and identically distributed. This assumption allows us to say the probability of seeing all of our data is just the product of each individual probability: $p(x1,…,xn|β)=p(x1|β)∗…∗p(xn|β).$ This is easy to do:

np.product(st.bernoulli.pmf(data_coin_flips, .5))

9.3326361850321888e-302

How does that number help us? Well by itself, it doesn’t really help too much. What we need to do now is get more of a distribution for our sampling model. Currently, we have only tested our model with p = .5, but what if p = .8? or .2? What would the probablility of our data look like then? This can be done by defining a grid of values for our p. Below I will make a grid of 100 values between 0 and 1 (because p has to be between 0 and 1) and then I will calculate the probability of seeing our data given each of these values:

import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
sns.set(style='ticks', palette='Set2')

params = np.linspace(0, 1, 100)
p_x = [np.product(st.bernoulli.pmf(data_coin_flips, p)) for p in params]
plt.plot(params, p_x)
sns.despine()

Now we are getting somewhere. We can see that the probablility of seeing our data peaks at p=.5 and almost certainly is between p=.4 and p=.6. Nice. So now we have a good idea of what p value generated our data assuming it was drawn from a Bernoulli distribution. We’re done, right? Not quite…

### Prior Distribution

Bayes theorem says that we need to think about both our sampling distribution and our prior distribution. What do I mean by prior distribution? It is the $p(θ)$ or the probability of seeing a specific value for our parameter. In our sampling distribution we defined 100 values from 0 to 1 for our parameter p. Now we must define the prior probability of seeing each of those values. That is the probability we would have assumed before seeing any data. Most likely, we would have assumed a fair coin, which looks like the distribution above. Lets see how we can do this:

fair_flips = bernoulli_flips = np.random.binomial(n=1, p=.5, size=1000)
p_fair = np.array([np.product(st.bernoulli.pmf(fair_flips, p)) for p in params])
p_fair = p_fair / np.sum(p_fair)
plt.plot(params, p_fair)
sns.despine()

Basically we created 1,000 fair coin flips and then generated the sampling distribution like we did before (except we divided by the sum of the sampling distribution to make the values sum to 1). Now we have a “fair coin” prior on our parameters. This basically means that before we saw any data we thought coin flips were fair. And we can see that assumption in our prior distribution by the fact that our prior distribution peaks at .5 and is almost all between .4 and .6.

I know what you are thinking – this is pretty boring. The sampling and prior distributions look exactly the same. So lets mix things up. Lets keep our fair prior but change our data to be an unfair coin:

unfair_flips = bernoulli_flips = np.random.binomial(n=1, p=.8, size=1000)
p_unfair = np.array([np.product(st.bernoulli.pmf(unfair_flips, p)) for p in params])
fig, axes = plt.subplots(2, 1, sharex=True)
axes[0].plot(params, p_unfair)
axes[0].set_title("Sampling Distribution")
axes[1].plot(params, p_fair)
axes[1].set_title("Prior Distribution")
sns.despine()
plt.tight_layout()

Ah – now this is interesting. We have strong data evidence of an unfair coin (since we generated the data we know it is unfair with p=.8), but our prior beliefs are telling us that coins are fair. How do we deal with this?

### Bayes Theorem (Posterior Distribution)

Bayes theorem is what allows us to go from our sampling and prior distributions to our posterior distribution. The posterior distribution is the $P(θ|X)$. Or in English, the probability of our parameters given our data. And if you think about it that is what we really want. We are typically given our data – from maybe a survey or web traffic – and we want to figure out what parameters are most likely given our data. So how do we get to this posterior distribution? Here comes some math (don’t worry it is not too bad):

By definition, we know that (If you don’t believe me, check out this page for a refresher):

• $P(A|B) = \dfrac{P(A,B)}{P(B)}$ Or in English, the probability of seeing A given B is the probability of seeing them both divided by the probability of B.
• $P(B|A) = \dfrac{P(A,B)}{P(A)}$ Or in English, the probability of seeing B given A is the probability of seeing them both divided by the probability of A.

You will notice that both of these values share the same numerator, so:

• $P(A,B) = P(A|B)*P(B)$
• $P(A,B) = P(A|B)*P(B)$

Thus:

$P(A|B)*P(B) = P(B|A)*P(A)$

Which implies:

$P(A|B) = \dfrac{P(B|A)*P(A)}{P(B)}$

And plug in $θ$ for $A$ and $X$ for $B$:

$P(\theta|X) = \dfrac{P(X|\theta)*P(\theta)}{P(X)}$

Nice! Now we can plug in some terminology we know:

$Posterior = \dfrac{likelihood * prior}{P(X)}$

But what is the $P(X)?$ Or in English, the probability of our data? That sounds weird… Let’s go back to some math and use B and A again:

We know that $P(B)=∑AP(A,B)$ (check out this page for a refresher)

And from our definitions above, we know that:

$P(A,B) = P(A|B)*P(A)$

Thus:

$P(B) = \sum_{A} P(A|B)*P(A)$

Plug in our $θ$ and $X$:

$P(X) = \sum_{\theta} P(\theta|X)*P(\theta)$

Plug in our terminology:

$P(X) = \sum_{\theta} likelihood * prior$

Wow! Isn’t that awesome! But what do we mean by $∑θ$. This means to sum over all the values of our parameters. In our coin flip example, we defined 100 values for our parameter p, so we would have to calculated the likelihood * prior for each of these values and sum all those anwers. That is our denominator for Bayes Theorem. Thus our final answer for Bayes is:

$Posterior = \dfrac{likelihood * prior}{\sum_{\theta} likelihood * prior}$

That was a lot of text. Let’s do some more coding and put everything together.

def bern_post(n_params=100, n_sample=100, true_p=.8, prior_p=.5, n_prior=100):
params = np.linspace(0, 1, n_params)
sample = np.random.binomial(n=1, p=true_p, size=n_sample)
likelihood = np.array([np.product(st.bernoulli.pmf(sample, p)) for p in params])
#likelihood = likelihood / np.sum(likelihood)
prior_sample = np.random.binomial(n=1, p=prior_p, size=n_prior)
prior = np.array([np.product(st.bernoulli.pmf(prior_sample, p)) for p in params])
prior = prior / np.sum(prior)
posterior = [prior[i] * likelihood[i] for i in range(prior.shape[0])]
posterior = posterior / np.sum(posterior)

fig, axes = plt.subplots(3, 1, sharex=True, figsize=(8,8))
axes[0].plot(params, likelihood)
axes[0].set_title("Sampling Distribution")
axes[1].plot(params, prior)
axes[1].set_title("Prior Distribution")
axes[2].plot(params, posterior)
axes[2].set_title("Posterior Distribution")
sns.despine()
plt.tight_layout()

return posterior
example_post = bern_post()

You will notice that I set 100 as the number of observations for the prior and likelihood. This increases the variance of our distributions. More data typically decreases the spread of a distribution. Also, as you get more data to estimate your likelihood, the prior distribution matters less.

moredata_post = bern_post(n_sample=1000)

You can see that effect in the graphs above. Because we have more data to help us estimate our likelihood our posterior distribution is closer to our likelihood. Pretty cool.

### Conclusion

There you have it. An introduction to Bayes Theorem. Now if you ever are doubting the fairness of a coin you know how to investigate that problem! Or maybe the probability of a population voting yes for a law? Or any other yes/no outcome. I know what you are thinking – not that impressive. You want to predict things or run a fancy algorithm… in time my friends. I actually would like to write a series of posts leading up to Bayesian Linear Regression. Hopefully this is the first post in that series :)

#### Some Side Notes

1. You will notice that the denominator for Bayes Theorem is just a constant. So if you only want to get the maximum posterior value, you don’t even need to calculate that constant. For this reason you will often see the posterior shown as proportional (or in math) to the likelihood * prior.
2. Frequentist statistics is focused on the likelihood. Or you could say that frequentists are bayesians with a non-informative prior (like a uniform distribution). But don’t hate on frequentists too much; most of bayesian inference in applied settings relies of frequentists statistics.
3. Now that you know frequentist statistics focuses on the likelihood it is much clearer why people often misinterpret frequentist confidence intervals. The likelihood is P(X|θ) – or the probability of our data given our parameters. That is a bit wierd because we are given our data, not our parameters. What most frequentists models do is take the maximum of the likelihood distribution (or Maximum Likelihood Estimation (MLE)). Basically find what parameters maximize the probability of seeing our data. The important point here is that they are treating the data as random, so what the frequentist confidence interval is saying is that if you were to keep getting new data (maybe some more surveys) and calculate confidence intervals on each of these new samples, 95% of these samples would have confidence intervals that contain the true parameter you are trying to estimate.

(Featured image credit: Mattbuck, via Wikimedia Commons)

format: plain


text:
# Adjusting variable distribution and exploring data using mass linear regression

Graphs and analysis using the #TidyTuesday data set for week 33 of 2021 (10/8/2021): “BEA Infrastructure Investment”

Ronan Harrington https://github.com/rnnh/
08-15-2021

## Introduction

In this post, the BEA Infrastructure Investment data set from the #TidyTuesday project is used to illustrate variable transformation and the exploreR::masslm() function. The variable for gross infrastructure investment adjusted for inflation is transformed to make it less skewed. Using these transformed investment values, multiple linear models are then created to quickly see which variables in the data set have the largest impact on infrastructure investment.

## Setup

Loading the R libraries and data set.

Show code
# Loading libraries
library(tidyverse)
library(tidytuesdayR)
library(exploreR)

Downloading file 1 of 3: ipd.csv
Downloading file 2 of 3: chain_investment.csv
Downloading file 3 of 3: investment.csv

## Plotting distribution of inflation-adjusted infrastructure investments

In this section, the gross infrastructure investment (chained 2021 dollars) in millions of USD are plotted with and without a $$\log{10}$$ transformation. From the histograms below, we can see that applying a $$\log{10}$$ transformation gives the variable a less skewed distribution. This transformation should be considered for statistical testing of inflation-adjusted infrastructure investments.

Show code
# Creating tbl_df with gross_inv_chain values
untransformed_tbl_df <- tibble(
gross_inv_chain = tt$chain_investment$gross_inv_chain,
transformation = "Untransformed"
)

# Creating tbl_df with log10(gross_inv_chain) values
log10_tbl_df <- tibble(
gross_inv_chain = log10(tt$chain_investment$gross_inv_chain),
transformation = "Log10"
)

# Combining the above tibbles into one tbl_df
gross_inv_chain_tbl_df <- rbind(untransformed_tbl_df, log10_tbl_df)

# Plotting distribution of inflation-adjusted infrastructure investments
gross_inv_chain_tbl_df %>%
ggplot(aes(x = gross_inv_chain, fill = transformation)) +
geom_histogram(show.legend = FALSE, position = "identity",
bins = 12, colour = "black") +
facet_wrap(~transformation, scales = "free") +
labs(fill.position = "none", y = NULL,
x = "Gross infrastructure investments adjusted for inflation (millions USD)",
title = "Distributions of untransformed and log transformed infrastructure investments",
subtitle = "Log transformed investments are more normally distributed") +
scale_fill_brewer(palette = "Set1") +
theme_classic()

## Exploring a data set using mass linear regression

In this section, exploreR::masslm() is applied to a copy of the data set with $$\log{10}$$ transformed investment values. The masslm() function from the exploreR package quickly produces a linear model of the dependent variable and every other variable in the data set. It then returns a data frame containing the features of each linear model that are useful when selecting predictor variables:

• R squared The proportion of variation in the dependent (response) variable that is explained by the independent (predictor) variable.
• p-value The statistical significance of the model. A p-value $$\lt 5\%$$ is typically considered significant.

This function is useful for quickly determining which variables should be included in predictive models. Note that the data set used should satisfy the assumptions of linear models, including a normally distributed response variable. In this case, the $$\log{10}$$ transformed investment variable is close to normal.

From this mass linear regression model, we can see that investment category is the single variable that explains the largest proportion of variation in $$\log{10}$$ investment; and the linear model with group number is the most significant, followed by year.

Show code
# Creating a copy of the chain_investment data set with log10 transformed
# gross investment values
chain_investment_df <- tt\$chain_investment %>%
# Creating a log10 transformed copy of gross_inv_chain
mutate(gross_inv_transformed = log10(gross_inv_chain)) %>%
# Removing -Inf values
filter(gross_inv_transformed != -Inf) %>%
# Selecting variables to include in the data frame
select(category, meta_cat, group_num, year, gross_inv_transformed)

# Applying mass linear regression
transformed_investment_masslm <- masslm(chain_investment_df,
dv.var = "gross_inv_transformed")

# Printing the masslm results in order of R squared values (decreasing)
transformed_investment_masslm %>%
arrange(-R.squared)

         IV Coefficient    P.value  R.squared
1  category   -0.579900  8.471e-10 0.63754622
2  meta_cat    0.349300  7.848e-10 0.37782201
3 group_num   -0.058750 3.625e-204 0.14695670
4      year    0.009507  7.007e-59 0.04377399
Show code
# Printing the masslm results in order of p-values
transformed_investment_masslm %>%
arrange(P.value)

         IV Coefficient    P.value  R.squared
1 group_num   -0.058750 3.625e-204 0.14695670
2      year    0.009507  7.007e-59 0.04377399
3  meta_cat    0.349300  7.848e-10 0.37782201
4  category   -0.579900  8.471e-10 0.63754622

## References

### Corrections

If you see mistakes or want to suggest changes, please create an issue on the source repository.

### Reuse

Text and figures are licensed under Creative Commons Attribution CC BY 4.0. Source code is available at https://github.com/rnnh/TidyTuesday/, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".

### Citation

For attribution, please cite this work as

Harrington (2021, Aug. 15). Ronan's #TidyTuesday blog: Adjusting variable distribution and exploring data using mass linear regression. Retrieved from https://tidytuesday.netlify.app/posts/2021-08-15-bea-infrastructure-investment/

BibTeX citation

@misc{harrington2021adjusting,
author = {Harrington, Ronan},
title = {Ronan's #TidyTuesday blog: Adjusting variable distribution and exploring data using mass linear regression},
url = {https://tidytuesday.netlify.app/posts/2021-08-15-bea-infrastructure-investment/},
year = {2021}
}


text:
A bolt pattern is arrangement of bolted joints, typically four or more, that connect two or more components together. When designing a bolt pattern, it is important to have an understanding of the loads that the pattern will need to resist in operation as well as the way in which those applied loads will distribute among the individual bolts in the pattern. Once the axial and shear loads have been calculated for the individual bolted joints in the pattern, the individual joints can be analyzed as discussed here.

The following sections detail the methodology to resolve forces and moments applied to a bolt pattern into axial and shear loads acting on the individual bolted joints in the pattern.

Contents

## Overview and Discussion

The distribution of a forces and moments over a bolt pattern is similar to the analysis of a beam or a shaft. Applied loads are translated to the centroid of the pattern (analagous to the neutral axis of a beam or shaft). The forces and moments at the centroid are then resolved into axial and shear forces acting at the individual bolted joints. Axial forces are distributed over a bolt pattern based on pattern's area, $$A$$, and moments of inertia, $$I_{c.x}$$ and $$I_{c.y}$$. Likewise, shear forces are distributed based on the pattern's area, $$A$$, and polar moment of inertia, $$I_{c.p}$$.

Note that we maintain the distinction that the forces are distributed to individual bolted joints rather than individual bolts. The reason for this is that not all of the load applied to a bolted joint is actually seen by the bolt, as discussed here.

Two common loading conditions are discussed in the next two sections, followed by a discussion of the generalized approach for distributing applied loads over a bolt pattern.

The figure below shows a pattern with an eccentric shear load that is applied in the plane of the pattern. The load is eccentric because it does not act through the centroid of the pattern. Therefore it induces a torsional moment about the Z-axis (perpendicular to the plane of the pattern) that will tend to rotate the pattern about its centroid. In this case the bolts are loaded in shear only (no axial loads), and the shear loads are due to a comination of the direct shear force and the induced moment. Loading of this type is analagous to a shaft under a combined shear and torsion load.

The figure below shows a pattern with an eccentric shear load that is applied out of the plane of the pattern. While this applied load passes over the centroid in the X-Y plane, its line of action is offset from the centroid in the Z-direction. Therefore the applied shear load induces a bending moment about the X-axis which results in axial loading on the bolts.

There are several standard approaches to distributing axial loads among the bolts in a case like this, all of which involve calculating the moment of inertia of the pattern about some bending axis and then using $$(Mr/I) \cdot A$$ to distribute the loads. (It should be noted that if all of the bolts in the pattern are the same size, then $$(Mr/I) \cdot A$$ simplifies to $$Mr / \sum r^2$$ ). The key difference in the standard approaches is the selection of the point about which the pattern is assumed to pivot:

• One approach is to inspect the component geometry and to assume that the pattern will "heel" about some reasonable point, say the base of the bracket in the figure above.
• Another approach is to treat this problem as if the pattern were a beam in bending and to calculate an appropriate location of the neutral axis.
• The methodology outlined in the AISC "Steel Construction Manual" recognizes that on the compression side of the beam, the compressed area is equal to the contact area of the components below the neutral axis. On the tension side of the beam, only the bolts are in tension. The location of the neutral axis is found through an iterative procedure by which the prospective location of the neutral axis is varied until the moment of inertia of the bolts in tension above the neutral axis is equal to the moment of inertia of the compressed plates below the neutral axis.
• A simpler approach considers only the bolt areas and assumes that the neutral axis of the beam lies at the centroid of the bolt pattern.

It is important to recognize that as the assumed pivot location moves farther from the pattern, the axial loads on the bolts are reduced. The reason for this is that the axial loads are proportional to $$Mr / \sum r^2$$, where $$r$$ is the distance between the pivot location and the bolt of interest. As $$r$$ increases, the affect of the $$r^2$$ term in the denominator outpaces the affect of the $$r$$ term in the numerator. Therefore, the most conservative approach is to consider that the pattern pivots about its centroid.

### Generalized Approach

If you decide to take the conservative approach of translating all applied loading to the centroid of the bolt pattern, then it becomes easy to generalize the analysis of any bolt pattern with any applied loading, as in the figure below.

The steps to distribute applied forces and moments to the individual bolted joints in the pattern are:

1. Calculate relevant pattern properties (area, centroid, moments of inertia).
2. Translate all applied forces and moments to the centroid of the pattern.
3. Calculate axial and shear loads acting on individual bolted joints in the pattern.

Details to perform each of these steps are provided in the following sections.

### Need a Bolt Pattern Calculator?

Try our bolt pattern force distribution calculator, which allows for applied forces to be distributed over bolts in a pattern.

## Bolt Pattern Properties

The same properties that are required when analyzing a beam or a shaft are also required when distributing forces over a bolt pattern. It should be noted that the equations presented in this section for calculating pattern properties do not require the bolts in the pattern to be the same size.

### Bolt Pattern Area

The combined area of all bolts in the pattern must be calculated in order to distribute direct forces among the bolts:

$$A_{cmb} = \sum_i A_i$$

where $$A_i$$ is the tensile stress area of an individual bolt.

### Bolt Pattern Centroid

Just as bending stresses in a beam and torsional stresses in a shaft are centered about the neutral axis, moments on a bolt pattern will tend to rotate the pattern about its centroid. The location of the pattern centroid is calculated as:

 $$x_c = { \sum_i x_i A_i \over \sum_i A_i }$$ $$y_c = { \sum_i y_i A_i \over \sum_i A_i }$$

where $$A_i$$ is the bolt area and $$x_i$$ and $$y_i$$ are the x- and y- bolt locations, respectively.

It should be noted that this calculation is directly analogous to the calculation of the centroid of a cross section.

### Bolt Pattern Moments of Inertia

The moments of inertia of a bolt pattern indicate the ability of the pattern to resist bending moments. We conservatively assume that moments will tend to cause the pattern to rotate about its centroid, so moments of inertia about the pattern centroid are of interest. The centroidal moments of inertia are calculated as:

 $$I_{c.x} = \sum_i r_{c.y,i}^2 A_i$$ centroidal moment of inertia about the X-axis $$I_{c.y} = \sum_i r_{c.x,i}^2 A_i$$ centroidal moment of inertia about the Y-axis

where $$A_i$$ is the bolt area and $$r_{c.x,i}$$ and $$r_{c.y,i}$$ are the x- and y- distances of the bolt from the pattern centroid, respectively.

The polar moment of inertia of the pattern indicates the pattern's ability to resist torsional moments (i.e. moments about the Z-axis perpendicular to the plane of the pattern), and it is calculated as:

$$I_{c.p} = \sum_i r_{c.xy,i}^2 A_i = \sum_i ( r_{c.x,i}^2 + r_{c.y,i}^2 ) A_i = I_{c.x} + I_{c.y}$$

where $$r_{c.xy,i}$$ is the shortest distance between the bolt and the centroid and is calculated as $$r_{c.xy,i} = \sqrt{ r_{c.x,i}^2 + r_{c.y,i}^2 }$$.

It should be noted that these calculations are directly analogous to the calculation of the moment of inertia of a cross section and polar moment of inertia of a cross section.

## Forces and Moments at Centroid

As discussed previously, all applied forces and moments are translated to the centroid of the bolt pattern. As indicated in the figure below, any number of forces can be applied to the bolt pattern at any location.

The forces at the centroid are calculated as the sum of all applied forces:

$$\overline{F}_c = \sum_i \overline{F}_i$$

The moments at the centroid are calculated as the sum of all applied moments, plus the sum of the cross product of each applied force with the vector from the centroid to the location of that applied force:

$$\overline{M}_c = \sum_i \overline{M}_i + \sum_i \left( \overline{R}_{c,i} \times \overline{F}_i \right)$$

In the equations above, the bars over the variables indicate that they are vectors. The variable $$\overline{F}$$ is a force vector composed of the force components in each direction: $$F_x$$, $$F_y$$, and $$F_z$$. Likewise, $$\overline{M}$$ is a moment vector composed of moments about each axis. $$\overline{R}$$ is a location vector specifying the location of an applied force with respect to the pattern centroid. The location vector $$\overline{R}$$ points from the centroid to the location of the applied force.

### Need a Bolt Pattern Calculator?

Try our bolt pattern force distribution calculator, which allows for applied forces to be distributed over bolts in a pattern.

## Forces on Individual Bolted Joints

Once the pattern properties are known and the applied forces and moments have been translated to the pattern centroid, it is possible to calculate the axial and shear forces on the individual bolted joints. The figure below shows an individual bolted joint with applied axial and shear loading:

### Axial Forces

The axial forces are a result of the direct force in the Z-direction, $$F_{c.z}$$, the centroidal moment about the X-axis, $$M_{c.x}$$, and the centroidal moment about the Y-axis, $$M_{c.y}$$, as shown in the figure below:

The direct force in the Z-direction, $$F_{c.z}$$, is divided between the individual bolted joints according to the bolt stiffnesses. Because the bolts are all assumed to have the same material and length, the stiffness is dependent only on area. The axial force on a bolted joint due to the direct force in Z is calculated as:

$$P_{z.FZ} = { F_{c.z} A \over \sum_i A_i }$$

where $$A$$ is the area of the bolt in question. If the bolt areas are the same, the equation above simplifies to $$P_{z.FZ} = F_{c.z} / n$$, where $$n$$ is the number of bolts in the pattern.

The axial forces on a bolt due to moments about X- and Y- axes are calculated as:

 $$P_{z.MX} = { M_{c.x} r_{c.y} \over I_{c.x} } \cdot A$$ axial force on bolt due to MX about centroid $$P_{z.MY} = { -M_{c.y} r_{c.x} \over I_{c.y} } \cdot A$$ axial force on bolt due to MY about centroid

where $$M_{c.x}$$ and $$M_{c.y}$$ are the centroidal moments about the X- and Y- axes, $$r_{c.x}$$ and $$r_{c.y}$$ are the bolt distances from the centroid in the X- and Y-directions, and $$I_{c.x}$$ and $$I_{c.y}$$ are the pattern moments of inertia about the X- and Y- axes.

If the bolt areas are the same, the equations above simplify to:

 $$P_{z.MX} = { M_{c.x} r_{c.y} \over \sum_i r_{c.y,i}^2 }$$ $$P_{z.MY} = { -M_{c.y} r_{c.x} \over \sum_i r_{c.x,i}^2 }$$

The total axial force on a bolt is the sum of the axial force components:

$$P_{axial} = P_{z.FZ} + P_{z.MX} + P_{z.MY}$$

### Shear Forces

The shear forces are a result of the direct force in the X-direction, $$F_{c.x}$$, the direct force in the Y-direction, $$F_{c.y}$$, and the centroidal moment about the Z-axis, $$M_{c.z}$$, as shown in the figure below:

The direct forces in the X- and Y- directions, $$F_{c.x}$$ and $$F_{c.y}$$, respectively, are divided between the bolts according to the bolt stiffnesses. Because the bolts are all assumed to have the same material and length, the stiffness is dependent only on area. The shear reactions on a bolt due to the direct forces in X- and Y- are calculated as:

 $$P_{x.FX} = { -F_{c.x} A \over \sum_i A_i }$$ X-reaction on bolt due to direct force in X $$P_{y.FY} = { -F_{c.y} A \over \sum_i A_i }$$ Y-reaction on bolt due to direct force in Y

where $$A$$ is the area of the bolt in question. If the bolt areas are the same, the equations above simplify to $$P_{x.FX} = F_{c.x} / n$$ and $$P_{y.FX} = F_{c.y} / n$$, where $$n$$ is the number of bolts in the pattern.

The shear reaction on a bolt due to moment about the Z-axis is calculated as:

$$P_{xy.MZ} = { M_{c.z} r_{c.xy} \over I_{c.p} } \cdot A$$

where $$M_{c.z}$$ is the centroidal moment about the Z-axis and $$I_{c.p}$$ is the pattern's polar moment of inertia. The value $$r_{c.xy}$$ is the shortest distance between the bolt and the centroid and is calculated as $$r_{c.xy} = \sqrt{ r_{c.x}^2 + r_{c.y}^2 }$$.

The shear reaction $$P_{xy.MZ}$$ is then resolved into X- and Y- components based on the angle $$\theta$$ (see the figure above):

 $$P_{x.MZ} = P_{xy.MZ} \cdot \sin{ \theta }$$ X-reaction on bolt due to MZ about centroid $$P_{y.MZ} = -P_{xy.MZ} \cdot \cos{ \theta }$$ Y-reaction on bolt due to MZ about centroid

The value $$\theta$$ is the angle between the bolt location and the positive X-axis and is calculated as $$\theta = \tan^{-1}({ r_{c.y} / r_{c.x} })$$.

The total shear reaction on a bolt is calculated as the vector sum of the X- components plus the Y- components:

$$P_{shear} = \sqrt{ (P_{x.FX} + P_{x.MZ})^2 + (P_{y.FY} + P_{y.MZ})^2 }$$

## What's Next?

Now that the axial and shear forces on the individual bolted joints have been calculated, the stresses in the bolted joints can be analyzed as discussed here.

### Need a Bolt Pattern Calculator?

Try our bolt pattern force distribution calculator, which allows for applied forces to be distributed over bolts in a pattern.


text:
# Tag Info

Accepted

### Why is the log in the big-O of binary search not base 2?

When you change the base of logarithm the resulting expression differs only by a constant factor which, by definition of Big-O notation, implies that both functions belong to the same class with ...
• 9,602
Accepted

### Why is binary search faster than ternary search?

If you apply binary search, you have $$\log_2(n)+O(1)$$ many comparisons. If you apply ternary search, you have $$2 \cdot \log_3(n) + O(1)$$ many comparisons, as in each step, you need to perform 2 ...
• 2,682
Accepted

### How is this sorting algorithm Θ(n³) and not Θ(n²), worst-case?

This algorithm can be re-written like this Scan A until you find an inversion. If you find one, swap and start over. If there is none, terminate. Now there can be ...
• 70.9k
Accepted

### Will hardware/implementation affect the time/space complexity of algorithms?

Sure. Certainly. Here's how to reconcile your discomfort. When we analyze the running time of algorithms, we do it with respect to a particular model of computation. The model of computation ...
• 141k

### How is algorithm complexity modeled for functional languages?

If you assume that the $\lambda$-calculus is a good model of functional programming languages, then one may think: the $\lambda$-calculus has a seemingly simple notion of time-complexity: just count ...
• 8,178

### Why is binary search faster than ternary search?

DCTLib is right, but forget the math for a second. By your logic then, n-ary should be the fastest. But if you think about it, n-ary is exactly equal to a regular iteration search (just iterating ...
• 421
Accepted

### How long does the Collatz recursion run?

This is Collatz conjecture - still open problem. Conjecture is about proof that this sequence stops for any input, since this is unresolved, we do not know how to solve this runtime recurrence ...
• 9,325

### Checking equality of integers: O(1) in C but O(log n) in Python 3?

Integers are just binary strings and, to determine equality, both languages will compare the strings bit-by-bit. Not quite. C ints are machine-word-sized and ...
• 3,107
Accepted

### Why does Randomized Quicksort have O(n log n) worst-case runtime cost

Both of your sources refer to the "worst-case expected running time" of $O(n \log n).$ I'm guessing this refers to the expected time requirement, which differs from the absolute worst case. Quicksort ...
• 737

### How is algorithm complexity modeled for functional languages?

Algorithm complexity is designed to be independent of lower level details. No, not really. We always count elementary operations in some machine model: Steps for Turing machines. Basic operations on ...
• 70.9k
Accepted

### Why doesn't Knuth's linear-time multiplication algorithm "count"?

While the algorithm you mention appears in Knuth's TAOCP, it is certainly not due to Knuth, and is more widely known as the Schönhage–Strassen algorithm; Knuth even attributes this algorithm to them ...
• 270k
Accepted

### Are there any algorithms or data structures that need to find the median value of a set?

if there are any practical applications of this algorithm in the domain of computer science besides being a theoretical improvement The application of this algorithm is trivial - you use it whenever ...
• 9,602

### Why aren't primality tests easily linear in time complexity?

It seems that the main sticking point of the question here is: Why express runtime in terms of the size of the input, rather than the numeric value that the input represents? And indeed in some cases ...

### Checking equality of integers: O(1) in C but O(log n) in Python 3?

Complexity is defined relative to a computation model. P and NP, for example, are defined in terms of Turing machines. For comparison, consider the word RAM model. In this model, memory is divided ...
• 18.9k
Accepted

### Matrix Max in less than O(n)

If you don't know anything about the contents of the matrix (such as some kind of monotonicity property), linear time is the best you can do for a one-off search with a deterministic algorithm by a ...
• 1,999
Accepted

### Why is the dynamic programming algorithm of the knapsack problem not polynomial?

When we say polynomial or exponential, we mean polynomial or exponential in some variable. $nW$ is polynomial in $n$ and $W$. However, we usually consider the running time of an algorithm as a ...

### How long does the Collatz recursion run?

You translated the code correctly. There are many methods for solving recurrences. However, it is currently unknown if collatz even halts for all ...
• 70.9k

### how to calculate time complexity of non terminating loops

If the program runs forever, its running time is infinite. So, if it always enters an infinite loop, its running time is infinite. This is a degenerate case. Normally we focus only on algorithms ...
• 141k
Accepted

### A* graph search time-complexity

These are basically two different perspectives or two different ways of viewing the running time. Both are valid (neither is incorrect), but $O(b^d)$ is arguably more useful in the settings that ...
• 141k
Accepted

### What's the Big O runtime of a DFS word search through a matrix?

The complexity will be $O(m*n*4^{s})$ where m is the no. of rows and n is the no. of columns in the 2D matrix and s is the length of the input string. When we start searching from a character we ...
• 1,135
Accepted

### Proving the lower bound of compares in comparison based sorting

A sorting algorithm using at most $h$ comparisons on all inputs corresponds to a tree of height at most $h$. Such a tree has at most $2^h$ leaves. On the other hand, each permutation of $1,\ldots,N$ ...
• 270k

### Is there a method for automatic runtime analysis of algorithms?

No algorithm can decide whether a given algorithm ever halts or not, so in particular no algorithm can tightly analyze the complexity of a given algorithm.
• 270k
Accepted

### Is there a method for automatic runtime analysis of algorithms?

The COSTA tool does just this, although it fails in many cases, as you can imagine, due to computability problems. There are many papers about this; Cost Analysis of Java Bytecode by E. Albert, P. ...
• 20k
Accepted

### If recursive Fibonacci is $O(2^N)$ then why do I get 15 calls for N=5?

You ask, I have an $O(2^n)$ runtime, why do I not observe $2^n$ recursive calls for $n=15$? There are many things wrong in the implied conclusion. $O(\_)$ only gives you an upper bound. The true ...
• 70.9k

### How long does the Collatz recursion run?

The time complexity function is \begin{cases} T(n)= O(1) \text{ for } n\le 1\\ T(n)=T(n/2) + O(1) \text{ for } n\text{ even}\\ T(n)=T(3n+1) + O(1)\text{ for } n\text{ odd}\\ \end{cases} which can be ...
• 4,747

### Are there any algorithms or data structures that need to find the median value of a set?

Median filtering is common in reduction of certain types of noise in image processing. Especially salt and pepper noise. It works by picking out the median value in each color channel in each local ...
• 231
Accepted

### How can merging two sorted arrays of N items require at least 2N - 1 comparisons in every case?

What is meant by "lower bound" in this case is a lower bound on the worst-case number of comparisons. In this case, it happens to also be an upper bound. The lower bound has to be something that is ...

### Time Complexity proof for Segment Tree implementation of the ranged sum problem

The claim is that there are at most $2$ nodes which are expanded at each level. We will prove this by contradiction. Consider the segment tree given below. Let's say that there are $3$ nodes that ...
• 243
Accepted

### Tree decomposition - Fastest algorithm in practise

The first caveat is that deciding whether the treewidth of a graph is at most $t$ is NP-complete, but it is FPT (which is what Bodlaender's paper shows). So for small $t$, we can (in principle) solve ...
• 17.7k
Consider the following (silly) function cow, which accepts a number $n$: ...


text:
# Has anyone studied spaces based on the number of centers circles have on them?

I have been away from math for a while so be gentle if this is not very rigorous or if I am redefining already defined objects. The essential question that started me on this was the following. In the plane every circle has one center. On the surface of the sphere every circle has two centers. Okay, so let's formalize this notion.

Let $$(X, d)$$ be a metric space and let $$C_r(x) = \{y \in X \ | \ d(x-y) = r \}$$. Let $$D \subseteq \mathbb{R}$$ be the set of all distances between any two points of $$X$$. Then consider the following function:

$$f: X \times D \rightarrow P(X)$$

$$(x, r) \mapsto \{ z \in X \ | \ \exists s \in D \text{ with } C_s(z) = C_r(x) \}$$

We know that $$x \in f(x,r), \forall (x,r) \in X \times D$$. So this set is never empty and contains the original point. It is also reflexive in the sense that if $$y \in f(x,r)$$ for some $$r \in D$$ then $$x \in f(y,s)$$ for some $$s \in D$$. It would be nice if we could make this an equivalence relation. We can say that a metric space $$(X,d)$$ is regular if $$f(x,r) = f(x,s) \ \forall r,s \in \mathbb{D}$$. For a regular space $$(X,d)$$ we can define the following equivalence relation:

$$x \sim y \iff \exists r,s \in D \text{ with } C_r(x) = C_s(y)$$

Returning to our examples:

a) $$\mathbb{R}^n$$ with the metric induced by any $$\mathbb{L}_p$$ norm. Every equivalence class contains only the point itself.

b) The surface of a sphere with the "shortest path on the surface" distance. The equivalence classes are the sets of antipodal points.

c) A cone. I believe the center of the base and the top of the cone are in the same equivalence class, while every other point is in its own class.

We can say that a regular metric space is uniform if every equivalence class is of equal size. For a regular uniform metric space we can call the size of the equivalence class the degree of it. So that the first two examples would be regular uniform metric spaces but the 3rd would be regular but not uniform. The degree of $$\mathbb{R}^n$$ would be 1, while the degree of the surface of the unit sphere would be 2.

So here are some obvious questions I haven't been able to find answers to.

1. Are there any other non-trivial regular uniform metric spaces?

2. Can we find a regular metric space for every degree $$n$$? @user44191 Provided an affirmative answer to this question. A discrete set of $$n$$ points such that all distances between them are unique is a regular uniform metric space with degree $$n$$.

3. What can we say about the p-adic numbers $$\mathbb{Q}_p$$? They wouldn't be regular in the above sense. Consider for example the set of points that are distance 1 from 0. Then every other number will be a center of these numbers. We can make the same argument for the set of points that are distance $$\frac{1}{p}$$ from 0. If this was an equivalence relation then from the above to sets we should get that all points are in the same equivalence class. But looking back at the set of points with distance 1 from 0, we see that those points are not the center of it's own set as for example $$|5|_3 = 1$$ and $$|2|_3 = 1$$ but $$|5-2|_3 \neq 1$$. But I guess they are uniform.

4. What can we say about the surface of $$S^n$$? What can we say about a torus? A torus of genus 2? etc

5. Are there any equivalent definitions of uniformity or regularity?

6. Is there another condition other than regularity that would make this notion an equivalence class?

7. What can we say about the set of spaces on which this notion is an equivalence class?

• For question 2: consider the $n$-point discrete space where the distances are all distinct. If I understand correctly, this should be a regular metric space with degree $n$. – user44191 Nov 10 '19 at 4:46
• Definition of $C_r(x)$ features $d(x-y)$ where $x$ belongs to an abstract space $X$ while $y\in\Bbb R^n$. Something is not exactly right. – Wlod AA Nov 10 '19 at 4:49
• If $X$ is regular then it seems to me (am I wrong?) that arbitrary $x\in X$ can be equivalent only to itself and not to any other point. – Wlod AA Nov 10 '19 at 4:58
• @WlodAA I believe $y \in X$ is the intended definition; the sphere seems to be regular, with antipodes equivalent, if I'm not mistaken. – user44191 Nov 10 '19 at 5:00
• Considering $\mathbb{T}^2$ as a product of the circle with itself, we can define a "taxicab metric" on it by adding the distances on both projections; I think this gives us a regular metric space, with $(a, b) \sim (a + \frac{1}{2}, b + \frac{1}{2})$. Not sure though - you may want to check. – user44191 Nov 12 '19 at 23:01

Assume $$X$$ is a connected, regular space. I claim there are no triples $$x \sim y \sim z$$.

Proof:

Assume otherwise. First, note that if $$x \sim y$$, then $$C_{d(x, y)}(x) = \{y\}$$. This is because $$C_{d(x, y)}(x) = C_r(y)$$ for some $$r$$, but as $$y \in C_{d(x, y)}(x)$$, we have that $$y \in C_r(y)$$, so $$r = 0$$, and $$C_0(y) = \{y\}$$.

Assume $$x \sim y$$. Then let $$f: X \rightarrow \mathbb{R}_{\geq 0}^2, f(z) = (d(x, z), d(y, z))$$. The image includes $$(d(x, y), 0)$$ and $$(0, d(x, y))$$; the above note implies that the image cannot otherwise cross the lines $$(d(x, y), a)$$ or $$(b, d(x, y))$$. Then $$\{(a, b): a, b > d(x, y)\}$$ is open, as is $$\{(a, b): a < d(x, y) \, \text{OR} \, b < d(x, y)\}$$; their union contains all of $$\mathbb{R}_{\geq 0}^2$$ except points on those two lines (and the latter includes the two points above), so their preimages are open and contain all of $$X$$. Therefore, as $$X$$ is connected, one of the preimages must be empty - and $$(d(x, y), 0)$$ and $$(0, d(x, y))$$ are both in the latter, so the former must be empty. Therefore, for any $$z \in X$$, we must have that $$d(x, z) < d(x, y)$$ or $$d(y, z) < d(x, y)$$ (or both).

Now assume $$x \sim y \sim z$$. Then by the first note, we have that $$d(x, y) \neq d(x, z) \neq d(y, z)$$; one of the three distances must therefore be smallest, say (without loss of generality) $$d(x, y)$$. But then $$d(x, z) > d(x, y)$$ and $$d(y, z) > d(x, y)$$, contradicting the above paragraph.

Using the above, I'm going to rephrase the question slightly. Given two points $$x, y \in X$$, define $$f_{x, y}: X \rightarrow \mathbb{R}_{\geq 0}^2, z \mapsto (d(x, z), d(y, z))$$. Then $$x \sim y$$ iff $$f(X)$$ meets each horizontal and vertical line at most once, that is, if $$d(x, z) = d(x, w)$$, then $$d(y, z) = d(y, w)$$, and vice versa.

Using this definition, I claim that if $$x \sim y$$, then for any $$z \in X$$, we have that $$d(x, z) \leq d(x, y)$$. In other words, $$y$$ is the most distant point from $$x$$.

Assume otherwise. Note that by the above proof, we have that $$d(y, z) < d(x, y)$$. Then consider the disjoint open sets $$A, B \subseteq \mathbb{R}_{\geq 0}^2$$:

$$A := \{(a, b): a < d(x, y), b > d(y, z)\}$$

$$B := \{(a, b): a > d(x, y) \, OR \, b < d(y, z)\}$$

Clearly, $$f(x) = (0, d(x, y)) \in A, f(y) = (d(x, y), 0) \in B$$, and both are in $$f(X)$$. But $$f(X)$$ is connected, so there must be some point $$f(w)$$ in neither of these open sets. Therefore either $$d(x, w) = d(x, y)$$ and $$d(y, w) > d(y, z)$$, or $$d(x, w) < d(x, y), d(y, w) = d(y, z)$$. In the former case, $$d(y, w) = d(y, y) = 0$$, but we assumed $$w \neq y$$. In the latter case, $$d(x, w) = d(x, z)$$, but we assumed $$d(x, z) > d(x, y)$$, and $$d(x, y) > d(x, w)$$. So this is a contradiction, meaning our assumption was wrong, and if $$x \sim y$$, then $$y$$ is the most distant point from $$x$$.

If $$X$$ is a uniform, connected regular space of degree $$2$$, then this gives us a nice "antipode" map. I strongly suspect this map is continuous.

We can then define a function $$g: [0, d(x, y)] \rightarrow [0, d(x, y)]$$ such that if $$d(x, z) = c$$, then $$d(y, z) = g(c)$$. The proof above should be relatively easy to adapt to prove that $$g(c)$$ is increasing and continuous.

It's not hard to see that if $$X, Y$$ are regular with degree $$2$$ and the circle distances are additive (i.e. $$C_r(x) = C_{c - r}(y)$$ for some $$c$$), then $$X \times Y$$ is also regular with degree $$2$$ if the metric used is $$d((x_1, y_1), (x_2, y_2)) = d(x_1, x_2) + d(y_1, y_2)$$. This applies to the torus - in fact, it gives an uncountable number of metrics on the torus, as the metric on either factor circle can be "stretched" by any constant.

• I am really glad you are still working on this. It seems this was another case of Maschler spaces. I followed your argument "the above note implies that the image cannot otherwise cross the lines". I will try to understand more of it tomorrow. – Heraiwa Nov 11 '19 at 5:28
• @Heraiwa If you draw the lines in the plane, the argument should be pretty clear: the image can't cross a certain horizontal line or a certain vertical line, except at the axis - so as the image is connected, it can't have anything on the other side of both of those lines. The way I worded it was somewhat formal because I thought it would be clearer. – user44191 Nov 11 '19 at 6:07


text:
## Main Question or Discussion Point

How can i tell compuler select the correct type to whether to be able to 'infile' into int or float vars declared as private member in a class. sorry for my bad english.
Code:
template<bool TypeFlag, typename iType, typename fType>
struct Selector{
typedef iType Result;
};

template<typename iType, typename fType>
struct Selector<false, iType, fType>{
typedef fType Result;
};
try{
std::ifstream infile(filename);

//////////////What can I do here///////////

catch(...){
exit(EXIT_FAILURE);
}
}
Regards,

Related Computing and Technology News on Phys.org
I'm not sure I understand you. Do you want to have the computer read something from an infile and then decide whether it is an int or a float? In that case, it's impossible, unless you impose additional requirements (eg. int variables should be between 0 and 1000). Suppose you write something to an output file using something like:

int data_1 = 10;
float data_2 = 20;
outfile << data_1 << data_2;

The problem is that no type information is saved - you have to keep track of what went in through some other method.

dduardo
Staff Emeritus
Here's an idea: When you create the file print a character before the number to signify which type it is. For example

i123 d12.42

i = int
d = double

When you read back the file you'll know what type to store the variable in.

You can do the same thing with classes, for example:

classnamebject_name {variable1 , classnamebject_name {variable 1, variable 2} , variable 3}

When you convert variables or objects using this technique it is called serializing data. When you read back the data it is called unserializing.

thank you, zefram_c and dduardo

i have a class that i would use templates, the types will be, say, typename T and typename U, and in this class, i have some other functions among which there will be some that are 100% the same as each other in implementation but differ in types, say X and Y, that will be passed later from code in main function. i am not checking types of variables, just looking for a way to 'switch' from type to type and ways to use that 'switcher' in actuality. i posted code above only as an example. so, do you know how to deal with this problem ? i think it is trivial to you. but i am quite a dumb guy and really busy these day, i can't think much, i need help.
thanks if you or someone could be a little bit more definitive about the problem i am having.

Regards,
---really sorry for my bad english and misunderstandings if any.

kronecker, i don't understand what you want to do. What exactly is it that you want to "switch"? With your example above you could do something like:

Code:
// modify this at compile time, e.g. thru preprocessing based on platform
const bool a = true;

template <typename T>
void f() {cout << typeid(T).name() << endl; }

int main()
{
typedef Selector<a, int, float>::Result T;
f<T>();
}
But I'm not sure where you could use something like this. Or even if I'm close to what you want.. BTW, the most important thing to notice in the code above is that "a" MUST be computable at compile time.

i write simple classes for numerical computation. i was looking for ways to use that selector more effectively or maybe some other ways to handling type problems, my previous posts should already meantion, thank you very much for your help anyway.

What type problems?


text:
# Mersenne Primes and the Lucas–Lehmer Test

September 30, 2016 — John McGee, Applications Developer, Wolfram Technology Group

## Introduction

A Mersenne prime is a prime number of the form Mp = 2p – 1, where the exponent p must also be prime. These primes take their name from the French mathematician and religious scholar Marin Mersenne, who produced a list of primes of this form in the first half of the seventeenth century. It has been known since antiquity that the first four of these, M2 = 3, M3 = 7,
M5 = 31 and M7 = 127, are prime.

Mersenne claimed that 2p – 1 was prime for primes p ≤ 257 only for
p ∈ {2,3,5,7,13,17,19,31,67,127,257}. It is easy to verify where he was correct and where he was not, using the Wolfram Language function PrimeQ. PrimeQ uses modern prime testing methods that do not require finding a factor to prove a number to be composite.

It is possible that his claim that M67 is prime was a typographical error of M61. However, it is not hard to understand why primality testing was difficult in Mersenne’s time, since trial division was one of the few tools available. For example, for M257, the smallest factor is a 15-digit number and, even with modern factoring methods, it is not easy to find. The Wolfram Language function FactorInteger uses advanced methods that enable it to factor large integers.

## Euler and the Mersenne Factor Theorem

Some of the first advances in primality testing were accomplished by the great mathematician Leonhard Euler, who verified that M31 is prime sometime before 1772. He did this by showing that any prime divisor of M31 must be congruent to 1 or 62 (mod 248).

Such a relatively short list could be checked by trial division (by hand) in a reasonable amount of time in Euler’s day. His was an application of the Mersenne factor theorem, which states that if q is a divisor of Mp, then q ≡ 1 or –1 (mod 8), q ≡ 1 (mod p) and q = 2kp + 1 for some positive integer k. These facts greatly limit the possible divisors of Mp. The functions below use this theorem to provide a list of possible factors of Mp that are less than .

We use these functions to quickly find a factor of 241 – 1. Note that q is a factor of 2p – 1 if and only if 2p ≡ 1 (mod q). This enables the use of PowerMod, which provides very efficient modular exponentiation.

The following is a Mersenne number with 161,649 digits

## Lucas and Lehmer

The next major advance was the discovery by Édouard Lucas of a clever method to test the primality of numbers of this form. He used his method in 1876 to verify that M127, the largest Mersenne prime discovered before the age of computers, is prime. In the early twentieth century, after the understanding of binary arithmetic and algebra became widely known, Derek Henry Lehmer refined Lucas’ method. The resulting Lucas–Lehmer primality test provides an efficient method of testing if a number of this form is prime. It does this by using the modular equivalence

This means that k is congruent to the number represented by its lowest-order p bits plus the number represented by the remaining bits. This relation can be applied recursively until k < 2p – 1.

Consider the example that follows. Here we show that for
k = 1234567891. Note that , the lowest-order 23 bits, and , the remaining bits shifted to the lowest position.

The function below encodes this method to compute k (mod 2p – 1) using bit operations only (no division). Notice that 2n – 1 has the binary form 111 … 1112, all 1s and no 0s, so it also serves as a mask for the lower-order p bits of k.

The following function encodes the Lucas–Lehmer primality test (LLT). We define the sequence s0 = 4, s2i – 1 – 2; i > 0. Then Mp = 2p – 1 is prime iff sp – 1 ≡ 0(modMp).

Note: Experiments have shown that the runtime of these functions is dominated by the large integer arithmetic.

To efficiently test if 2p – 1 is prime, it is better to first check for small prime divisors and to perform other basic primality testing. We first use the Mersenne prime divisor theorem encoded in checkMPDivisors and then the Wolfram Language function PrimeQ. If neither of these complete in a short time, then we apply the Lucas–Lehmer test.

Here we present an extended version of PrimeQ that applies the Lucas–Lehmer test for large integers of the form 2p – 1.

## M13 to M20

The first Mersenne prime discovered by a computer running the Lucas–Lehmer test was M521, found by Raphael M. Robinson on January 30, 1952, using the early vacuum tube-based computer SWAC (Standards Western Automatic Computer). The Williams tube memory unit of this computer, holding 256 words of 37 bits each, is shown below.

The 20th Mersenne prime was discovered by Alexander Hurwitz in November of 1961 by running the Lucas–Lehmer test for about 50 minutes on an IBM 7090. We reproduce these early results below, using about 151 seconds of single-core computing time on a modern laptop.

One feature of the Wolfram Language that makes it suitable for this kind of work is its fast, large-integer arithmetic. This was a real challenge in the early days of computerized Mersenne prime searching. Researchers quickly adopted fast Fourier transform methods to convert the problem of multiplying two large integers, essentially a convolution of two lists of digits, into a simple element-by-element product of transformed digits. Fast integer multiplication is needed for the squaring step in the Lucas–Lehmer test. The Wolfram Language uses the latest platform-optimized algorithms to work with exact integer numbers with up to billions of digits. By way of example, we verify that the last of these, M4423, is indeed a Mersenne prime and show all of its digits.

## Perfect Numbers

There is an interesting connection between Mersenne primes and perfect numbers. A perfect number is a number that is equal to the sum of all of its divisors (other than the number itself). Euclid suspected, and Euler finally proved, that all even, perfect numbers have the form P = 2p – 1(2p – 1) = 2p – 1Mp. The Wolfram Language function PerfectNumberQ checks if a number is perfect. We demonstrate this property directly for M31.

## The 21st, 22nd and 23rd Mersenne Primes

We proceed to rediscover #21 = M9689, #22 = M9941 and #23 = M11213. These were all discovered by Donald B. Gillies running the LLT on an ILLIAC II during the spring 1963 (the article can be found here). We use nearly 6 minutes of elapsed time to test all of the numbers of the form 2p – 1 for primes 7,927 ≤ p ≤ 17,389.

## The 24th, 25th and 26th Mersenne Primes

We next extend the search to find #24 = M19937,
#25 = M21701 and #26 = M23209. The last of these was discovered in February of 1979 by Landon Curt Noll and Laura Nickel. They searched the range M21001 to M24499 using 6,000 CPU hours on a CDC Cyber 174 (that article can be found here). Our computations are becoming sufficiently intense to warrant the use of parallel processing. Since the tests of the candidate prime factors are independent, we can use ParallelMap to speed up the work. We check the range
17,393 ≤ p ≤ 27,449 in about three and a half minutes using 4 cores.

Notice how the specialized Lucas–Lehmer test is significantly faster than the more general function PrimeQ for these Mersenne primes.

## The 27th and 28th Mersenne Primes

We next test the range 27457 ≤ p ≤ 48611 to locate
#27 = M44497. This was discovered in April 1979 on a Cray-1 by Harry Nelson and his team. Our search of this range runs in about 15 minutes.

The next Mersenne prime is #28 = M86243. It was discovered in September of 1982 by David Slowinski, also on a Cray-1. The Cray-1 supercomputer weighed about 5 tons, consumed about 115 kilowatts of power and delivered 160 MFLOPS of computing performance. It was supplied with 1 million 64-bit words of memory (8 megabytes), and cost about $16 million in today’s dollars. A detail of its significant cooling system is shown below. By comparison, a Raspberry Pi weighs a few ounces, runs on 4 watts, delivers about 410 MFLOPS and is provided with 1 gigabyte of RAM, all for about$40, and it comes with Mathematica.

The number M86243 has 25,962 digits. In 1 hour and 14 minutes we were able to find this value (on my laptop, not on a Raspberry Pi) by testing over the range 48,619 ≤ p ≤ 87,533.

## The 29th Mersenne Prime

Since we are now using serious computer time, we also produce a timestamp for each run. We now check the range 87,557 ≤ p ≤ 110,597. In 1 hour and 44 minutes, this reveals #29 = M110503, first discovered on January 29, 1988, by Walker Colquitt and Luke Welsh running the LLT on an NEC DX-2 supercomputer (the article can be found here).

## The 30th and 31st Mersenne Primes

The next two Mersenne primes, M132049, the 30th, and M216091, the 31st, were actually discovered before #29, by the same team that discovered #28. They used a Cray X-MP to find #30 in September of 1983 and #31 in September 1985. We verify #30 by searching the range 110,603 ≤ p ≤ 139,901. It took nearly 4 hours and 8 minutes to check each Mp in this range.

## The Great Internet Mersenne Prime Search

The discovery of the 34th Mersenne prime, M1257787, in September 1996 ended the reign of the supercomputer in the search for Mersenne primes. The next 15 were found by volunteers of the Great Internet Mersenne Prime Search (GIMPS), which runs a variant of the Lucas–Lehmer test as a background process on personal computers. This large-scale distributed computing project currently achieves a performance equivalent to approximately 300 TFLOPS per second, harnessing the otherwise idle time of more than 1.3 million computers.

We verify the 34th Mersenne prime by directly using the Lucas–Lehmer test. We are reaching the limits of personal computer capability. Testing thousands of Mersenne numbers in this range would take many days. It is interesting to note that the Lucas–Lehmer test is often used as a stress test for the reliability of computer hardware and software, as even one arithmetic error among the billions of computations needed for testing one large prime will produce an incorrect conclusion, miss a true Mersenne prime or falsely report that a composite is prime. The fact that we have tested every Mp for primes between 2 and 139,901 is strong evidence for the reliability of large integer arithmetic and binary operations in Mathematica.

## Factoring Mersenne Numbers

As we have seen, the possible factors of numbers of the form
2p – 1 are limited by the Mersenne factor theorem. This has enabled an efficient computerized search for the factors of large integers of this form. The integer q divides 2p – 1 only if 2p ≡ 1 (mod q), which can be quickly checked with PowerMod. At the time of this writing, all Mersenne numbers up to 21201 – 1 have been fully factored. A Wolfram Demonstration that illustrates these can be found here. The GIMPS project has also led to the discovery of large factors of many Mersenne numbers as a byproduct of its primality testing. A recent paper that presents a modern approach to this problem, along with 17 new factorizations, can be found here. The largest number factorized was 21199 – 1; its smallest newly found factor has 76 decimal digits. Its smallest factor is actually 23. The total computer time for this effort was nearly 7,500 core years.

We can quickly find the first few factors of 21201 – 1 using the Wolfram Language function FactorInteger.

The Wolfram Language has cataloged of all of the Mersenne primes discovered to date, with ordering up to #44. Access to this information is provided by the functions MersennePrimeExponent and MersennePrimeExponentQ.

If you find this subject interesting, you can find more details at the following websites.

 Thanks John This is a very interesting blog. But besides trying to numerically check for primality, what is the status of attempts at proving that all Mersenne numbers are prime? How does it relate to the distribution of primes and the Riemann hypothesis about the zeros of the Zeta function? Regards Michael Posted by Michael    September 30, 2016 at 1:48 pm
 Michael, the Lucas-Lehmer test itself has been mathematically proven to return true, only for a Mersenne prime. The results produced by our version exactly match publish results produced by other researchers. As far as the density of Mersenne primes, it is not even known if this is an infinite set. See a discussion here: http://primes.utm.edu/mersenne/heuristic.html concerning conjectures about the density of these primes. Posted by John McGee    October 7, 2016 at 8:28 am


text:
# How much math is needed to understand Fourier Transform methods for option pricing?

I know of FDM and MC methods for option pricing, but have little to no experience with Fourier Methods.

I am intending to dive into the literature on using Fourier methods for option pricing, starting e.g. with the paper by Carr and Madan in 1999 I think (haven't read it yet, but have been recommended to start there).

My question is, is this particular literature demanding from a mathematical perspective? How much is required to properly understand and apply the methods that are currently being developed?

Should I read a book on Fourier Analysis before attempting to digest these finance-oriented papers?

My current mathematical background is mostly in real analysis and probability theory.

• At the end of the day FT is just a (complex) integral, a little complex analysis would help but if you are comfortable with integration theory, you should be able to digest the paper. Key insight is that density is the Fourier transform of the characteristic function. Also notice that one needs to think of the existence of the FT, so in the paper we need the dampening factor $\exp ( \alpha k )$ throughout to ensure the existence of the FT. – starovoitovs Jul 10 at 23:26
• Perhaps this can be helpful: quant.stackexchange.com/questions/9241/… – sets Jul 15 at 9:25

You do not need to know too much Fourier theory. Of course it helps, but it is not necessary. There are many different applications of Fourier methods e.g.

1. Carr Madan (1999): you damp the option price as a function of the log-strike price and compute the fourier transform of the entire option price
2. Bakshi and Madan (2000), Duffie, Pan and Singleton (2000): You Fourier transform the density and obtain a formula similar to the Black-Scholes solution. A decomposition of the call option price in Delta and exercise probability.
3. Lewis (2001): you Fourier transform the option payoff (these have typically polynomial growth, so you need a generalised Fourier transform which extends to the complex domain. This approach requires some complex analysis (residue theorem) and is in fact equivalent to Carr Madan (1999): Choosing a contour to integrate along or an optimal damping factor is the same question.
4. CONV and COS method: These are stable numerical algorithms which, for instance, allow very fast computation of European, Bermudan and (as a limit) American options.

In general, keep in mind, that Fourier methods do not apply to strongly path dependent options (asians, look-backs etc.) Furthermore, the main idea is always to replace the integral (expectation) which occurs from risk-neutral pricing by another integral which contains the characteristic function.

Regarding the maths, there are two things I believe are particular helpful.

1. Gil-Paelz inversion theorem What is option pricing about? Assuming a model for the distribution of $$S_T$$ and compute the (discounted) expectation of the payoff. This is an integral of payoff times density. Many processes and models (e.g. SVJ, VG, NIG, CGMY etc) have complicated densities but easy characteristic functions. Note that \begin{align*} \varphi_{\ln(S_T)}(u) = \mathbb{E}^\mathbb{Q}\left[e^{iu\ln(S_T)}\right] = \int_\mathbb{R} e^{iux} f_{\ln(S_T)}(x) \mathrm{d}x, \end{align*} where $$f_{\ln(S_T)}$$ is the risk-neutral density of $$\ln(S_T)$$. Thus, the characterisitc function of $$\ln(S_T)$$ is the Fourier transform of its density $$f_{\ln(S_T)}$$. Hence, $$\varphi_{\ln(S_T)}$$ captures the distribution of $$\ln(S_T)$$ and we can show that \begin{align*} F_{\ln(S_T)}(x) &= \frac{1}{2}+\frac{1}{2\pi} \int_0^\infty \frac{e^{iux}\varphi_{\ln(S_T)}(-u)-e^{-iux}\varphi_{\ln(S_T)}(u)}{iu}\mathrm{d}u. \end{align*} This will help a lot in deriving different option pricing formulae. Note that we typically consider $$\ln(S_T)$$ rather than $$S_T$$ since the characteristic function already involves somehow the exponential function.

2. Uncertainty principle. This is a concept from physics and states that if $$f$$ ''spreads out widely'', then its Fourier transform $$\hat{f}$$ is rather peaked and has a ''small'' support. This means that if you have an option with a short maturity, the density function will be peaked since there is no time for large movements. After all, a stock won't move a lot in a week or two. Thus, the chracterisitic function, $$\hat{f}$$, will spread out a lot and due to the oscillating behaviour, it may be challenging to integrate numerically.

Schmelzle wrote a nice survey paper but his work includes some errors in the section about the Lewis (2001) approach. If you're interested in the implementation of these models, have a look at Hirsa.

P.S. When I looked the first time at Fourier methods, my background was real analysis and probability theory. Just like yours. And I could follow it and understand it just fine. Don't you worry :)


text:
# Is there a name for this notion of “radius of compactness” in a metric space?

I was proving some result about Riemannian manifolds that led me to introduce the following definition:

Let $$M$$ be a metric space and $$x \in M$$. Define the "radius of compactness" $$RC(x)$$ to be the $$\sup$$ of $$R \geq 0$$ such that the closed ball $$\bar B(x, R)$$ is compact.

I like to think of it as a measure of how far a point is from the "boundary" of $$M$$.

It satisfies nice properties such as:

• $$|RC(y) - RC(x)| \leq d(x,y)$$
• If $$M$$ is locally compact and $$K \subset M$$ compact, then there exists $$\epsilon > 0$$ such that $$RC(x) \geq \epsilon$$ for all $$x \in K$$.

Is there a standard name for this $$RC$$?


text:
Casino On The Internet Betting System - Good Development System

If you discuss about the on line casino on the internet betting system, you will uncover there are several individuals who will discourage you. They will say that betting on the internet truly is not a great useful resource to make cash. But I will say that it is very simple to make from on the internet on line casino game titles, if you know the casino on-line betting techniques. Really money management information is what most of the gamblers deficiency. 우리카지노 ’s why some are currently bankrupt whereas some are enjoying an affluent lifestyle. Did any person of you listen to about “Positive Development System”, this is one particular of the extremely properly identified casino on the web betting technique.

You can say this is a logic that tells you the possibilities of profitable four occasions in a solitary row. At the starting or just at the preliminary phase the guess is of one unit, the 2nd wager is of 3 units, the third guess is of two units and the fourth bet is of 6 units. Therefore it is also known as the 1-three-2-six program.

I will illustrate this casino on the internet betting method in detail, to give you a distinct understanding. For occasion you location your initial wager of $10. The 2nd guess is intended to be$thirty - when you get the very first guess, your $ten gets extra up with the$20 already placed on the table. The complete arrives to $thirty. So the 2nd bet you area would be of$thirty. The grand total ahead of you perform the 3rd wager will be of $60 total (the$30 guess put by you in the 2nd guess blended collectively with the next guess winning presently put on the table). From the $60 you just take absent$forty and the 3rd wager is of $twenty. Your third wager will be of$20 and soon after winning the third wager you will win $forty. Now, for the fourth wager you will insert$20 a lot more to the whole $40 to make it a$sixty bet for the forth wager you spot.

Profitable the fourth wager you will be remaining with $one hundred twenty. This is the net earnings you make from this casino online betting method. To proceed the match you will yet again spot a guess of$ten and stick to the “Constructive Development Technique” after yet again. Following ending the forth wager, you begin in excess of again. Additionally, each and every time you unfastened a guess, begin once again with first $10 bet. Compute your loss at every single level. If you get rid of in the course of the: - first or very first bet, then your overall reduction is of$10.
- second bet, then your overall decline will be of $twenty (this since you have additional$10)
- third wager, then your gathered decline will also give you a overall profit of $twenty (this is since you have presently taken away$forty).
- fourth bet, then your loss and earnings is equal (just due to the fact you have taken absent $40 and additional$twenty).

The good facet of this on line casino on the web betting system is you just stake $twenty with a chance of creating a revenue that is 6 moments far more than the stake. In other words and phrases you are generating$a hundred and twenty by just jeopardizing the volume \$20.

It is really crucial to adhere to a casino on the internet betting system to engage in properly without going bankrupt.


text:
#### Change prompt

from octave:1> to >>

octave:1> PS1('>> ');
>>

#### Help command

>> help help
'help' is a function from the file C:\Octave\OCTAVE~1.0\mingw64\share\octave\5.2.0\m\help\help.m

-- help NAME
-- help --list
-- help .
-- help
Display the help text for NAME.

For example, the command 'help help' prints a short message
describing the 'help' command.

Given the single argument '--list', list all operators, keywords,
built-in functions, and loadable functions available in the current
session of Octave.

Given the single argument '.', list all operators available in the
current session of Octave.

If invoked without any arguments, 'help' displays instructions on
how to access help from the command line.

The help command can provide information about most operators, but
NAME must be enclosed by single or double quotes to prevent the
Octave interpreter from acting on NAME.  For example, 'help "+"'
displays help on the addition operator.

Additional help for built-in functions and operators is
available in the online version of the manual.  Use the command
'doc <topic>' to search the manual index.

Help and information about Octave is also available on the WWW
at https://www.octave.org and via the help@octave.org
mailing list.

# Basic Operations

>> 5+6
ans = 11

>> 3-2
ans = 1

>> 5*8
ans = 40

>> 1/2
ans = 0.50000

>> 2^6
ans = 64

>> 1 == 2
ans = 0

>> 1 ~= 2
ans = 1

>> 1 && 0
ans = 0

>> 1 || 0
ans = 1

>> xor(1,0)
ans = 1

## Variables

>> a =3;
>> a
a = 3

>> b = 'hi';
>> b
b = 'hi'

>> c = (3>=1);
>> c
c = 1

>> a = pi;
>> a
>> a = 3.1416

### Displaying Variables

>> a = pi;
>> disp(a)
3.1416

#### strings

>> a = pi;
>> disp(sprintf('2 decimals: %0.2f', a));
2 decimals: 3.14
>> disp(sprintf('6 decimals: %0.6f', a));
6 decimals: 3.141593

#### format

>> a = pi;
>> format long
>> a
a = 3.14159165358979
>> format short
>> a
a = 3.1416

## Matrices

#### 3×2 matrix

>> A = [1 2; 3 4; 5 6]
A =
1   2
3   4
5   6

or

>> A = [1 2;
> 3 4;
> 5 6]
A =
1   2
3   4
5   6

>> V = [1 2 3 ]
V =
1   2   3

#### 3×1 Vector

>> V = [1; 2; 3;]
V =
1
2
3

### Creating vector from 1 to 2 increasing next element by 0.1

>> V = 1:0.1:2
V =
1.0000    1.1000    1.2000    1.3000    1.4000    1.5000    1.6000    1.7000    1.8000    1.9000    2.0000

### Creating vector from 1 to 6 increasing next element by 1

>> V = 1:6
V =
1   2   3   4   5   6

>> ones(2,3)
ans =
1   1   1
1   1   1

### Matrix 2×3 with all values eqauls 2

>> C = 2*ones(2,3)
C =
2   2   2
2   2   2

>> w = ones(1,3)
w =
1   1   1

### Row vector with 0

>> w = zeros(1,3)
w =
0   0   0

### Random matrix (numbers between 0 and 1)

>> rand(1,3)
ans =
0.32784   0.49321   0.17449
>> rand(3,3)
ans =
0.60369   0.15891   0.75957
0.47148   0.11786   0.86370
0.46272   0.51486   0.71124

### Gaussian random variables

>> randn(1,3)
ans =
-0.67605  -0.54028   0.42763

### More complicated

>> w = -6 + sqrt(10)*randn(1,100)
w =
Columns 1 through 8:
-4.5865653   -8.8356479   -8.1333614   -2.3439273   -7.3621630   -9.0450181   -2.1489934   -4.2785162
Columns 9 through 16:
-10.7109066   -7.5124041   -6.6674847  -14.4887638   -0.4808289   -2.1354753   -7.8977581   -8.9314986
Columns 17 through 24:
-8.9893458   -1.5510572  -10.3796268  -10.2004881   -1.9258233  -10.7401729   -1.4114864  -11.9918739
Columns 25 through 32:
-1.4288672   -9.1350859   -4.0339111   -7.1057222   -5.5066692   -5.3943613   -9.9372854  -12.3273985
Columns 33 through 40:
-9.1448272   -3.5632437   -9.8825123   -3.6862963   -0.9314631   -7.1585643   -7.5249914   -6.7338681
Columns 41 through 48:
-4.0692604   -3.1736484   -7.1340077   -5.6141781   -7.1375287  -10.1554003   -7.3884758   -8.2993959
Columns 49 through 56:
-7.9971092   -3.7486347   -4.9565946  -10.1077235   -3.7186839   -7.6239998   -9.3567788   -7.0285528
Columns 57 through 64:
-5.7498064   -5.1235503    2.6385315   -7.0795552   -5.1845500   -2.7447164   -7.7785950    0.8816560
Columns 65 through 72:
-6.5335031    0.0066014   -3.9586199   -2.1445580   -3.4484948  -12.8979017   -1.5233025   -4.4102424
Columns 73 through 80:
0.9831482   -4.7574813   -2.7595514  -11.2178766   -5.3229671   -1.8543523   -3.7577500   -3.8678680
Columns 81 through 88:
-3.8637871   -7.2037168  -11.7877396   -7.1542280   -7.5999564    0.2206099   -6.7140798   -1.3302397
Columns 89 through 96:
-3.6058713   -5.5103303   -5.8951690   -1.6396029   -8.2389941   -4.5821969   -4.4162341   -6.9883700
Columns 97 through 100:
-3.4562746   -3.8389149   -0.8080181   -4.6449049

#### Histagram
r
>> hist(w)

>> hist(w,50)

>> eye(4)
ans =
Diagonal Matrix
1   0   0   0
0   1   0   0
0   0   1   0
0   0   0   1

# Moving Data Around

## Matrices size

>> A = [1 2; 3 4; 5 6]
A =
1   2
3   4
5   6

>> size(A)
ans =
3   2

>> size(A,1)
ans =  3

>> size(A,2)
ans =  2

V =
4   3   1   5
>> length(V)
ans =  4

>> length(A)
ans =  3

#### Show the current directory

>> pwd
ans = C:\Users\Gosia

#### Changing directory

>> cd 'D:\Projects'
>> pwd
ans = D:\Projects

#### display data from the file

>> data
data =
6.11010   17.59200
5.52770    9.13020
8.51860   13.66200
7.00320   11.85400
5.85980    6.82330
8.38290   11.88600
7.47640    4.34830
8.57810   12.00000
6.48620    6.59870
5.05460    3.81660
5.71070    3.25220
14.16400   15.50500
5.73400    3.15510
8.40840    7.22580
5.64070    0.71618
5.37940    3.51290
6.36540    5.30480
5.13010    0.56077
6.42960    3.65180
7.07080    5.38930
6.18910    3.13860

>> size(data)
ans =
97    2

#### check variables that Octave has in memory

>> who
Variables in the current scope:
A     V     a     ans   b     c     data  w

#### check variables with details that Octave has in memory

>> whos
Variables in the current scope:

Attr Name        Size                     Bytes  Class
==== ====        ====                     =====  =====
A           3x2                         48  double
V           1x4                         32  double
a           1x1                          8  double
ans         1x2                         16  double
b           1x2                          2  char
c           1x1                          1  logical
data       97x2                       1552  double
w           1x100                      800  double

Total is 310 elements using 2459 bytes

#### assign only selected data from one variable to another

>> m = data(1:10)
m =

6.1101   5.5277   8.5186   7.0032   5.8598   8.3829   7.4764   8.5781   6.4862   5.0546

>> whos
Variables in the current scope:

Attr Name        Size                     Bytes  Class
==== ====        ====                     =====  =====
A           3x2                         48  double
V           1x4                         32  double
a           1x1                          8  double
ans         1x2                         16  double
b           1x2                          2  char
c           1x1                          1  logical
data       97x2                       1552  double
m           1x10                        80  double
w           1x100                      800  double

Total is 320 elements using 2539 bytes

#### save a variable to a file (in binary format)

>> save hello.mat m;
>> ls
Volume in drive D is Data
Volume Serial Number is 9C15-10B4

Directory of D:\machine-learning-ex\test

[.]         [..]        data.txt    hello.mat
2 File(s)          1,691 bytes
2 Dir(s)  88,129,142,784 bytes free

#### save a variable to a file (in ascii format)

>> save hello.txt m -ascii

hello.txt

6.11010000e+00 5.52770000e+00 8.51860000e+00 7.00320000e+00 5.85980000e+00 8.38290000e+00 7.47640000e+00 8.57810000e+00 6.48620000e+00 5.05460000e+00

#### clear selected variable

>> clear data
>> whos
Variables in the current scope:

Attr Name        Size                     Bytes  Class
==== ====        ====                     =====  =====
A           3x2                         48  double
V           1x4                         32  double
a           1x1                          8  double
ans         1x2                         16  double
b           1x2                          2  char
c           1x1                          1  logical
w           1x100                      800  double

Total is 116 elements using 907 bytes

#### clear all variables

>> clear
>> whos
>>

>> whos
Variables in the current scope:

Attr Name        Size                     Bytes  Class
==== ====        ====                     =====  =====
m           1x10                        80  double

Total is 10 elements using 80 bytes

#### Get specific element from the variable (example: 3rd column and 2nd row)

>> A=[1 2; 3 4; 5 6]
A =

1   2
3   4
5   6

>> A(3,2)
ans =  6

>> A(2,:)
ans =
3   4

### Get everything in a second column

>> A(:,2)
ans =

2
4
6

>> A([1 3],:)
ans =

1   2
5   6

>> A([2 3],:)
ans =

3   4
5   6

>> A([1 2 3],:)
ans =

1   2
3   4
5   6

#### Change data in matrix

>> A
A =

1   2
3   4
5   6

>> A(:,2) = [10;11;12]
A =

1   10
3   11
5   12

#### Add (append) another column vector data to the right

>> A = [A, [100; 101; 102]];
>> A
A =

1    10   100
3    11   101
5    12   102

>> A(:)
ans =

1
3
5
10
11
12
100
101
102

#### Combine two matrices

>> A = [1 2; 3 4; 5 6];
>> B = [11 12; 13 14; 15 16];
>> C = [A B]
C =

1    2   11   12
3    4   13   14
5    6   15   16

or

>> A = [1 2; 3 4; 5 6];
>> B = [11 12; 13 14; 15 16];
>> C = [A; B]
C =

1    2
3    4
5    6
11   12
13   14
15   16

# Computing on Data

### Multiply

>> A = [1 2; 3 4; 5 6];
>> C = [1 1; 2 2];
>> A*C
ans =

5    5
11   11
17   17

### Element-wise multiplication

>> A = [1 2; 3 4; 5 6];
>> B = [11 12; 13 14; 15 16];
>> A .*B
ans =

11   24
39   56
75   96

### Element-wise squaring

>> A = [1 2; 3 4; 5 6];
>> A .^ 2
ans =

1    4
9   16
25   36

### Element-wise Reciprocal

>> V = [1; 2; 3];
>> 1 ./ V
ans =

1.00000
0.50000
0.33333

### Element-wise inverse of A (matrix)

>> A = [1 2; 3 4; 5 6];
>> 1 ./A
ans =

1.00000   0.50000
0.33333   0.25000
0.20000   0.16667

### Element-wise logarithm

>> V = [1; 2; 3];
>> log(V)
ans =

0.00000
0.69315
1.09861

### Exponentiation

>> V = [1; 2; 3];
>> exp(V)
ans =

2.7183
7.3891
20.0855

### Element-wise absolute value

>> V = [1; 2; 3];
>> abs(V)
ans =

1
2
3

>> abs([-1; -3; -10])
ans =

1
3
10

### negative

>> V = [1; 2; 3];
>> -V
ans =

-1
-2
-3

or

>> V = [1; 2; 3];
>> -1*V
ans =

-1
-2
-3

### Increment all values by one

>> V = [1; 2; 3];
>> V + ones(length(V),1)
ans =

2
3
4

or

>> V = [1; 2; 3];
>> V + 1
ans =

2
3
4

### Matrix transport

>> A = [1 2; 3 4; 5 6];
>> A'
ans =

1   3   5
2   4   6

### Max value

>> a = [1 15 2 0.5];
>> max(a)
ans =  15

### Max value and index of this max value

>> val = max(a)
val =  15
>> [val, ind] = max(a)
val =  15
ind =  2

### Column-wise maximum (on matrices)

>> A = [1 2; 3 4; 5 6];
>> max(A)
ans =

5   6

or

>> A = [1 2; 3 4; 5 6];
>> max(A,[],1)
ans =

5   6

### Row-wise maximum (on matrices)

>> A = [1 2; 3 4; 5 6];
>> max(A,[],2)
ans =

2
4
6

### The max element in whole matrix

>> A = [1 2; 3 4; 5 6];
>> max(max(A))
ans =  6

or turn into vector

>> A = [1 2; 3 4; 5 6];
>> max(A(:))
ans =  6

### Element-wise comparison (showing true/false)

>> a = [1 15 2 0.5];
>> a < 3
ans =

1  0  1  1

### Element-wise comparison (showing with element)

>> a = [1 15 2 0.5];
>> find(a < 3)
ans =

1   3   4

### Magic (N-by-N magic square)

>> A = magic(3)
A =

8   1   6
3   5   7
4   9   2
>> A = magic(9)
A =

47   58   69   80    1   12   23   34   45
57   68   79    9   11   22   33   44   46
67   78    8   10   21   32   43   54   56
77    7   18   20   31   42   53   55   66
6   17   19   30   41   52   63   65   76
16   27   29   40   51   62   64   75    5
26   28   39   50   61   72   74    4   15
36   38   49   60   71   73    3   14   25
37   48   59   70   81    2   13   24   35

#### sum columns

>> sum(A,1)
ans =

369   369   369   369   369   369   369   369   369

>> sum(A,2)
ans =

369
369
369
369
369
369
369
369
369

#### sum diagonal

>> eye(9)
ans =

Diagonal Matrix

1   0   0   0   0   0   0   0   0
0   1   0   0   0   0   0   0   0
0   0   1   0   0   0   0   0   0
0   0   0   1   0   0   0   0   0
0   0   0   0   1   0   0   0   0
0   0   0   0   0   1   0   0   0
0   0   0   0   0   0   1   0   0
0   0   0   0   0   0   0   1   0
0   0   0   0   0   0   0   0   1

>> A .* eye(9)
ans =

47    0    0    0    0    0    0    0    0
0   68    0    0    0    0    0    0    0
0    0    8    0    0    0    0    0    0
0    0    0   20    0    0    0    0    0
0    0    0    0   41    0    0    0    0
0    0    0    0    0   62    0    0    0
0    0    0    0    0    0   74    0    0
0    0    0    0    0    0    0   14    0
0    0    0    0    0    0    0    0   35

>> sum(sum(A.*eye(9)))
ans =  369

>> flipud(eye(9))
ans =

Permutation Matrix

0   0   0   0   0   0   0   0   1
0   0   0   0   0   0   0   1   0
0   0   0   0   0   0   1   0   0
0   0   0   0   0   1   0   0   0
0   0   0   0   1   0   0   0   0
0   0   0   1   0   0   0   0   0
0   0   1   0   0   0   0   0   0
0   1   0   0   0   0   0   0   0
1   0   0   0   0   0   0   0   0

>> sum(sum(A.*flipud(eye(9))))
ans =  369

### finding by row and column

>> A = magic(3);
>> [r,c] = find(A >= 7)
r =

1
3
2

c =

1
2
3

>> [row,column] = find(A >= 7)
row =

1
3
2

column =

1
2
3

### Pseudo-inverse

>> A = magic(3)
A =

8   1   6
3   5   7
4   9   2

>> pinv(A)
ans =

0.147222  -0.144444   0.063889
-0.061111   0.022222   0.105556
-0.019444   0.188889  -0.102778

>> a = [1 15 2 0.5];
>> sum(a)
ans =  18.500

### Multiply all elements

>> a = [1 15 2 0.5];
>> prod(a)
ans =  15

### Rounds down elements ( Return the largest integer not greater than X)

>> a = [1 15 2 0.5];
>> floor(a)
ans =

1   15    2    0

### Rounds up elements ( Return the smallest integer not less than X)

>> a = [1 15 2 0.5];
>> ceil(a)
ans =

1   15    2    1

# Plotting Data

Sample generated data for the plot (numbers from 0 to up 0.98) and sin function

>> t=[0:0.01:0.98];
>> y1 = sin(2*pi*4*t);

Generating plot, where horizontal axis = t and verical axis = y1:

>> plot(t,y1);

When run once again the plot with different argument, it will be replaced with new plot:

>> y2 = cos(2*pi*4*t);
>> plot(t,y2);

To have both (or more) plots on one designer it's needed to use the hold on; function:

>> plot(t,y1);
>> hold on;
>> plot(t,y2);

### Adding labels, legend and title

>> xlabel('time');
>> ylabel('value');
>> legend('sin','cos');
>> title('my plot');

### Saving the figure

>> cd D:\; print -dpng 'myPlot.png'

>> close;

### To open two figures in seperate windows:

>> figure(1); plot(t,y1);
>> figure(2); plot(t,y2);

### Dividing plot by 1x2 grid and access first element and next second element

>> subplot(1,2,1);
>> plot(t,y1);

>> clf;

### Visualize matrix

>> A = magic(5)
A =

17   24    1    8   15
23    5    7   14   16
4    6   13   20   22
10   12   19   21    3
11   18   25    2    9

>> imagesc(A)

#### To show with grey colour with colour bar

>> imagesc(A), colorbar, colormap gray;